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   25/11/24    21:15:28     12406          
  2.       SUBROUTINE VOMSI2(MELEME,IELDEB,IELFIN,NKPVIR,XVOL,XVOLS,XVOLV)
  3.       IMPLICIT REAL*8 (A-H,O-Z)
  4.       IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM         : VOMSI2
  7. C DESCRIPTION : VOlume d'un Maillage de SIMplexes
  8. C               MELEME est supposé actif.
  9. C
  10. C     Comme VOMSIM mais on précise les numéros d'éléments de début
  11. C     et de fin à prendre en compte.
  12. C
  13. C     On sort aussi le volume signe
  14. C
  15. C LANGAGE     : ESOPE
  16. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  17. C               mél : gounand@semt2.smts.cea.fr
  18. C***********************************************************************
  19. C VERSION    : v1, 03/11/2017, version initiale
  20. C HISTORIQUE : v1, 03/11/2017, création
  21. C HISTORIQUE :
  22. C HISTORIQUE :
  23. C***********************************************************************
  24. -INC CCGEOME
  25. -INC PPARAM
  26. -INC CCOPTIO
  27. -INC SMCOORD
  28. -INC SMELEME
  29. -INC SMLREEL
  30.       LOGICAL LVIRT
  31. * Statement functions
  32. *      DETTRI(A11,A12,A21,A22)=A11*A22-A21*A12
  33.       DETTET(A11,A12,A13,A21,A22,A23,A31,A32,A33)=
  34.      &  A11*(A22*A33-A23*A32)
  35.      &       +A12*(A23*A31-A21*A33)
  36.      &       +A13*(A21*A32-A22*A31)
  37. *
  38. * Executable statements
  39. *
  40.       XVOL=0.D0
  41.       XVOLS=0.D0
  42.       XVOLV=0.D0
  43. *Maillage vide
  44. *      SEGACT MELEME
  45.       IF (LISOUS(/1).EQ.0.AND.NUM(/2).EQ.0) THEN
  46.          WRITE(IOIMP,*) 'ITLOC coucou vomsi2'
  47.          write(ioimp,*) 'LISOU1=',LISOUS(/1)
  48.          write(ioimp,*) 'NUM1=',NUM(/1)
  49.          write(ioimp,*) 'NUM2=',NUM(/2)
  50.          write(ioimp,*) 'MELEME=',MELEME
  51.          CALL ECMAI1(MELEME,0)
  52.          CALL ERREUR(5)
  53.          RETURN
  54.       ENDIF
  55.       IDIMP1=IDIM+1
  56.       NBNN=NUM(/1)
  57.       IF (NBNN.NE.IDIMP1.OR.LISOUS(/1).NE.0) THEN
  58.          WRITE(IOIMP,*) 'NBNN=',NBNN
  59.          WRITE(IOIMP,*) 'IDIMP1=',IDIMP1
  60.          WRITE(IOIMP,*) 'NBSOUS=',LISOUS(/1)
  61.          WRITE(IOIMP,*) 'ITLOC'
  62.          CALL ECMAI1(MELEME,0)
  63.          CALL ERREUR(5)
  64.          RETURN
  65.       ENDIF
  66.       IF (.NOT.(IELDEB.GE.1.AND.IELFIN.GE.IELDEB.AND.NUM(
  67.      $     /2).GE.IELFIN))  THEN
  68.          WRITE(IOIMP,*) 'ITLOC coucou vomsi2'
  69.          write(ioimp,*) 'IELDEB=',IELDEB
  70.          write(ioimp,*) 'IELFIN=',IELFIN
  71.          write(ioimp,*) 'NUM2=',NUM(/2)
  72.          write(ioimp,*) 'MELEME=',MELEME
  73.          CALL ECMAI1(MELEME,0)
  74.          CALL ERREUR(5)
  75.          RETURN
  76.       ENDIF
  77. *
  78.       DO 10 IBELEM=IELDEB,IELFIN
  79. *         WRITE(IOIMP,*) 'IBELEM=',IBELEM
  80. * Volume d'un triangle
  81.          lvirt=.false.
  82.          IF (IDIM.EQ.2) THEN
  83.             I0=NUM(1,IBELEM)
  84.             I1=NUM(2,IBELEM)
  85.             I2=NUM(3,IBELEM)
  86.             IF (NKPVIR.NE.0) THEN
  87. *               IF (I0.LE.NKPVIR.OR.I1.LE.NKPVIR.OR.I2.LE.NKPVIR) goto 10
  88.                IF (I0.LE.NKPVIR.OR.I1.LE.NKPVIR.OR.I2.LE.NKPVIR)
  89.      $              lvirt=.true.
  90.             ENDIF
  91.             IP0=(I0-1)*IDIMP1
  92.             IP1=(I1-1)*IDIMP1
  93.             IP2=(I2-1)*IDIMP1
  94.             X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  95.             Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  96.             X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  97.             Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  98.             XVOLIO=(X10*Y20-X20*Y10)/2.D0
  99. * Volume d'un tétraèdre
  100.          ELSEIF (IDIM.EQ.3) THEN
  101.             I0=NUM(1,IBELEM)
  102.             I1=NUM(2,IBELEM)
  103.             I2=NUM(3,IBELEM)
  104.             I3=NUM(4,IBELEM)
  105.             IF (NKPVIR.NE.0) THEN
  106. *               IF (I0.LE.NKPVIR.OR.I1.LE.NKPVIR.OR.I2.LE.NKPVIR.OR.
  107. *     $             I3.LE.NKPVIR) goto 10
  108.                IF (I0.LE.NKPVIR.OR.I1.LE.NKPVIR.OR.I2.LE.NKPVIR.OR.
  109.      $             I3.LE.NKPVIR) lvirt=.true.
  110.             ENDIF
  111.             IP0=(I0-1)*IDIMP1
  112.             IP1=(I1-1)*IDIMP1
  113.             IP2=(I2-1)*IDIMP1
  114.             IP3=(I3-1)*IDIMP1
  115.             X10=XCOOR(IP1+1)-XCOOR(IP0+1)
  116.             Y10=XCOOR(IP1+2)-XCOOR(IP0+2)
  117.             Z10=XCOOR(IP1+3)-XCOOR(IP0+3)
  118.             X20=XCOOR(IP2+1)-XCOOR(IP0+1)
  119.             Y20=XCOOR(IP2+2)-XCOOR(IP0+2)
  120.             Z20=XCOOR(IP2+3)-XCOOR(IP0+3)
  121.             X30=XCOOR(IP3+1)-XCOOR(IP0+1)
  122.             Y30=XCOOR(IP3+2)-XCOOR(IP0+2)
  123.             Z30=XCOOR(IP3+3)-XCOOR(IP0+3)
  124.             XVOLIO=(DETTET(X10,X20,X30,Y10,Y20,Y30,Z10,Z20
  125.      $           ,Z30))/6.D0
  126.          ENDIF
  127.          xvoli=abs(xvolio)
  128.          if (.not.lvirt) then
  129.             XVOL=XVOL+xvoli
  130.          else
  131.             XVOLV=XVOLV+xvoli
  132.          endif
  133.          XVOLS=XVOLS+xvolio
  134.  10   CONTINUE
  135.       RETURN
  136. *
  137. * End of subroutine VOMSI2
  138. *
  139.       END
  140.  
  141.  

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