Télécharger la plaquette Cast3M


fr en

Télécharger vomsi3.eso

Retour à la liste

Numérotation des lignes : 

vomsi3

  1. C VOMSI3    SOURCE    GOUNAND   21/04/06    21:15:45     10940          
  2.       SUBROUTINE VOMSI3(MELEMX,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         : VOMSI3
  8. C DESCRIPTION : VOlume d'un Maillage de SIMplexes
  9. C               MELEMX est supposé actif.
  10. C
  11. C     Comme VOMSI2 mais avec un MELEMX au lieu de MELEME
  12. C
  13. C LANGAGE     : ESOPE
  14. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  15. C               mél : gounand@semt2.smts.cea.fr
  16. C***********************************************************************
  17. C VERSION    : v1, 03/11/2017, version initiale
  18. C HISTORIQUE : v1, 03/11/2017, création
  19. C HISTORIQUE :
  20. C HISTORIQUE :
  21. C***********************************************************************
  22. -INC CCGEOME
  23. -INC PPARAM
  24. -INC CCOPTIO
  25. -INC SMCOORD
  26. -INC TMATOP1
  27. *-INC SMELEMX
  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 MELEMX
  48. *      IF (LISOUS(/1).EQ.0.AND.NUMX(/2).EQ.0) THEN
  49.       IF (NLCOU.EQ.0) THEN
  50.          WRITE(IOIMP,*) 'ITLOC coucou vomsi3'
  51.          write(ioimp,*) 'NLCOU=',NLCOU
  52.          write(ioimp,*) 'MELEMX=',MELEMX
  53.          CALL ECMELX(MELEMX,0)
  54.          CALL ERREUR(5)
  55.          RETURN
  56.       ENDIF
  57.       IDIMP1=IDIM+1
  58. *      NBNN=NUMX(/1)
  59.       NBNN=NNCOU
  60.       IF (NBNN.NE.IDIMP1) THEN
  61.          WRITE(IOIMP,*) 'NBNN=',NBNN
  62.          WRITE(IOIMP,*) 'IDIMP1=',IDIMP1
  63.          write(ioimp,*) 'MELEMX=',MELEMX
  64.          WRITE(IOIMP,*) 'ITLOC'
  65.          CALL ECMELX(MELEMX,0)
  66.          CALL ERREUR(5)
  67.          RETURN
  68.       ENDIF
  69.       IF
  70.      $     (.NOT.(IELDEB.GE.1.AND.IELFIN.GE.IELDEB.AND.NLCOU.GE.IELFIN
  71.      $     .AND.NUMX(/2).GE.NLCOU))  THEN
  72.          WRITE(IOIMP,*) 'ITLOC coucou vomsi3'
  73.          write(ioimp,*) 'IELDEB=',IELDEB
  74.          write(ioimp,*) 'IELFIN=',IELFIN
  75.          write(ioimp,*) 'NLCOU=',NLCOU
  76.          write(ioimp,*) 'NUM2=',NUMX(/2)
  77.          write(ioimp,*) 'MELEMX=',MELEMX
  78.          CALL ECMELX(MELEMX,0)
  79.          CALL ERREUR(5)
  80.          RETURN
  81.       ENDIF
  82. *
  83.       IF (IMET.EQ.0) THEN
  84.          DO 10 IBELEM=IELDEB,IELFIN
  85. *         WRITE(IOIMP,*) 'IBELEM=',IBELEM
  86. * Volume d'un triangle
  87.             lvirt=.false.
  88.             IF (IDIM.EQ.2) THEN 
  89.                I0=NUMX(1,IBELEM)
  90.                I1=NUMX(2,IBELEM)
  91.                I2=NUMX(3,IBELEM)
  92.                IF (IPVIRT.NE.0) THEN
  93.                   IF (IPVIRT.EQ.I0.OR.IPVIRT.EQ.I1.OR.IPVIRT.EQ.I2) 
  94.      $                 lvirt=.true.
  95. *     $                 goto 10
  96.                ENDIF
  97.                IP0=(I0-1)*IDIMP1
  98.                IP1=(I1-1)*IDIMP1
  99.                IP2=(I2-1)*IDIMP1
  100.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  101.                Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  102.                X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  103.                Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  104. *     XVOL=ABS(DETTRI(X10,Y10,X20,Y20))/2.D0
  105.                XVOLIO=(X10*Y20-X20*Y10)/2.D0
  106. *tst               sigvol=sign(1.d0,xvolio)
  107. *tst               if (sigvol.lt.0.d0) xvoli=xvoli*2
  108. * Volume d'un tétraèdre
  109.             ELSEIF (IDIM.EQ.3) THEN
  110.                I0=NUMX(1,IBELEM)
  111.                I1=NUMX(2,IBELEM)
  112.                I2=NUMX(3,IBELEM)
  113.                I3=NUMX(4,IBELEM)
  114.                IF (IPVIRT.NE.0) THEN
  115.                   IF (IPVIRT.EQ.I0.OR.IPVIRT.EQ.I1.OR.IPVIRT.EQ.I2
  116.      $                 .OR.IPVIRT.EQ.I3)
  117.      $                 lvirt=.true.
  118. *     $                 goto 10
  119.                ENDIF
  120.                IP0=(I0-1)*IDIMP1
  121.                IP1=(I1-1)*IDIMP1
  122.                IP2=(I2-1)*IDIMP1
  123.                IP3=(I3-1)*IDIMP1
  124.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  125.                Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  126.                Z10=XCOOR(IP1+3)-XCOOR(IP0+3)
  127.                X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  128.                Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  129.                Z20=XCOOR(IP2+3)-XCOOR(IP0+3)
  130.                X30=XCOOR(IP3+1)-XCOOR(IP0+1)
  131.                Y30=XCOOR(IP3+2)-XCOOR(IP0+2)
  132.                Z30=XCOOR(IP3+3)-XCOOR(IP0+3)
  133.                XVOLIO=(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20
  134.      $              ,Z30))/6.D0
  135. *               XVOL=XVOL+ABS(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20
  136. *     $              ,Z30))
  137.             ENDIF
  138.             xvoli=abs(xvolio)
  139.             if (.not.lvirt) then
  140.                XVOL=XVOL+xvoli
  141.             else
  142.                XVOLV=XVOLV+xvoli
  143.             endif
  144.  10      CONTINUE
  145.       ELSEIF (IMET.EQ.1) THEN
  146.          XMETD=1.D0/DENSIT
  147.          XDETMD=XMETD**IDIM
  148.          DO 20 IBELEM=IELDEB,IELFIN
  149. *         WRITE(IOIMP,*) 'IBELEM=',IBELEM
  150. *     Volume d'un triangle
  151.             lvirt=.false.
  152.             IF (IDIM.EQ.2) THEN 
  153.                I0=NUMX(1,IBELEM)
  154.                I1=NUMX(2,IBELEM)
  155.                I2=NUMX(3,IBELEM)
  156.                IF (IPVIRT.NE.0) THEN
  157.                   IF (IPVIRT.EQ.I0.OR.IPVIRT.EQ.I1.OR.IPVIRT.EQ.I2) 
  158.      $                 lvirt=.true.
  159. *     $                 goto 20
  160.                ENDIF
  161.                IP0=(I0-1)*IDIMP1
  162.                IP1=(I1-1)*IDIMP1
  163.                IP2=(I2-1)*IDIMP1
  164.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  165.                Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  166.                X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  167.                Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  168. *            XVOL=ABS(DETTRI(X10,Y10,X20,Y20))/2.D0
  169. *     XVOLI=ABS(X10*Y20-X20*Y10)*XDETMD
  170.                XVOLIO=((X10*Y20-X20*Y10)/2.D0)*XDETMD
  171. * Volume d'un tétraèdre
  172.             ELSEIF (IDIM.EQ.3) THEN
  173.                I0=NUMX(1,IBELEM)
  174.                I1=NUMX(2,IBELEM)
  175.                I2=NUMX(3,IBELEM)
  176.                I3=NUMX(4,IBELEM)
  177.                IF (IPVIRT.NE.0) THEN
  178.                   IF (IPVIRT.EQ.I0.OR.IPVIRT.EQ.I1.OR.IPVIRT.EQ.I2
  179.      $                 .OR.IPVIRT.EQ.I3)
  180.      $                 lvirt=.true.
  181. *     $                 goto 20
  182.                ENDIF
  183.                IP0=(I0-1)*IDIMP1
  184.                IP1=(I1-1)*IDIMP1
  185.                IP2=(I2-1)*IDIMP1
  186.                IP3=(I3-1)*IDIMP1
  187.                X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  188.                Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  189.                Z10=XCOOR(IP1+3)-XCOOR(IP0+3)
  190.                X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  191.                Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  192.                Z20=XCOOR(IP2+3)-XCOOR(IP0+3)
  193.                X30=XCOOR(IP3+1)-XCOOR(IP0+1)
  194.                Y30=XCOOR(IP3+2)-XCOOR(IP0+2)
  195.                Z30=XCOOR(IP3+3)-XCOOR(IP0+3)
  196. *               XVOLI=ABS(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  197. *     $              *XDETMD
  198.                XVOLIO=((DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20,Z30))
  199.      $              /6.D0)*XDETMD
  200.             ENDIF
  201.             xvoli=abs(xvolio)
  202.             if (.not.lvirt) then
  203.                XVOL=XVOL+xvoli
  204.             else
  205.                XVOLV=XVOLV+xvoli
  206.             endif
  207. *            XVOL=XVOL+xvoli
  208.  20      CONTINUE
  209.       ELSE
  210.          WRITE(IOIMP,*) 'vomsi3 imet=',IMET
  211.          CALL ERREUR(5)
  212.          RETURN
  213.       ENDIF
  214.       RETURN
  215. *
  216. * End of subroutine VOMSI3
  217. *
  218.       END
  219.  
  220.  
  221.  

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