C MESU      SOURCE    FANDEUR   12/01/24    21:15:15     7258           

C   fourni la mesure d'un maillage

      SUBROUTINE mesu

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

-INC CCOPTIO
-INC CCGEOME
-INC SMCOORD
-INC SMELEME

      SEGMENT /FER/ (NFI(ITT),MAI(IPP),ITOUR)
      CHARACTER*4 mcle(3)
      DATA mcle / 'LONG','SURF','VOLU' /

      ncle=IDIM
      icle=0
      CALL LIRMOT(mcle,ncle,icle,0)
      IF (IERR.NE.0) RETURN
      CALL LIROBJ('MAILLAGE',meleme,1,iOK)
      IF (IERR.NE.0) RETURN

      idimp1=IDIM+1
      SEGACT,meleme
      ipt1=meleme
      itp=0
      xmesu=0.
      nsous=lisous(/1)
      DO isou=1,MAX(nsous,1)
        IF (nsous.NE.0) THEN
          ipt1=lisous(isou)
          SEGACT,ipt1
        ENDIF
        jtyp=3
        IF (ksurf(ipt1.itypel).EQ.ipt1.itypel) jtyp=2
        IF (kdegre(ipt1.itypel).EQ.ipt1.itypel) jtyp=1
        IF (jtyp.EQ.1.AND.nsous.GT.1) CALL ERREUR(25)
        IF (itp.EQ.0) itp=jtyp
        IF (itp.NE.jtyp) CALL ERREUR(25)
        IF (icle.EQ.0) icle=itp
        IF (icle.EQ.1.AND.itp.NE.1) CALL ERREUR(21)
        IF (icle.EQ.2.AND.itp.EQ.3) CALL ERREUR(21)
        IF (icle.EQ.3.AND.itp.EQ.1) CALL ERREUR(21)
        IF (IERR.NE.0) RETURN
        GOTO (100,200,300),icle
C  calcul de la longueur d'une ligne
 100    IF (IIMPI.NE.0) WRITE(6,*) ' Mesure d''une ligne'
        nbnn=ipt1.num(/1)
        IF (IDIM.EQ.1) THEN
          DO iel=1,ipt1.num(/2)
            ip1=idimp1*(ipt1.num(1,iel)-1)
            ip2=idimp1*(ipt1.num(nbnn,iel)-1)
            xp12=XCOOR(ip2+1)-XCOOR(ip1+1)
            xmesu=xmesu+ABS(xp12)
          ENDDO
        ELSE IF (IDIM.EQ.2) THEN
          DO iel=1,ipt1.num(/2)
            ip1=idimp1*(ipt1.num(1,iel)-1)
            ip2=idimp1*(ipt1.num(nbnn,iel)-1)
            xp12=XCOOR(ip2+1)-XCOOR(ip1+1)
            yp12=XCOOR(ip2+2)-XCOOR(ip1+2)
            xmesu=xmesu+SQRT(xp12*xp12+yp12*yp12)
          ENDDO
        ELSE IF (IDIM.EQ.3) THEN
          DO iel=1,ipt1.num(/2)
            ip1=idimp1*(ipt1.num(1,iel)-1)
            ip2=idimp1*(ipt1.num(nbnn,iel)-1)
            xp12=XCOOR(ip2+1)-XCOOR(ip1+1)
            yp12=XCOOR(ip2+2)-XCOOR(ip1+2)
            zp12=XCOOR(ip2+3)-XCOOR(ip1+3)
            xmesu=xmesu+SQRT(xp12*xp12+yp12*yp12+zp12*zp12)
          ENDDO
        ENDIF
        GOTO 500
C  calcul de la surface d'une aire
 200    GOTO (220,240),jtyp
C  cas d'une ligne
 220    IF (IIMPI.NE.0) WRITE(6,*) ' Mesure de surface dans un contour'
        CALL AVTRSF(ipt1,fer,ipt5)
        SEGSUP,ipt5
        ib=0
        xb=0.
        yb=0.
        zb=0.
        DO it=1,itour
          DO im=mai(it)+1,mai(it+1)
            ib=ib+1
            ip=idimp1*(nfi(im)-1)
            xb=xb+XCOOR(ip+1)
            yb=yb+XCOOR(ip+2)
            IF (IDIM.EQ.3) zb=zb+XCOOR(ip+3)
          ENDDO
        ENDDO
        xb=xb/ib
        yb=yb/ib
        zb=zb/ib
        IF (IIMPI.NE.0) WRITE(6,*) ' Barycentre ',xb,yb,zb
        xs=0.D0
        ys=0.D0
        zs=0.D0
        zv1=0.D0
        zv2=0.D0
        DO it=1,itour
          DO im=mai(it)+1,mai(it+1)
            im1=im+1
            IF (im1.GT.mai(it+1)) im1=mai(it)+1
            ip1=idimp1*(nfi(im)-1)
            ip2=idimp1*(nfi(im1)-1)
            xv1=XCOOR(ip1+1)-xb
            yv1=XCOOR(ip1+2)-yb
            xv2=XCOOR(ip2+1)-xb
            yv2=XCOOR(ip2+2)-yb
            IF (IDIM.GE.3) THEN
              zv1=XCOOR(ip1+3)-zb
              zv2=XCOOR(ip2+3)-zb
            ENDIF
            xs=xs+(yv1*zv2-zv1*yv2)
            ys=ys+(zv1*xv2-xv1*zv2)
            zs=zs+(xv1*yv2-yv1*xv2)
          ENDDO
        ENDDO
        xmesu=xmesu+0.5D0*SQRT(xs*xs+ys*ys+zs*zs)
        SEGSUP,fer
        GOTO 500
C  cas d'une surface
 240    IF (IIMPI.NE.0) WRITE(6,*) ' Mesure d''une surface'
        iad=LTEL(2,ipt1.itypel)-1
        jad=LDEL(2,iad+1)-1
        ityp=LDEL(1,iad+1)
        npfac=KDFAC(1,ityp)
        idep=KDFAC(2,ityp)
        ifep=idep+3*(KDFAC(3,ityp)-1)
        xs=0.
        ys=0.
        zs=0.
        Cc=0.
        DO iel=1,ipt1.num(/2)
          DO i=idep,ifep,3
            iafa=LFAC(jad+KFAC(i))
            ibfa=LFAC(jad+KFAC(i+1))
            icfa=LFAC(jad+KFAC(i+2))
            ia=ipt1.num(iafa,iel)
            ib=ipt1.num(ibfa,iel)
            ic=ipt1.num(icfa,iel)
            IF (IIMPI.EQ.5) WRITE(6,*) ' Triangle ',ia,ib,ic
            ia=idimp1*(ia-1)
            ib=idimp1*(ib-1)
            ic=idimp1*(ic-1)
            xb=XCOOR(ib+1)
            yb=XCOOR(ib+2)
            xv1=XCOOR(ia+1)-xb
            yv1=XCOOR(ia+2)-yb
            xv2=XCOOR(ic+1)-xb
            yv2=XCOOR(ic+2)-yb
            IF (IDIM.GE.3) THEN
              zb=XCOOR(ib+3)
              zv1=XCOOR(ia+3)-zb
              zv2=XCOOR(ic+3)-zb
              xs=yv1*zv2-zv1*yv2
              ys=zv1*xv2-xv1*zv2
            ENDIF
            zs=xv1*yv2-yv1*xv2
            Cc=Cc+SQRT(xs*xs+ys*ys+zs*zs)
          ENDDO
        ENDDO
        xmesu=xmesu+0.5D0*Cc
        GOTO 500
C  calcul d'un volume
 300    GOTO (320,320,340),jtyp
C  cas du volume interieur a une surface
C  du fait du probleme d'orientation des elements on ne peut pas
C  traiter le probleme sous-maillage par sous maillage
 320    CALL MESUVO(meleme,xmesu)
        GOTO 11
C  cas d'un volume. On va utiliser le decoupage en triangle de la
C  surface des elements mais en le dedoublant de facon a etre sur
C  de correspondre d'un element au voisin
 340    nbfac=LTEL(1,ipt1.itypel)
        IF (nbfac.EQ.0) GOTO 500
        iad=LTEL(2,ipt1.itypel)-1
        kdeg=KDEGRE(ipt1.itypel)-1
        Cc=0.
        DO iel=1,ipt1.num(/2)
C centre de gravite du volume
          xb=0.
          yb=0.
          zb=0.
          nbnn=ipt1.num(/1)
          DO i=1,nbnn
            ipt=idimp1*(ipt1.num(i,iel)-1)
            xb=xb+XCOOR(ipt+1)
            yb=yb+XCOOR(ipt+2)
            zb=zb+XCOOR(ipt+3)
          ENDDO
          xb=xb/nbnn
          yb=yb/nbnn
          zb=zb/nbnn
          DO ifac=1,nbfac
            v=0.
            ityp=LDEL(1,iad+ifac)
            npfac=KDFAC(1,ityp)
            jad=LDEL(2,iad+ifac)-1
            idep=KDFAC(2,ityp)
            ifep=idep+3*(KDFAC(3,ityp)-1)
            DO i=idep,ifep,3
              iafa=LFAC(jad+KFAC(i))
              ibfa=LFAC(jad+KFAC(i+1))
              icfa=LFAC(jad+KFAC(i+2))
              ia=ipt1.num(iafa,iel)
              ib=ipt1.num(ibfa,iel)
              ic=ipt1.num(icfa,iel)
              v=v+VOLTRI(xb,yb,zb,ia,ib,ic)
            ENDDO
            DO i=idep,ifep,3
              ifa=KFAC(i)+kdeg
              IF (ifa.GT.npfac) ifa=ifa-npfac
              ifb=KFAC(i+1)+kdeg
              IF (ifb.GT.npfac) ifb=ifb-npfac
              ifc=KFAC(i+2)+kdeg
              IF (ifc.GT.npfac) ifc=ifc-npfac
              iafa=LFAC(jad+ifa)
              ibfa=LFAC(jad+ifb)
              icfa=LFAC(jad+ifc)
              ia=ipt1.num(iafa,iel)
              ib=ipt1.num(ibfa,iel)
              ic=ipt1.num(icfa,iel)
              v=v+VOLTRI(xb,yb,zb,ia,ib,ic)
            ENDDO
            IF (IIMPI.EQ.1) WRITE(6,*) 'v-2 ',v
            Cc=Cc+ABS(v)
          ENDDO
        ENDDO
        xmesu=xmesu+Cc/12.D0
C*      GOTO 500
 500    SEGDES,ipt1
      ENDDO

 11   CONTINUE
      CALL ECRREE(xmesu)
      SEGDES,meleme

      RETURN
      END