C MESUVO    SOURCE    CHAT      07/10/29    21:24:21     5940           
C   mesure du volume contenu dans une surface fermee
C
      subroutine mesuvo(meleme,xmesu)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
-INC CCOPTIO
-INC SMELEME
-INC CCGEOME
-INC SMCOORD
C Segments locaux 
C icpr tableau de compression des noeuds concernes 
C kon (numero point compresse par icpr,rang d'element,1) => pointeur su meleme 
C kon (numero point compresse par icpr,rang d'element,2) => numero element 
C ifut (rang d'element,1) => pointeur sur meleme
C ifut (rang d'element,2) => element traite
C ifut (rang d'element,3) => orientation a utiliser
C
      segment kon(ite,nbk,3)
      segment ikon(ite)
      segment icpr(xcoor(/1)/(idim+1))
      segment ifut(nbel,3)
      xmesu=0.
      nbel=0
      segini icpr
      ite=0
      segact meleme*mod
      ipt1=meleme
      do 10 isous=1,max(1,lisous(/1))
       if (lisous(/1).ne.0) then
        ipt1=lisous(isous)
        segact ipt1*mod
        if (ipt1.itypel.ne.ksurf(ipt1.itypel)) call erreur(25)
        if (ierr.ne.0) return
       endif
       nbel=nbel+ipt1.num(/2)
      do 15 ipt=1,ipt1.num(/1)
      do 15 iel=1,ipt1.num(/2)
       if (icpr(ipt1.num(ipt,iel)).ne.0) goto 15
       ite=ite+1
       icpr(ipt1.num(ipt,iel))=ite
  15  continue
  10  continue
      nbk=8
C on remplit kon qui donne les elements touchant un noeud (compacte)
      segini kon,ikon
      do 20 isous=1,max(1,lisous(/1))
       if (lisous(/1).ne.0) ipt1=lisous(isous)
      do 25 ipt=1,ipt1.num(/1)
      do 25 iel=1,ipt1.num(/2)
       ip=icpr(ipt1.num(ipt,iel))
       ikon(ip)=ikon(ip)+1
       if (ikon(ip).gt.nbk) then
        nbk=nbk+1
        segadj kon
       endif
       kon(ip,ikon(ip),1)=ipt1
       kon(ip,ikon(ip),2)=iel
  25  continue
  20  continue
C on commence a travailler
      segini ifut
      do 30 i=1,nbel
       ifut(i,3)=1
  30  continue
      ilcou=0
      iltop=1
      ifut(1,1)=ipt1
      ifut(1,2)=1
      ifut(1,3)=1
      ipt1.num(1,1)=-ipt1.num(1,1)
 100  continue
      ilcou=ilcou+1
      if (ilcou.gt.iltop) goto 500
      ipt2=ifut(ilcou,1)
      iel=ifut(ilcou,2)
      isens=ifut(ilcou,3)
C  on cherche les elements contactant iel et on les oriente correctement
C  on note qu'ils sont traites en mettant en neg le premier noeud
      do 110 i=1,ipt2.num(/1)
       ip1=abs(ipt2.num(i,iel))
       i1=i+1
       if (i1.gt.ipt2.num(/1)) i1=1
       i0=i-1
       if (i0.lt.1) i0=ipt2.num(/1)
       if (isens.eq.1) ip2=abs(ipt2.num(i1,iel))
       if (isens.eq.-1) ip2=abs(ipt2.num(i0,iel))
       do 120 j=1,ikon(icpr(ip1))
        ipt3=kon(icpr(ip1),j,1)
        jel=kon(icpr(ip1),j,2)
        jsens=isign(1,ipt3.num(2,jel))
*        write (6,*) ' isens,jsens,iel,jel ',isens,jsens,iel,jel,ip1,ip2
        if (ipt3.eq.ipt2.and.jel.eq.iel) goto 120
         do 130 in=1,ipt3.num(/1)
          if (abs(ipt3.num(in,jel)).ne.ip1) goto 130
          in1=in+1
          if (in1.gt.ipt3.num(/1)) in1=1
          in0=in-1
          if (in0.lt.1) in0=ipt3.num(/1)
          if ((abs(ipt3.num(in1,jel)).eq.ip2.and.jsens.eq.1).or.
     *     (abs(ipt3.num(in0,jel)).eq.ip2.and.jsens.eq.-1)) then
C  Il faut reorienter l'element
      if (iimpi.eq.5) write(6,*) ' element ',ipt3,jel,' change de sens'
           if (ipt3.num(1,jel).lt.0) goto 999
           if (ierr.ne.0) return
           ipt3.num(1,jel)=-abs(ipt3.num(1,jel))
           iltop=iltop+1
           ifut(iltop,1)=ipt3
           ifut(iltop,2)=jel
           ifut(iltop,3)=-jsens
           ipt3.num(2,jel)=-abs(ipt3.num(2,jel))
           goto 110
          elseif ((abs(ipt3.num(in1,jel)).eq.ip2.and.jsens.eq.-1)
     *   .or.(abs(ipt3.num(in0,jel)).eq.ip2.and.jsens.eq.1)) then
C  Element OK     
      if (iimpi.eq.5) write(6,*) ' element ',ipt3,jel, ' bon sens '
           if (ipt3.num(1,jel).lt.0) goto 110
           ipt3.num(1,jel)=-abs(ipt3.num(1,jel))
           iltop=iltop+1
           ifut(iltop,1)=ipt3
           ifut(iltop,2)=jel
           ifut(iltop,3)=jsens
           goto 110
           else
          endif
 130     continue
 120   continue
         goto 999
 110   continue
       goto 100
 500   continue
C  il ne reste plus qu'a calculer le volume et remettre les elements corrects
       if (iimpi.eq.5) write (6,*) ' iltop nbel ',iltop,nbel
       if (iltop.ne.nbel) goto 999
       if (ierr.ne.0) return
       do 600 i=1,nbel
        ipt1=ifut(i,1)
        iel=ifut(i,2)
        isens=ifut(i,3)
        IAD=LTEL(2,IPT1.ITYPEL)-1
        kdeg=kdegre(IPT1.ITYPEL)-1
        ITYP=LDEL(1,IAD+1)
        NPFAC=KDFAC(1,ITYP)
        JAD=LDEL(2,IAD+1)-1
        IDEP=KDFAC(2,ITYP)
        IFEP=IDEP+3*(KDFAC(3,ITYP)-1)
        DO 610 ITRIAN=IDEP,IFEP,3
         IAFA=LFAC(JAD+KFAC(ITRIAN))
         IBFA=LFAC(JAD+KFAC(ITRIAN+1))
         ICFA=LFAC(JAD+KFAC(ITRIAN+2))
         IA=abs(IPT1.NUM(IAFA,IEL))
         IB=abs(IPT1.NUM(IBFA,IEL))
         IC=abs(IPT1.NUM(ICFA,IEL))
*         if (iimpi.eq.5) write (6,*) ' triangle ',ia,ib,ic
           xmesu=xmesu+isens*voltri(0.d0,0.d0,0.d0,ia,ib,ic)/12
 610    continue
        DO 611 ITRIAN=IDEP,IFEP,3
         ifa=KFAC(ITRIAN)+kdeg
         if (ifa.gt.npfac) ifa=ifa-npfac
         ifb=KFAC(ITRIAN+1)+kdeg
         if (ifb.gt.npfac) ifb=ifb-npfac
         ifc=KFAC(ITRIAN+2)+kdeg
         if (ifc.gt.npfac) ifc=ifc-npfac
         IAFA=LFAC(JAD+ifa)
         IBFA=LFAC(JAD+ifb)
         ICFA=LFAC(JAD+ifc)
         IA=abs(IPT1.NUM(IAFA,IEL))
         IB=abs(IPT1.NUM(IBFA,IEL))
         IC=abs(IPT1.NUM(ICFA,IEL))
           xmesu=xmesu+isens*voltri(0.d0,0.d0,0.d0,ia,ib,ic)/12
611    continue
 600   continue
       goto 998
 999   continue
        call erreur(127)
 998   continue
C  maintenant on desactive les segments et on remet les elements correctement
       ipt1=meleme
       do 620 isous=1,max(1,lisous(/1))
        if (lisous(/1).ne.0) ipt1=lisous(isous)
       do 630 iel=1,ipt1.num(/2)
         ipt1.num(1,iel)=abs(ipt1.num(1,iel))
         ipt1.num(2,iel)=abs(ipt1.num(2,iel))
 630   continue
        segdes ipt1
 620   continue
       segdes meleme
       segsup icpr,kon,ikon,ifut
       xmesu=abs(xmesu)
       end