Télécharger la plaquette Cast3M


fr en

Télécharger mesuvo.eso

Retour à la liste

Numérotation des lignes : 

mesuvo

  1. C MESUVO    SOURCE    PV        20/03/30    21:21:09     10567          
  2. C   mesure du volume contenu dans une surface fermee
  3. C
  4.       subroutine mesuvo(meleme,xmesu)
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8(A-H,O-Z)
  7.  
  8. -INC PPARAM
  9. -INC CCOPTIO
  10. -INC SMELEME
  11. -INC CCGEOME
  12. -INC SMCOORD
  13. C Segments locaux
  14. C icpr tableau de compression des noeuds concernes
  15. C kon (numero point compresse par icpr,rang d'element,1) => pointeur su meleme
  16. C kon (numero point compresse par icpr,rang d'element,2) => numero element
  17. C ifut (rang d'element,1) => pointeur sur meleme
  18. C ifut (rang d'element,2) => element traite
  19. C ifut (rang d'element,3) => orientation a utiliser
  20. C
  21.       segment kon(ite,nbk,3)
  22.       segment ikon(ite)
  23.       segment icpr(nbpts)
  24.       segment ifut(nbel,3)
  25.       xmesu=0.d0
  26.       nbel=0
  27.       segini icpr
  28.       ite=0
  29.       segact meleme*mod
  30.       ipt1=meleme
  31.       do 10 isous=1,max(1,lisous(/1))
  32.          if (lisous(/1).ne.0) then
  33.             ipt1=lisous(isous)
  34.             segact ipt1*mod
  35.             if (ipt1.itypel.ne.ksurf(ipt1.itypel)) call erreur(25)
  36.             if (ierr.ne.0) return
  37.          endif
  38. * cas des faces quaf tri7/qua9 : pas de prise en compte du noeud central
  39.          if (ipt1.itypel.eq.7.or.ipt1.itypel.eq.11) then 
  40.             idec1=1
  41.          else
  42.             idec1=0
  43.          endif
  44.          nbel=nbel+ipt1.num(/2)
  45.          do 15 ipt=1,ipt1.num(/1)-idec1
  46.             do 151 iel=1,ipt1.num(/2)
  47.                if (icpr(ipt1.num(ipt,iel)).ne.0) goto 151
  48.                ite=ite+1
  49.                icpr(ipt1.num(ipt,iel))=ite
  50.  151        continue
  51.  15      continue
  52.  10   continue
  53.       nbk=8
  54. C on remplit kon qui donne les elements touchant un noeud (compacte)
  55.       segini kon,ikon
  56.       do 20 isous=1,max(1,lisous(/1))
  57.          if (lisous(/1).ne.0) ipt1=lisous(isous)
  58. * cas des faces quaf tri7/qua9 : pas de prise en compte du noeud central
  59.          if (ipt1.itypel.eq.7.or.ipt1.itypel.eq.11) then 
  60.             idec1=1
  61.          else
  62.             idec1=0
  63.          endif
  64.          do 25 ipt=1,ipt1.num(/1)-idec1
  65.             do 251 iel=1,ipt1.num(/2)
  66.                ip=icpr(ipt1.num(ipt,iel))
  67.                ikon(ip)=ikon(ip)+1
  68.                if (ikon(ip).gt.nbk) then
  69.                   nbk=nbk+1
  70.                   segadj kon
  71.                endif
  72.                kon(ip,ikon(ip),1)=ipt1
  73.                kon(ip,ikon(ip),2)=iel
  74.  251        continue
  75.  25      continue
  76.  20   continue
  77. C on commence a travailler
  78.       segini ifut
  79.       do 30 i=1,nbel
  80.          ifut(i,3)=1
  81.   30  continue
  82.       ilcou=0
  83.       iltop=1
  84.       ifut(1,1)=ipt1
  85.       ifut(1,2)=1
  86.       ifut(1,3)=1
  87.       ipt1.num(1,1)=-ipt1.num(1,1)
  88.  100  continue
  89.       ilcou=ilcou+1
  90.       if (ilcou.gt.iltop) goto 500
  91.       ipt2=ifut(ilcou,1)
  92.       iel=ifut(ilcou,2)
  93.       isens=ifut(ilcou,3)
  94. * cas des faces quaf tri7/qua9 : pas de prise en compte du noeud central
  95.       if (ipt2.itypel.eq.7.or.ipt2.itypel.eq.11) then 
  96.          idec2=1
  97.       else
  98.          idec2=0
  99.       endif
  100.  
  101. C  on cherche les elements contactant iel et on les oriente correctement
  102. C  on note qu'ils sont traites en mettant en neg le premier noeud
  103.       do 110 i=1,ipt2.num(/1)-idec2
  104.          ip1=abs(ipt2.num(i,iel))
  105.          i1=i+1
  106.          if (i1.gt.ipt2.num(/1)-idec2) i1=1
  107.          i0=i-1
  108.          if (i0.lt.1) i0=ipt2.num(/1)-idec2
  109.          if (isens.eq.1) ip2=abs(ipt2.num(i1,iel))
  110.          if (isens.eq.-1) ip2=abs(ipt2.num(i0,iel))
  111.          do 120 j=1,ikon(icpr(ip1))
  112.             ipt3=kon(icpr(ip1),j,1)
  113.             jel=kon(icpr(ip1),j,2)
  114.             jsens=isign(1,ipt3.num(2,jel))
  115. *        write (6,*) ' isens,jsens,iel,jel ',isens,jsens,iel,jel,ip1,ip2
  116.             if (ipt3.eq.ipt2.and.jel.eq.iel) goto 120
  117. * cas des faces quaf tri7/qua9 : pas de prise en compte du noeud central
  118.             if (ipt3.itypel.eq.7.or.ipt3.itypel.eq.11) then 
  119.                idec3=1
  120.             else
  121.                idec3=0
  122.             endif
  123.             do 130 in=1,ipt3.num(/1)-idec3
  124.                if (abs(ipt3.num(in,jel)).ne.ip1) goto 130
  125.                in1=in+1
  126.                if (in1.gt.ipt3.num(/1)-idec3) in1=1
  127.                in0=in-1
  128.                if (in0.lt.1) in0=ipt3.num(/1)-idec3
  129.                if ((abs(ipt3.num(in1,jel)).eq.ip2.and.jsens.eq.1).or.
  130.      *              (abs(ipt3.num(in0,jel)).eq.ip2.and.jsens.eq.-1))
  131.      $              then
  132. C  Il faut reorienter l'element
  133.                   if (iimpi.eq.5) write(6,*) ' element ',ipt3,jel
  134.      $                 ,' change de sens'
  135.                   if (ipt3.num(1,jel).lt.0) goto 999
  136.                   if (ierr.ne.0) return
  137.                   ipt3.num(1,jel)=-abs(ipt3.num(1,jel))
  138.                   iltop=iltop+1
  139.                   ifut(iltop,1)=ipt3
  140.                   ifut(iltop,2)=jel
  141.                   ifut(iltop,3)=-jsens
  142.                   ipt3.num(2,jel)=-abs(ipt3.num(2,jel))
  143.                   goto 110
  144.                elseif ((abs(ipt3.num(in1,jel)).eq.ip2.and.jsens.eq.-1)
  145.      *                 .or.(abs(ipt3.num(in0
  146.      $                 ,jel)).eq.ip2.and.jsens.eq.1)) then
  147. C  Element OK
  148.                   if (iimpi.eq.5) write(6,*) ' element ',ipt3,jel,
  149.      $                 ' bon sens '
  150.                   if (ipt3.num(1,jel).lt.0) goto 110
  151.                   ipt3.num(1,jel)=-abs(ipt3.num(1,jel))
  152.                   iltop=iltop+1
  153.                   ifut(iltop,1)=ipt3
  154.                   ifut(iltop,2)=jel
  155.                   ifut(iltop,3)=jsens
  156.                   goto 110
  157.                else
  158.                endif
  159.  130        continue
  160.  120     continue
  161.          goto 999
  162.  110  continue
  163.       goto 100
  164.  500  continue
  165. C  il ne reste plus qu'a calculer le volume et remettre les elements corrects
  166.       if (iimpi.eq.5) write (6,*) ' iltop nbel ',iltop,nbel
  167.       if (iltop.ne.nbel) goto 999
  168.       if (ierr.ne.0) return
  169.       do 600 i=1,nbel
  170.          ipt1=ifut(i,1)
  171.          iel=ifut(i,2)
  172.          isens=ifut(i,3)
  173.          IAD=LTEL(2,IPT1.ITYPEL)-1
  174.          kdeg=kdegre(IPT1.ITYPEL)-1
  175.          ITYP=LDEL(1,IAD+1)
  176.          NPFAC=KDFAC(1,ITYP)
  177.          JAD=LDEL(2,IAD+1)-1
  178.          IDEP=KDFAC(2,ITYP)
  179.          IFEP=IDEP+3*(KDFAC(3,ITYP)-1)
  180.          dxmesu=0.d0
  181.          DO 610 ITRIAN=IDEP,IFEP,3
  182.             IAFA=LFAC(JAD+KFAC(ITRIAN))
  183.             IBFA=LFAC(JAD+KFAC(ITRIAN+1))
  184.             ICFA=LFAC(JAD+KFAC(ITRIAN+2))
  185.             IA=abs(IPT1.NUM(IAFA,IEL))
  186.             IB=abs(IPT1.NUM(IBFA,IEL))
  187.             IC=abs(IPT1.NUM(ICFA,IEL))
  188. *         if (iimpi.eq.5) write (6,*) ' triangle ',ia,ib,ic
  189.             dxmesu=dxmesu+isens*voltri(0.d0,0.d0,0.d0,ia,ib,ic)/12
  190.  610     continue
  191. C     SG 20160712 : comme dans mesu.eso, pour les faces TRI7, QUA9 (ityp
  192. C     =7,8), il n'y a pas besoin de dedoubler car le decoupage de la
  193. C     face est symetrique (barycentrique) 
  194.          IF (ITYP.NE.7.AND.ITYP.NE.8) THEN
  195.             DO 611 ITRIAN=IDEP,IFEP,3
  196.                ifa=KFAC(ITRIAN)+kdeg
  197.                if (ifa.gt.npfac) ifa=ifa-npfac
  198.                ifb=KFAC(ITRIAN+1)+kdeg
  199.                if (ifb.gt.npfac) ifb=ifb-npfac
  200.                ifc=KFAC(ITRIAN+2)+kdeg
  201.                if (ifc.gt.npfac) ifc=ifc-npfac
  202.                IAFA=LFAC(JAD+ifa)
  203.                IBFA=LFAC(JAD+ifb)
  204.                ICFA=LFAC(JAD+ifc)
  205.                IA=abs(IPT1.NUM(IAFA,IEL))
  206.                IB=abs(IPT1.NUM(IBFA,IEL))
  207.                IC=abs(IPT1.NUM(ICFA,IEL))
  208.                dxmesu=dxmesu+isens*voltri(0.d0,0.d0,0.d0,ia,ib,ic)/12
  209.  611        continue
  210.          ELSE
  211.             dxmesu=dxmesu*2
  212.          ENDIF
  213.          xmesu=xmesu+dxmesu
  214.  600  continue
  215.       goto 998
  216.  999  continue
  217.       call erreur(127)
  218.  998  continue
  219. C  maintenant on desactive les segments et on remet les elements correctement
  220.       ipt1=meleme
  221.       do 620 isous=1,max(1,lisous(/1))
  222.          if (lisous(/1).ne.0) ipt1=lisous(isous)
  223.          do 630 iel=1,ipt1.num(/2)
  224.             ipt1.num(1,iel)=abs(ipt1.num(1,iel))
  225.             ipt1.num(2,iel)=abs(ipt1.num(2,iel))
  226.  630     continue
  227.          segdes ipt1
  228.  620  continue
  229.       segdes meleme
  230.       segsup icpr,kon,ikon,ifut
  231.       xmesu=abs(xmesu)
  232.       end
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales