Télécharger la plaquette Cast3M


fr en

Télécharger vomsi2.eso

Retour à la liste

Numérotation des lignes : 

vomsi2

  1. C VOMSI2    SOURCE    GOUNAND   21/04/06    21:15:44     10940          
  2.       SUBROUTINE VOMSI2(MELEME,IELDEB,IELFIN,IMET,IPVIRT,
  3.      $        XVOL,XVOLV)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5.       IMPLICIT INTEGER (I-N)
  6. C***********************************************************************
  7. C NOM         : VOMSI2
  8. C DESCRIPTION : VOlume d'un Maillage de SIMplexes
  9. C               MELEME est supposé actif.
  10. C
  11. C     Comme VOMSIM mais on précise les numéros d'éléments de début
  12. C     et de fin à prendre en compte.
  13. C
  14. C LANGAGE     : ESOPE
  15. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  16. C               mél : gounand@semt2.smts.cea.fr
  17. C***********************************************************************
  18. C VERSION    : v1, 03/11/2017, version initiale
  19. C HISTORIQUE : v1, 03/11/2017, création
  20. C HISTORIQUE :
  21. C HISTORIQUE :
  22. C***********************************************************************
  23. -INC CCGEOME
  24. -INC PPARAM
  25. -INC CCOPTIO
  26. -INC SMCOORD
  27. -INC SMELEME
  28. -INC SMLREEL
  29.       logical lvirt
  30. * Statement functions
  31. *      DISTA(A,B,C,D)=SQRT((A-C)*(A-C)+(B-D)*(B-D))
  32. *      DISTB(A,B,C,D,E,F)=SQRT((A-D)*(A-D)+(B-E)*(B-E)+(C-F)*(C-F))
  33. *      DISTA(A,B)=SQRT(A*A+B*B)
  34. *      DISTB(A,B,C)=SQRT(A*A+B*B+C*C)
  35. *      DETTRI(A11,A12,A21,A22)=A11*A22-A21*A12
  36.       DETTET(A11,A12,A13,A21,A22,A23,A31,A32,A33)=
  37.      &  A11*(A22*A33-A23*A32)
  38.      &       +A12*(A23*A31-A21*A33)
  39.      &       +A13*(A21*A32-A22*A31)
  40.  
  41. *     
  42. * Executable statements
  43. *
  44.       XVOL=0.D0
  45.       XVOLV=0.D0
  46. *Maillage vide
  47. *      SEGACT MELEME
  48.       IF (LISOUS(/1).EQ.0.AND.NUM(/2).EQ.0) THEN
  49.          WRITE(IOIMP,*) 'ITLOC coucou vomsi2'
  50.          write(ioimp,*) 'LISOU1=',LISOUS(/1)
  51.          write(ioimp,*) 'NUM1=',NUM(/1)
  52.          write(ioimp,*) 'NUM2=',NUM(/2)
  53.          write(ioimp,*) 'MELEME=',MELEME
  54.          CALL ECMAI1(MELEME,0)
  55.          CALL ERREUR(5)
  56.          RETURN
  57.       ENDIF
  58.       IDIMP1=IDIM+1
  59.       NBNN=NUM(/1)
  60.       IF (NBNN.NE.IDIMP1.OR.LISOUS(/1).NE.0) THEN
  61.          WRITE(IOIMP,*) 'NBNN=',NBNN
  62.          WRITE(IOIMP,*) 'IDIMP1=',IDIMP1
  63.          WRITE(IOIMP,*) 'NBSOUS=',LISOUS(/1)
  64.          WRITE(IOIMP,*) 'ITLOC'
  65.          CALL ECMAI1(MELEME,0)
  66.          CALL ERREUR(5)
  67.          RETURN
  68.       ENDIF
  69.       IF (.NOT.(IELDEB.GE.1.AND.IELFIN.GE.IELDEB.AND.NUM(
  70.      $     /2).GE.IELFIN))  THEN
  71.          WRITE(IOIMP,*) 'ITLOC coucou vomsi2'
  72.          write(ioimp,*) 'IELDEB=',IELDEB
  73.          write(ioimp,*) 'IELFIN=',IELFIN
  74.          write(ioimp,*) 'NUM2=',NUM(/2)
  75.          write(ioimp,*) 'MELEME=',MELEME
  76.          CALL ECMAI1(MELEME,0)
  77.          CALL ERREUR(5)
  78.          RETURN
  79.       ENDIF
  80. *
  81.       IF (IMET.EQ.0) THEN
  82.          DO 10 IBELEM=IELDEB,IELFIN
  83. *         WRITE(IOIMP,*) 'IBELEM=',IBELEM
  84. * Volume d'un triangle
  85.             lvirt=.false.
  86.             IF (IDIM.EQ.2) THEN 
  87.                I0=NUM(1,IBELEM)
  88.                I1=NUM(2,IBELEM)
  89.                I2=NUM(3,IBELEM)
  90.                IF (IPVIRT.NE.0) THEN
  91.                   IF (IPVIRT.EQ.I0.OR.IPVIRT.EQ.I1.OR.IPVIRT.EQ.I2) 
  92.      $                 lvirt=.true.
  93. *     $                 goto 10
  94.                ENDIF
  95.                IP0=(I0-1)*IDIMP1
  96.                IP1=(I1-1)*IDIMP1
  97.                IP2=(I2-1)*IDIMP1
  98.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  99.                Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  100.                X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  101.                Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  102. *     XVOL=ABS(DETTRI(X10,Y10,X20,Y20))/2.D0
  103.                XVOLIO=(X10*Y20-X20*Y10)/2.D0
  104. *tst               sigvol=sign(1.d0,xvolio)
  105. *tst               if (sigvol.lt.0.d0) xvoli=xvoli*2
  106. * Volume d'un tétraèdre
  107.             ELSEIF (IDIM.EQ.3) THEN
  108.                I0=NUM(1,IBELEM)
  109.                I1=NUM(2,IBELEM)
  110.                I2=NUM(3,IBELEM)
  111.                I3=NUM(4,IBELEM)
  112.                IF (IPVIRT.NE.0) THEN
  113.                   IF (IPVIRT.EQ.I0.OR.IPVIRT.EQ.I1.OR.IPVIRT.EQ.I2
  114.      $                 .OR.IPVIRT.EQ.I3)
  115.      $                 lvirt=.true.
  116. *     $                 goto 10
  117.                ENDIF
  118.                IP0=(I0-1)*IDIMP1
  119.                IP1=(I1-1)*IDIMP1
  120.                IP2=(I2-1)*IDIMP1
  121.                IP3=(I3-1)*IDIMP1
  122.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  123.                Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  124.                Z10=XCOOR(IP1+3)-XCOOR(IP0+3)
  125.                X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  126.                Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  127.                Z20=XCOOR(IP2+3)-XCOOR(IP0+3)
  128.                X30=XCOOR(IP3+1)-XCOOR(IP0+1)
  129.                Y30=XCOOR(IP3+2)-XCOOR(IP0+2)
  130.                Z30=XCOOR(IP3+3)-XCOOR(IP0+3)
  131.                XVOLIO=(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20
  132.      $              ,Z30))/6.D0
  133. *               XVOL=XVOL+ABS(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20
  134. *     $              ,Z30))
  135.             ENDIF
  136.             xvoli=abs(xvolio)
  137.             if (.not.lvirt) then
  138.                XVOL=XVOL+xvoli
  139.             else
  140.                XVOLV=XVOLV+xvoli
  141.             endif
  142.  10      CONTINUE
  143.       ELSEIF (IMET.EQ.1) THEN
  144.          XMETD=1.D0/DENSIT
  145.          XDETMD=XMETD**IDIM
  146.          DO 20 IBELEM=IELDEB,IELFIN
  147. *         WRITE(IOIMP,*) 'IBELEM=',IBELEM
  148. *     Volume d'un triangle
  149.             lvirt=.false.
  150.             IF (IDIM.EQ.2) THEN 
  151.                I0=NUM(1,IBELEM)
  152.                I1=NUM(2,IBELEM)
  153.                I2=NUM(3,IBELEM)
  154.                IF (IPVIRT.NE.0) THEN
  155.                   IF (IPVIRT.EQ.I0.OR.IPVIRT.EQ.I1.OR.IPVIRT.EQ.I2) 
  156.      $                 lvirt=.true.
  157. *     $                 goto 20
  158.                ENDIF
  159.                IP0=(I0-1)*IDIMP1
  160.                IP1=(I1-1)*IDIMP1
  161.                IP2=(I2-1)*IDIMP1
  162.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  163.                Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  164.                X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  165.                Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  166. *            XVOL=ABS(DETTRI(X10,Y10,X20,Y20))/2.D0
  167. *     XVOLI=ABS(X10*Y20-X20*Y10)*XDETMD
  168.                XVOLIO=((X10*Y20-X20*Y10)/2.D0)*XDETMD
  169. * Volume d'un tétraèdre
  170.             ELSEIF (IDIM.EQ.3) THEN
  171.                I0=NUM(1,IBELEM)
  172.                I1=NUM(2,IBELEM)
  173.                I2=NUM(3,IBELEM)
  174.                I3=NUM(4,IBELEM)
  175.                IF (IPVIRT.NE.0) THEN
  176.                   IF (IPVIRT.EQ.I0.OR.IPVIRT.EQ.I1.OR.IPVIRT.EQ.I2
  177.      $                 .OR.IPVIRT.EQ.I3)
  178.      $                 lvirt=.true.
  179. *     $                 goto 20
  180.                ENDIF
  181.                IP0=(I0-1)*IDIMP1
  182.                IP1=(I1-1)*IDIMP1
  183.                IP2=(I2-1)*IDIMP1
  184.                IP3=(I3-1)*IDIMP1
  185.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  186.                Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  187.                Z10=XCOOR(IP1+3)-XCOOR(IP0+3)
  188.                X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  189.                Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  190.                Z20=XCOOR(IP2+3)-XCOOR(IP0+3)
  191.                X30=XCOOR(IP3+1)-XCOOR(IP0+1)
  192.                Y30=XCOOR(IP3+2)-XCOOR(IP0+2)
  193.                Z30=XCOOR(IP3+3)-XCOOR(IP0+3)
  194. *               XVOLI=ABS(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  195. *     $              *XDETMD
  196.                XVOLIO=((DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  197.      $              /6.D0)*XDETMD
  198.             ENDIF
  199.             xvoli=abs(xvolio)
  200.             if (.not.lvirt) then
  201.                XVOL=XVOL+xvoli
  202.             else
  203.                XVOLV=XVOLV+xvoli
  204.             endif
  205. *            XVOL=XVOL+xvoli
  206.  20      CONTINUE
  207.       ELSE
  208.          WRITE(IOIMP,*) 'vomsi2 imet=',IMET
  209.          CALL ERREUR(5)
  210.          RETURN
  211.       ENDIF
  212.       RETURN
  213. *
  214. * End of subroutine VOMSI2
  215. *
  216.       END
  217.  
  218.  
  219.  

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