Télécharger la plaquette Cast3M


fr en

Télécharger isova2.eso

Retour à la liste

Numérotation des lignes : 

isova2

  1. C ISOVA2    SOURCE    CB215821  21/06/10    21:15:31     11029          
  2.       SUBROUTINE ISOVA2(MELEME,MCHAML,XISO,XTOL,IOPT,NEWNOD,ELEMS,ITYPL)
  3.       IMPLICIT REAL*8 (A-H,O-Z)
  4.       IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM         : ISOVA2
  7. C DESCRIPTION : Construit le maillage d'une isovaleur d'un champ par
  8. C               éléments
  9. C
  10. C
  11. C
  12. C
  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, 17/12/2008, version initiale
  19. C HISTORIQUE : v1, 17/12/2008, création
  20. C HISTORIQUE : 19/12/2011, bp : ajout cub8
  21. C HISTORIQUE : 10/09/2014, refonte + ajout du traitement 'EGSUP'
  22. C                          'EGINF'  sauf cub8
  23. C HISTORIQUE :
  24. C***********************************************************************
  25.  
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC CCGEOME
  29. -INC SMCOORD
  30. -INC SMELEME
  31. -INC SMCHAML
  32. -INC SMLENTI
  33. *
  34. * Segments ajustables contenant les noeuds des éléments créés ainsi que
  35. *      leur couleur
  36. * ELEM(1) contient des POI1
  37. * ELEM(2) contient des SEG2
  38. * ELEM(3) contient des TRI3
  39. * ELEM(4) contient des TET4
  40. * ELEM(5) contient des PYR5
  41. * ELEM(6) contient des PRI6
  42. * ELEM(7) contient des QUA4
  43. *
  44.       PARAMETER (NTYEL=7)
  45.       SEGMENT ELEMS
  46.       POINTEUR ELEM(NTYEL).MLENTI
  47.       ENDSEGMENT
  48. * Défini dans isova1
  49.       INTEGER ITYPL(NTYEL)
  50. * Pile des nouveaux noeuds
  51.       SEGMENT NEWNOD
  52.       INTEGER NNOD
  53.       INTEGER NOEINF(MAXNOD)
  54.       INTEGER NOESUP(MAXNOD)
  55.       REAL*8  COEINF(MAXNOD)
  56.       ENDSEGMENT
  57. *
  58.       LOGICAL LEGISO,LOKISO
  59.       LOGICAL LPAIR
  60.       PARAMETER (NMAXNO=8)
  61.       REAL*8 VAL(NMAXNO),VALP(NMAXNO)
  62.       LOGICAL LVALP(NMAXNO),LTRI
  63.       INTEGER IPERM(NMAXNO),NUMP(NMAXNO)
  64.       CHARACTER*8 MCHA
  65.       INTEGER IMPR,IRET
  66.  
  67. *     tableaux descriptifs des cas et des aretes
  68. *     pour le Marching Cubes (CUB8)
  69.       PARAMETER (NMAXPT=16,NMAXCA=256,npt2=2,NMAXAR=12)
  70.       INTEGER MATCAS(NMAXPT,NMAXCA),MATARE(npt2,NMAXAR)
  71. *     tableaux descriptifs des cas
  72. *     pour les éléments POI1, SEG2, TRI3, TET4
  73. *
  74. *  ikas                       pile elem
  75. * *************************** Commun à tous ********************
  76. *  1 : rien                      0
  77. *  2 : élément total            9999
  78. * ***************************** POI1 ***************************
  79. *  3 : noeud 1                   1,  1,0
  80. * ***************************** SEG2 ***************************
  81. *= 3 : noeud 1                   1,  1,0
  82. *  4 : noeud 2                   1,  2,0
  83. *  5 : noeud (1,2)               1,  1,2
  84. *  6 : segment [1, k5]           2,  1,0, 5,5,
  85. *  7 : segment [k5, 2]           2,  5,5, 2,0,
  86. * ***************************** TRI3 ***************************
  87. *
  88. *                 3                            3
  89. *               .                            .
  90. *              / \                          / \
  91. *             /   \                        /   \
  92. *            /     \                      /     \
  93. *     (1,3) /       \              (1,3) /       \ (2,3)
  94. *          .-        \                  .---------.
  95. *         /  \--      \                / \         \
  96. *        /      \-     \              /   \         \
  97. *       /         \--   \            /     \         \
  98. *      /             \-- \          /       \         \
  99. *   1 /                 \-\ 2    1 /         \         \ 2
  100. *    .--------------------\.      .-----------.---------.
  101. *                                           (1,2)
  102. *
  103. *= 3 : noeud 1                   1,  1,0
  104. *  8 : noeud 3                   1,  3,0
  105. *  9 : segment [1,2]             2,  1,0,   2,0
  106. * 10 : segment [3,2]             2,  3,0,   2,0
  107. * 11 : segment [(1,3),2]         2,  1,3,   2,0
  108. * 12 : triangle [1, k11]         3,  1,0,   2,0,  1,3,
  109. * 13 : triangle [ik11, 3]        3,  1,3,   2,0,  3,0,
  110. * 14 : segment [(1,3),(1,2)]     2,  1,3,   1,2
  111. * 15 : triangle [1, ik14]        3,  1,0,   1,2,  1,3,
  112. * 16 : quadrangle [k14, 2,3]     7,  1,3,   1,2,  2,0,  3,0,
  113. * 17 : segment [(1,3),(2,3)]     2,  1,3,   2,3
  114. * 18 : triangle [k17, 3]         3,  1,3,   2,3,  3,0,
  115. * 19 : quadrangle [1,2,ik17]     7,  1,0,   2,0,  2,3,  1,3,
  116. * ***************************** TET4 ***************************
  117. *   Dur de faire des graphiques ASCII clairs en 3D.
  118. *
  119. *                 4  .-
  120. *                  /-| \--
  121. *                 /   \   \--
  122. *               /-    |      \---
  123. *              /       \         \--
  124. *            /-        |            \--
  125. *           /           \           ---\. 2
  126. *         /-            |   -------/   /
  127. *        /         ------\-/          /
  128. *      /-  -------/      |           /
  129. *   3 .---/              |          /
  130. *        \---             \       /-
  131. *            \---         |      /
  132. *                \--       \    /
  133. *                   \---   |   /
  134. *                       \---\ /
  135. *                            . 1
  136. *
  137. *= 3 : noeud 1                   1,  1,0
  138. * 20 : noeud 4                   1,  4,0,
  139. *= 9 : segment [1,2]             2,  1,0,   2,0
  140. * 21 : segment [4,3]             2,  4,0,   3,0
  141. * 22 : triangle [1,2,3]          3,  1,0,   2,0,   3,0,
  142. * 23 : triangle [4,2,3]          3,  4,0,   2,0,   3,0,
  143. * 24 : triangle [(1,4),2,3]      3,  1,4,   2,0,   3,0,
  144. * 25 : tétra    [1, k24]         4,  1,0,   2,0,   3,0,  1,4,
  145. * 26 : tétra    [k24, 4]         4,  1,4,   2,0,   3,0,  4,0,
  146. ******************************
  147. * 27 : triangle [(1,4),2,(1,3)]  3,  1,4,   2,0,   1,3,
  148. * 28 : tétra    [1,k27]          4,  1,0,   2,0,   1,3,   1,4,
  149. * 29 : pyra     [4,3,k27]        5,  1,4,   1,3,   3,0,   4,0,   2,0,
  150. * 30 : triangle [(1,4),(2,4),3]  3,  1,4,   2,4,   3,0,
  151. * 31 : tétra    [k30,4]          4,  1,4,   2,4,   3,0,   4,0,
  152. * 32 : pyra     [1,2,k30]        5,  1,4,   2,4,   2,0,   1,0,   3,0,
  153. ***************************************
  154. * 33 : triangle [(1,2),(1,3),(1,4)]  3,   1,2,   1,3,   1,4,
  155. * 34 : tétra    [1,k33]
  156. * 35 : prisme   [k33,2,3,4]
  157. * 36 : triangle [(1,4),(2,4),(3,4)]
  158. * 37 : tétra    [k36,4]
  159. * 38 : prisme   [1,2,3,k36]
  160. * 39 : quadrangle [(1,4),(2,4),(2,3),(1,3)]
  161. * 40 : prisme   [(1,3),(2,3),3,(1,4),(2,4),4]
  162. * 41 : prisme   [1,(1,3),(1,4),2,(2,3),(2,4)]
  163.       PARAMETER (NLI1=7,NLI2=12,NLI3=7,NLI4=15)
  164.       PARAMETER (NLI1P=NLI1+1,NLI2P=NLI1+NLI2+1,NLI3P=NLI1+NLI2+NLI3+1)
  165.       PARAMETER (NCOMAX=1+(2*6),NLIMAX=NLI1+NLI2+NLI3+NLI4)
  166.       INTEGER TABDES(NCOMAX,NLIMAX)
  167.       INTEGER TABDE1(NCOMAX,NLI1)
  168.       INTEGER TABDE2(NCOMAX,NLI2)
  169.       INTEGER TABDE3(NCOMAX,NLI3)
  170.       INTEGER TABDE4(NCOMAX,NLI4)
  171.       EQUIVALENCE (TABDES(1,1),TABDE1(1,1)),
  172.      $     (TABDES(1,NLI1P),TABDE2(1,1)),
  173.      $     (TABDES(1,NLI2P),TABDE3(1,1)),
  174.      $     (TABDES(1,NLI3P),TABDE4(1,1))
  175. *
  176.       DATA TABDE1/
  177.      $     0,   0,0,   0,0,   0,0,   0,0,   0,0,   0,0,
  178.      $     9999,   0,0,   0,0,   0,0,   0,0,   0,0,   0,0,
  179.      $     1,   1,0,   0,0,   0,0,   0,0,   0,0,   0,0,
  180.      $     1,   2,0,   0,0,   0,0,   0,0,   0,0,   0,0,
  181.      $     1,   1,2,   0,0,   0,0,   0,0,   0,0,   0,0,
  182.      $     2,   1,0,   1,2,   0,0,   0,0,   0,0,   0,0,
  183.      $     2,   1,2,   2,0,   0,0,   0,0,   0,0,   0,0/
  184.       DATA TABDE2/
  185.      $     1,   3,0,   0,0,   0,0,   0,0,   0,0,   0,0,
  186.      $     2,   1,0,   2,0,   0,0,   0,0,   0,0,   0,0,
  187.      $     2,   3,0,   2,0,   0,0,   0,0,   0,0,   0,0,
  188.      $     2,   1,3,   2,0,   0,0,   0,0,   0,0,   0,0,
  189.      $     3,   1,0,   2,0,   1,3,   0,0,   0,0,   0,0,
  190.      $     3,   1,3,   2,0,   3,0,   0,0,   0,0,   0,0,
  191.      $     2,   1,3,   1,2,   0,0,   0,0,   0,0,   0,0,
  192.      $     3,   1,0,   1,2,   1,3,   0,0,   0,0,   0,0,
  193.      $     7,   1,3,   1,2,   2,0,   3,0,   0,0,   0,0,
  194.      $     2,   1,3,   2,3,   0,0,   0,0,   0,0,   0,0,
  195.      $     3,   1,3,   2,3,   3,0,   0,0,   0,0,   0,0,
  196.      $     7,   1,0,   2,0,   2,3,   1,3,   0,0,   0,0/
  197.       DATA TABDE3/
  198.      $     1,   4,0,   0,0,   0,0,   0,0,   0,0,   0,0,
  199.      $     2,   4,0,   3,0,   0,0,   0,0,   0,0,   0,0,
  200.      $     3,   1,0,   2,0,   3,0,   0,0,   0,0,   0,0,
  201.      $     3,   4,0,   2,0,   3,0,   0,0,   0,0,   0,0,
  202.      $     3,   1,4,   2,0,   3,0,   0,0,   0,0,   0,0,
  203.      $     4,   1,0,   2,0,   3,0,   1,4,   0,0,   0,0,
  204.      $     4,   1,4,   2,0,   3,0,   4,0,   0,0,   0,0/
  205.       DATA TABDE4/
  206.      $     3,   1,4,   2,0,   1,3,   0,0,   0,0,   0,0,
  207.      $     4,   1,0,   2,0,   1,3,   1,4,   0,0,   0,0,
  208.      $     5,   1,4,   1,3,   3,0,   4,0,   2,0,   0,0,
  209.      $     3,   1,4,   2,4,   3,0,   0,0,   0,0,   0,0,
  210.      $     4,   1,4,   2,4,   3,0,   4,0,   0,0,   0,0,
  211.      $     5,   1,4,   2,4,   2,0,   1,0,   3,0,   0,0,
  212.      $     3,   1,2,   1,3,   1,4,   0,0,   0,0,   0,0,
  213.      $     4,   1,0,   1,2,   1,3,   1,4,   0,0,   0,0,
  214.      $     6,   1,2,   1,3,   1,4,   2,0,   3,0,   4,0,
  215.      $     3,   1,4,   2,4,   3,4,   0,0,   0,0,   0,0,
  216.      $     4,   1,4,   2,4,   3,4,   4,0,   0,0,   0,0,
  217.      $     6,   1,0,   2,0,   3,0,   1,4,   2,4,   3,4,
  218.      $     7,   1,4,   2,4,   2,3,   1,3,   0,0,   0,0,
  219.      $     6,   1,3,   2,3,   3,0,   1,4,   2,4,   4,0,
  220.      $     6,   1,0,   1,3,   1,4,   2,0,   2,3,   2,4/
  221. *
  222. * Fonctions locales
  223. *
  224. * Cette fonction dit si on est sur une isovaleur
  225. * Cette fonction dit ok si on vérifie .GT.XISO si IOPT=1
  226. *                                ou   .LT.XISO si IOPT différent de 1
  227. *                                              (typiquement -1)
  228.       LEGISO(X)=(ABS(X-XISO).LT.XTOL)
  229. * Suite à un bug du traducteur Esope, on ne peut déclarer
  230. * deux statement functions dans une subroutine. On remplace
  231. * donc celle-ci par sa valeur :
  232. *      LOKISO(X)=((X.GT.XISO).EQV.(IOPT.EQ.1))
  233. *
  234. * Executable statements
  235. *
  236.       IMPR=0
  237.       IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans isova2.eso'
  238. *
  239. C***********************************************************************
  240. c     Boucle sur les zones elementaires
  241. C***********************************************************************
  242. c
  243.       NNO=NUM(/1)
  244.       NEL=NUM(/2)
  245.       N2=IELVAL(/1)
  246.       DO I2=1,N2
  247.  
  248.          MELVAL=IELVAL(I2)
  249.          SEGACT MELVAL
  250.  
  251. * Vérification de la cohérence entre la champ et la géométrie
  252.          N1PTEL=VELCHE(/1)
  253.          N1EL=VELCHE(/2)
  254.          IF ((N1PTEL.NE.1.AND.N1PTEL.NE.NNO).OR.
  255.      $        (N1EL.NE.1.AND.N1EL.NE.NEL)) THEN
  256. * 148 2
  257. *Le champ par élément et l'objet géométrique correspondant n'ont pas le
  258. * 148 2
  259. *même nombre de points par élément : contacter votre support
  260.             CALL ERREUR(148)
  261.             RETURN
  262.          ENDIF
  263. * Au maximum, on va créer quatre nouveaux noeuds par élément
  264.          MAXNOC=NOEINF(/1)
  265.          MAXNOD=NNOD+4*NEL
  266.          IF (MAXNOD.GT.MAXNOC) SEGADJ NEWNOD
  267. C
  268. c******* Aiguillage selon type d element *******************************
  269. c
  270.          CALL PLACE2(ITYPL,NTYEL,ITYLOC,ITYPEL)
  271. *dbg         WRITE(IOIMP,*) 'ITYLOC=',ITYLOC
  272.          IF (ITYLOC.GE.1.AND.ITYLOC.LE.4) GOTO 101
  273. *         IF(ITYPEL.eq.1)  GOTO 101
  274. *         IF(ITYPEL.eq.2)  GOTO 102
  275. *         IF(ITYPEL.eq.4)  GOTO 104
  276. * Cub8
  277.          if(ITYPEL.eq.14)  goto 114
  278. *         if(ITYPEL.eq.23)  goto 123
  279.  
  280. ******** Cas non traité *********************************************
  281.          MOTERR(1:4)=NOMS(ITYPEL)
  282. *  44 2
  283. *Type d'élément inconnu %m1:4
  284.          CALL ERREUR(44)
  285.          RETURN
  286. ******** Cas des POI1, SEG2, TRI3 et TET4
  287.  101     CONTINUE
  288.          NNODEL=NBNNE(ITYPEL)
  289. *dbg         WRITE(IOIMP,*) 'NNODEL=',NNODEL
  290.          DO IEL=1,NEL
  291. * Par defaut, il n'y a pas d'isovaleur a construire
  292. * et l'orientation est bonne
  293.             IKAS=1
  294.             LPAIR=.TRUE.
  295. * Remplissage du tableau VAL, calcul du min et de max
  296.             DO INODEL=1,NNODEL
  297.                VAL(INODEL)=VELCHE(MIN(N1PTEL,INODEL),MIN(N1EL,IEL))
  298.             ENDDO
  299.             VMIN=VAL(1)
  300.             VMAX=VAL(1)
  301.             DO INODEL=2,NNODEL
  302.                VMIN=MIN(VMIN,VAL(INODEL))
  303.                VMAX=MAX(VMAX,VAL(INODEL))
  304.             ENDDO
  305. *            WRITE(IOIMP,*) 'XISO=',XISO,' VMIN=',VMIN
  306. *     $           ,' VMAX=',VMAX
  307. *     Y a-t-il une isovaleur ?
  308.             IF (((VMIN-XISO).LE.XTOL).AND.((VMAX-XISO).GE.-XTOL))
  309.      $           THEN
  310. * Permutation permettant d'avoir VAL en ordre croissant
  311. * VAL(IPERM(1)).LE.VAL(IPERM(2)).LE.VAL(IPERM(3))
  312. * On garde la parité de la permutation car on va essayer d'avoir une
  313. * orientation consistante des éléments crées :
  314. * - Conserver l'orientation des éléments volumiques
  315. * - Normale aux éléments surfaciques dans le sens du gradient
  316. *
  317. *obs  On ne permute pas dans le cas des SEG2 car on souhaite garder
  318. *obs  l'orientation. Pour les éléments de dimension plus haute, on
  319. *obs  ne garde pas l'orientation (sinon cela semble multiplier le nombre de
  320. *obs  cas).
  321.                DO INODEL=1,NNODEL
  322.                   IPERM(INODEL)=INODEL
  323.                ENDDO
  324. *               IF (ITYLOC.NE.2) THEN
  325.                CALL ISHELR(NNODEL,IPERM,NNODEL,VAL,NINV,IMPR,IRET)
  326.                LPAIR=(MOD(NINV,2).EQ.0)
  327. *               ENDIF
  328. *               Write(ioimp,*) 'Iperm=',(IPERM(I),I=1,NNODEL)
  329. *               Write(ioimp,*) 'Ninv=',ninv
  330. *
  331. * Nombre de points de l'élément ayant "exactement" la valeur XISO
  332. * demandee
  333.                NPEQ1=0
  334.                DO INODEL=1,NNODEL
  335.                   VALP(INODEL)=VAL(IPERM(INODEL))
  336.                   NUMP(INODEL)=NUM(IPERM(INODEL),IEL)
  337.                   LVALP(INODEL)=(LEGISO(VALP(INODEL)))
  338.                   IF (LVALP(INODEL)) THEN
  339.                      NPEQ1=NPEQ1+1
  340.                   ENDIF
  341.                ENDDO
  342. *dbg               WRITE(IOIMP,*) 'NPEQ1=',NPEQ1,' XISO=',XISO,' XTOL=',XTOL
  343. ************************************************************************
  344. **************                                                 *********
  345. **************    Calcul du IKAS pour chaque type d'élément    *********
  346. **************                                                 *********
  347. ************************************************************************
  348. **************                                                 *********
  349. **************                 POI1                            *********
  350.                IF (ITYLOC.EQ.1) THEN
  351.                   IF (NPEQ1.EQ.1) THEN
  352.                      IKAS=2
  353.                   ENDIF
  354. **************                                                 *********
  355. **************                                                 *********
  356. **************                 SEG2                            *********
  357. **************                                                 *********
  358.                ELSEIF (ITYLOC.EQ.2) THEN
  359.                   IF (NPEQ1.EQ.2) THEN
  360. * Le segment entier
  361.                      IKAS=2
  362.                   ELSEIF (NPEQ1.EQ.1) THEN
  363. * Une des extrémités
  364.                      IF (LVALP(1)) THEN
  365. * Le premier point
  366. * ou le segment entier (IOPT.EQ.1)  ET (VALP(2).GT.XISO)
  367. *                   ou (IOPT.EQ.-1) ET (VALP(2).LT.XISO)
  368. *                        IF (LOKISO(VALP(2))) THEN
  369.                         IF ((IOPT.NE.0).AND.
  370.      $                       ((VALP(2).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  371.                            IKAS=2
  372.                         ELSE
  373.                            IKAS=3
  374.                         ENDIF
  375.                      ELSE
  376. * Le deuxième point
  377. * ou le segment entier (IOPT.EQ.1)  ET (VALP(1).GT.XISO)
  378. *                   ou (IOPT.EQ.-1) ET (VALP(1).LT.XISO)
  379. *                        IF (LOKISO(VALP(1))) THEN
  380.                         IF ((IOPT.NE.0).AND.
  381.      $                       ((VALP(1).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  382.                            IKAS=2
  383.                         ELSE
  384.                            IKAS=4
  385.                         ENDIF
  386.                      ENDIF
  387.                   ELSE
  388. * Un point milieu (IOPT=0)
  389. * ou un bout de segment (IOPT.NE.0)
  390.                      IF (IOPT.EQ.0) THEN
  391.                         IKAS=5
  392.                      ELSE
  393. *     On garde l'orientation du segment
  394. *                     IF (LOKISO(VALP(1))) THEN
  395.                         IF ((VALP(1).GT.XISO).EQV.(IOPT.EQ.1)) THEN
  396.                            IKAS=6
  397.                         ELSE
  398.                            IKAS=7
  399.                         ENDIF
  400.                      ENDIF
  401.                   ENDIF
  402. **************                                                 *********
  403. **************                 TRI3                            *********
  404. **************                                                 *********
  405.                ELSEIF (ITYLOC.EQ.3) THEN
  406.                   IF (NPEQ1.EQ.3) THEN
  407. * Le triangle entier
  408.                      IKAS=2
  409.                   ELSEIF (NPEQ1.EQ.2) THEN
  410. * Un des côtés
  411. * ou le triangle entier (IOPT.EQ.1)  ET (VALP(3 ou 1).GT.XISO)
  412. *                       (IOPT.EQ.-1) ET (VALP(3 ou 1).LT.XISO)
  413.                      IF (LVALP(1)) THEN
  414.                         IF ((IOPT.NE.0).AND.
  415.      $                       ((VALP(3).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  416. * Le triangle entier
  417.                            IKAS=2
  418.                         ELSE
  419. * Le segment [1,2]
  420.                            IKAS=9
  421.                         ENDIF
  422.                      ELSE
  423.                         IF ((IOPT.NE.0).AND.
  424.      $                       ((VALP(1).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  425. * Le triangle entier
  426.                            IKAS=2
  427.                         ELSE
  428. * Le segment [2,3]
  429.                            IKAS=10
  430.                         ENDIF
  431.                      ENDIF
  432.                   ELSEIF (NPEQ1.EQ.1) THEN
  433. *     Un des points                  ou un segment
  434. * ou  un des points ou le triangle   ou un triangle contenant le segment
  435.  
  436.                      IF (LVALP(1)) THEN
  437. *  Le point P1
  438. *  ou le triangle entier
  439.                         IF ((IOPT.NE.0).AND.
  440.      $                       ((VALP(2).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  441. *  Le triangle entier
  442.                            IKAS=2
  443.                         ELSE
  444. *  3 : noeud 1
  445.                            IKAS=3
  446.                         ENDIF
  447.                      ELSEIF (LVALP(2)) THEN
  448. * Un segment reliant P2 à l'intersection sur (P1, P3)
  449. * ou un triangle contenant ce segment (IOPT.NE.0)
  450.                         IF (IOPT.EQ.0) THEN
  451. * 11 : segment [2,(1,3)]
  452.                            IKAS=11
  453.                         ELSE
  454.                            IF ((VALP(1).GT.XISO).EQV.(IOPT.EQ.1)) THEN
  455. * 12 : triangle [1, k11]
  456.                               IKAS=12
  457.                            ELSE
  458. * 13 : triangle [k11, 3]
  459.                               IKAS=13
  460.                            ENDIF
  461.                         ENDIF
  462.                      ELSE
  463. *  Le point P3
  464. *  ou le triangle entier
  465.                         IF ((IOPT.NE.0).AND.
  466.      $                       ((VALP(2).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  467. *  Le triangle entier
  468.                            IKAS=2
  469.                         ELSE
  470. *  8 : noeud 3
  471.                            IKAS=8
  472.                         ENDIF
  473.                      ENDIF
  474.                   ELSE
  475. * Un segment reliant (P1,P2) à (P1,P3)
  476. *                 ou (P2,P3) à (P1,P3)
  477. * ou un triangle ou un quadrangle contenant ces segments
  478.                      IF (IOPT.EQ.0) THEN
  479.                         MLENTI=ELEM(2)
  480.                         IF (VALP(2).GT.XISO) THEN
  481. * Un segment reliant (P1,P2) à (P1,P3)
  482. * 14 : segment [(1,2),(1,3)]
  483.                            IKAS=14
  484.                         ELSE
  485. * Un segment reliant (P2,P3) à (P1,P3)
  486. * 17 : segment [(2,3),(1,3)]
  487.                            IKAS=17
  488.                         ENDIF
  489.                      ELSE
  490. * Un triangle ou un quadrangle s'appuyant sur le
  491. * segment reliant (P1,P2) à (P1,P3)
  492.                         IF (VALP(2).GT.XISO) THEN
  493.                            IF ((VALP(1).GT.XISO).EQV.(IOPT.EQ.1)) THEN
  494. * 15 : triangle [1, k14]
  495.                               IKAS=15
  496.                            ELSE
  497. * 16 : quadrangle [k14, 3,2]
  498.                               IKAS=16
  499.                            ENDIF
  500.                         ELSE
  501. * Un triangle ou un quadrangle s'appuyant sur le
  502. * segment reliant (P2,P3) à (P1,P3)
  503.                            IF ((VALP(3).GT.XISO).EQV.(IOPT.EQ.1)) THEN
  504. * 18 : triangle [3, k17]
  505.                               IKAS=18
  506.                            ELSE
  507. * 19 : quadrangle [1,2,k17]
  508.                               IKAS=19
  509.                            ENDIF
  510.                         ENDIF
  511.                      ENDIF
  512.                   ENDIF
  513. **************                                                 *********
  514. **************                 TET4                            *********
  515. **************                                                 *********
  516.                ELSEIF (ITYLOC.EQ.4) THEN
  517.                IF (NPEQ1.EQ.4) THEN
  518. *  Isovaleur = tetrahedre
  519.                   IKAS=2
  520.                ELSEIF (NPEQ1.EQ.3) THEN
  521. *  Isovaleur = 1 face du tetrahedre
  522. * ou le tetraedre entier (IOPT.EQ.1)  ET (VALP(4 ou 1).GT.XISO)
  523. *                        (IOPT.EQ.-1) ET (VALP(4 ou 1).LT.XISO)
  524.                   IF (LVALP(1)) THEN
  525.                      IF ((IOPT.NE.0).AND.
  526.      $                    ((VALP(4).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  527. * le tetraedre entier
  528.                         IKAS=2
  529.                      ELSE
  530. * 22 : triangle [1,2,3]
  531.                         IKAS=22
  532.                      ENDIF
  533.                   ELSE
  534.                      IF ((IOPT.NE.0).AND.
  535.      $                    ((VALP(1).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  536. * le tetraedre entier
  537.                         IKAS=2
  538.                      ELSE
  539. * 23 : triangle [2,3,4]
  540.                         IKAS=23
  541.                      ENDIF
  542.                   ENDIF
  543.                ELSEIF (NPEQ1.EQ.2) THEN
  544. *  Isovaleur = 1 segment (P1,P2) ou (P3,P4) ou 1 triangle
  545. *        appuye sur 2 points du tetrahedre (P2,P3)
  546. *  ou (IOPT.NE.0)  le segment ou le tétra complet
  547. *                  ou un tétra appuye sur le triangle appuye sur (P2,P3)
  548. *                  IF (LVALP(1).OR.LVALP(4)) THEN
  549. *  Isovaleur = 1 segment  (P1,P2) ou (P3,P4)
  550. *  ou le tétra complet
  551.                   IF (LVALP(1)) THEN
  552.                      IF ((IOPT.NE.0).AND.
  553.      $                    ((VALP(3).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  554. * le tetraedre entier
  555.                         IKAS=2
  556.                      ELSE
  557. *  9 : segment [1,2]
  558.                         IKAS=9
  559.                      ENDIF
  560.                   ELSEIF (LVALP(4)) THEN
  561.                      IF ((IOPT.NE.0).AND.
  562.      $                    ((VALP(2).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  563. * le tetraedre entier
  564.                         IKAS=2
  565.                      ELSE
  566. * 21 : segment [3,4]
  567.                         IKAS=21
  568.                      ENDIF
  569.                   ELSE
  570. *  Isovaleur = 1 triangle appuye sur 2 points tetra (P2,P3)
  571. * ou un tétra contenant ce triangle (IOPT.NE.0)
  572.                      IF (IOPT.EQ.0) THEN
  573. * 24 : triangle [2,3,(1,4)]
  574.                         IKAS=24
  575.                      ELSE
  576.                         IF ((VALP(1).GT.XISO).EQV.(IOPT.EQ.1)) THEN
  577. * 25 : tétra    [1, k24]
  578.                            IKAS=25
  579.                         ELSE
  580. * 26 : tétra    [k24, 4]
  581.                            IKAS=26
  582.                         ENDIF
  583.                      ENDIF
  584.                   ENDIF
  585.                ELSEIF (NPEQ1.EQ.1) THEN
  586. *  Isovaleur = 1 point ou 1 triangle appuye sur 1 point du tetrahedre
  587. *  ou (IOPT.NE.0) 1 point ou le tétraèdre
  588. *                 ou  1 tétra ou 1 pyra
  589. *                  IF (LVALP(1).OR.LVALP(4)) THEN
  590. *  Isovaleur = 1 point ou le tétraèdre
  591.                   IF (LVALP(1)) THEN
  592.                      IF ((IOPT.NE.0).AND.
  593.      $                    ((VALP(2).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  594. * le tetraedre entier
  595.                         IKAS=2
  596.                      ELSE
  597. *  3 : noeud 1
  598.                         IKAS=3
  599.                      ENDIF
  600.                   ELSEIF (LVALP(4)) THEN
  601.                      IF ((IOPT.NE.0).AND.
  602.      $                    ((VALP(3).GT.XISO).EQV.(IOPT.EQ.1))) THEN
  603. * le tetraedre entier
  604.                         IKAS=2
  605.                      ELSE
  606. * 20 : noeud 4
  607.                         IKAS=20
  608.                      ENDIF
  609.                   ELSE
  610. *  Isosurface = 1 triangle appuye sur 1 point du tetrahedre
  611. *      ou 1 tétra ou 1 pyra
  612.                      IF (LVALP(2)) THEN
  613. *  Le point appuye est P2
  614. *  les autres points intersection (P1,P3) et (P1,P4)
  615.                         IF (IOPT.EQ.0) THEN
  616. * 27 : triangle [(1,3),(1,4),2]
  617.                            IKAS=27
  618.                         ELSE
  619. *   tétra contenant le triangle précédent et P1
  620.                            IF ((VALP(1).GT.XISO).EQV.(IOPT.EQ.1)) THEN
  621. * 28 : tétra    [1,k27]
  622.                               IKAS=28
  623.                            ELSE
  624. *   pyra contenant le triangle précédent et (P3,P4) (sommet P2)
  625. * 29 : pyra     [4,3,k27]
  626.                               IKAS=29
  627.                            ENDIF
  628.                         ENDIF
  629.                      ELSE
  630. *  Le point appuye est P3
  631. *  les autres points intersection (P2,P4) et (P1,P4)
  632.                         IF (IOPT.EQ.0) THEN
  633. * 30 : triangle [(1,4),(2,4),3]
  634.                            IKAS=30
  635.                         ELSE
  636. *   tétra contenant le triangle précédent et P4
  637.                            IF ((VALP(4).GT.XISO).EQV.(IOPT.EQ.1)) THEN
  638. * 31 : tétra    [k30,4]
  639.                               IKAS=31
  640.                            ELSE
  641. *   pyra contenant le triangle précédent et (P1,P2) (sommet P3)
  642. * 32 : pyra     [1,2,k30]
  643.                               IKAS=32
  644.                            ENDIF
  645.                         ENDIF
  646.                      ENDIF
  647.                   ENDIF
  648.                ELSE
  649. *  Cas général : pas de point du tétra
  650. *  Isosurface = 1 triangle quelconque section
  651. *   du tetrahedre ou un quadrangle   (= 2 triangles)
  652. *  ou (IOPT.NE.0)  1 tétra ou 1 prisme s'appuyant sur le triangle
  653. *              ou  1 prisme s'appuyant sur le quadrangle
  654. *
  655. *  Il y toujours un point intersection
  656. *  entre P1 et P4
  657.                   IF (VALP(2).GT.XISO) THEN
  658. *  Isosurface entre V1 et V2
  659. *  Points intersection (P1,P2) (P1,P3) (P1,P4)
  660. *  Un seul triangle
  661.                      IF (IOPT.EQ.0) THEN
  662. * 33 : triangle [(1,2),(1,3),(1,4)]
  663.                         IKAS=33
  664.                      ELSE
  665. * 1 tétra sommet P1 base (P1,P2) (P1,P3) (P1,P4)
  666.                         IF ((VALP(1).GT.XISO).EQV.(IOPT.EQ.1)) THEN
  667. * 34 : tétra    [k33,1]
  668.                            IKAS=34
  669.                         ELSE
  670. * 1 prisme base (P1,P2) (P1,P3) (P1,P4) et P2 P3 P4
  671. * 35 : prisme   [k33,2,3,4]
  672.                            IKAS=35
  673.                         ENDIF
  674.                      ENDIF
  675.                   ELSEIF (VALP(3).LT.XISO) THEN
  676. *  Isosurface entre V3 et V4
  677. *  Points intersection (P3,P4) (P2,P4) (P1,P4)
  678. *  Un seul triangle
  679.                      IF (IOPT.EQ.0) THEN
  680. * 36 : triangle [(1,4),(2,4),(3,4)]
  681.                         IKAS=36
  682.                      ELSE
  683. * 1 tétra sommet P4 base (P3,P4) (P2,P4) (P1,P4)
  684.                         IF ((VALP(4).GT.XISO).EQV.(IOPT.EQ.1)) THEN
  685. * 37 : tétra    [k36,4]
  686.                            IKAS=37
  687.                         ELSE
  688. * 1 prisme base (P3,P4) (P2,P4) (P1,P4) et P3 P2 P1
  689. * 38 : prisme   [k36,1,2,3]
  690.                            IKAS=38
  691.                         ENDIF
  692.                      ENDIF
  693.                   ELSE
  694. *  Isosurface entre V2 et V3 ;
  695. *  Points intersection (P1,P3) (P2,P3) (P2,P4) (P1,P4)
  696. *  Un quadrangle = 2 triangles
  697.                      IF (IOPT.EQ.0) THEN
  698. * 39 : quadrangle [(1,4),(2,4),(2,3),(1,3)]
  699.                         IKAS=39
  700.                      ELSE
  701. * 1 prisme base P4 (P1,P4) (P2,P4) et P3 (P1,P3) (P2,P3)
  702.                         IF ((VALP(4).GT.XISO).EQV.(IOPT.EQ.1)) THEN
  703. * 40 : prisme   [4,(1,4),(2,4),3,(1,3),(2,3)]
  704.                            IKAS=40
  705.                         ELSE
  706. * 1 prisme base P1 (P1,P3) (P1,P4) et P2 (P2,P3) (P2,P4)
  707. * 41 : prisme   [1,(1,3),(1,4),2,(2,3),(2,4)]
  708.                            IKAS=41
  709.                         ENDIF
  710.                      ENDIF
  711.                   ENDIF
  712.                ENDIF
  713. **************                                                 *********
  714. **************             Element inconnu                     *********
  715. **************                                                 *********
  716.                ELSE
  717.                   WRITE(IOIMP,*) 'ITYLOC=',ITYLOC
  718.                   GOTO 9999
  719.                ENDIF
  720. *
  721. * Doit-on garder quand même l'élément entier ?
  722. *            ELSEIF (IOPT.NE.0.AND.LOKISO(VAL(1))) THEN
  723.             ELSEIF (IOPT.NE.0.AND.((VAL(1).GT.XISO).EQV.(IOPT.EQ.1)))
  724.      $              THEN
  725.                IKAS=2
  726.             ENDIF
  727.  
  728. **************************************************************************
  729. *
  730. *           Décodage du cas (IKAS) et création des éléments voulus
  731. *
  732. **************************************************************************
  733. *dbg            WRITE(IOIMP,*) 'IKAS=',IKAS
  734. *         type d'élément à créer
  735.             ITYL=TABDES(1,IKAS)
  736. *dbg            WRITE(IOIMP,*) '  IEL=',IEL,' ITYL=',ITYL
  737. * Rien à faire
  738.             IF (ITYL.EQ.0) THEN
  739. * L'élément entier, on garde la même description (orientation...)
  740. * que l'element initial (utilisation de NUM et non NUMP)
  741.             ELSEIF (ITYL.EQ.9999) THEN
  742.                CALL PLACE2(ITYPL,NTYEL,ITYL2,ITYPEL)
  743.                IF (ITYL2.GE.1.AND.ITYL2.LE.NTYEL) THEN
  744. *                  WRITE(IOIMP,*) '  ITYPEL=',ITYPEL,' ITYL2=',ITYL2
  745.                   MLENTI=ELEM(ITYL2)
  746.                   DO INOD=1,nbnne(itypel)
  747. *                     WRITE(IOIMP,*) '    INOD=',INOD,
  748. *     $                    ' NUM=',NUM(INOD,IEL)
  749.                      LECT(**)=NUM(INOD,IEL)
  750.                   ENDDO
  751.                   LECT(**)=ICOLOR(IEL)
  752.                ELSE
  753.                   GOTO 9999
  754.                ENDIF
  755. * Ce qu'il est dit de faire dans le tableau TABDES
  756. * On utilise la description permutée (VALP, NUMP) qui a servi a calculer
  757. * IKAS
  758. * Si la permutation avait changé l'orientation LPAIR=.FALSE.,
  759. * on change l'orientation de l'élément résultant
  760.             ELSEIF (ITYL.GE.1.AND.ITYL.LE.NTYEL) THEN
  761.                MLENTI=ELEM(ITYL)
  762.                IDX=1
  763.                ITYPG=ITYPL(ITYL)
  764.                NNO=NBNNE(ITYPG)
  765. *dbg               write(ioimp,*) 'Element ',noms(itypg),' nb no=',nno
  766.                DO INOD=1,NNO
  767.                   IDX=IDX+1
  768.                   INO1=TABDES(IDX,IKAS)
  769.                   IDX=IDX+1
  770.                   INO2=TABDES(IDX,IKAS)
  771. *dbg                  WRITE(IOIMP,*) '    INOD=',INOD,' INO1=',INO1,
  772. *dbg     $                 ' INO2=',INO2
  773.                   IF (INO2.EQ.0) THEN
  774.                      LECT(**)=NUMP(INO1)
  775.                   ELSE
  776.                      VAL1=VALP(INO1)
  777.                      VAL2=VALP(INO2)
  778.                      NUM1=NUMP(INO1)
  779.                      NUM2=NUMP(INO2)
  780. *dbg                     WRITE(IOIMP,*) '    NUM1=',NUM1,' NUM2=',NUM2
  781. *dbg                    WRITE(IOIMP,*) '    VAL1=',VAL1,' VAL2=',VAL2
  782.                      CALL ISOVA3(XISO,VAL1,VAL2,NUM1,NUM2,MLENTI,NEWNOD)
  783.                      IF (IERR.NE.0) RETURN
  784.                   ENDIF
  785.                ENDDO
  786.                LECT(**)=ICOLOR(IEL)
  787. *dbg               WRITE(IOIMP,*) 'LPAIR=',LPAIR
  788.                IF (.NOT.LPAIR) THEN
  789. *dbg                  write(ioimp,*) 'lect avant =',(lect(i),i=1,lect(/1))
  790. * Pour les éléments simplex ou qua4, il suffit d'inverser le sens de parcours
  791.                   IF ((ITYL.GE.2.AND.ITYL.LE.4).OR.(ITYL.EQ.7)) THEN
  792.                      idfin=lect(/1)-1
  793.                      iddeb=idfin-nno+1
  794.                      do iswap=1,nno/2
  795.                         idf=idfin-(iswap-1)
  796.                         idd=iddeb+(iswap-1)
  797.                         itmp=lect(idf)
  798.                         lect(idf)=lect(idd)
  799.                         lect(idd)=itmp
  800.                      enddo
  801. * Pour la pyramide, il suffit d'inverser le sens de parcours de la base carrée
  802.                   ELSEIF (ITYL.EQ.5) THEN
  803.                      idfin=lect(/1)-1-1
  804.                      iddeb=idfin-(nno-1)+1
  805.                      do iswap=1,(nno-1)/2
  806.                         idf=idfin-(iswap-1)
  807.                         idd=iddeb+(iswap-1)
  808.                         itmp=lect(idf)
  809.                         lect(idf)=lect(idd)
  810.                         lect(idd)=itmp
  811.                      enddo
  812. * Pour le prisme, il faut inverser le sens de parcours des deux bases
  813.                   ELSEIF (ITYL.EQ.6) THEN
  814.                      idf2=lect(/1)-1
  815.                      idd2=idf2-2
  816.                      idf1=idd2-1
  817.                      idd1=idf1-2
  818.                      itmp=lect(idf2)
  819.                      lect(idf2)=lect(idd2)
  820.                      lect(idd2)=itmp
  821.                      itmp=lect(idf1)
  822.                      lect(idf1)=lect(idd1)
  823.                      lect(idd1)=itmp
  824.                   ENDIF
  825. *dbg                  write(ioimp,*) 'lect apres =',(lect(i),i=1,lect(/1))
  826.                ENDIF
  827.             ELSE
  828.                GOTO 9999
  829.             ENDIF
  830.          ENDDO
  831.          goto 900
  832.  
  833. ******** Cas des CUB8 ***********************************************
  834.  114     CONTINUE
  835.          IF (IOPT.NE.0) THEN
  836. *  16 2
  837. * Type d'élément incorrect
  838.             CALL ERREUR(16)
  839.             RETURN
  840.          ENDIF
  841.  
  842. *   On remplit les tableaux MATCAS et MATARE decrivant les cas et aretes
  843. c          write(6,*) 'ISOVA2: cub8: appel a ISOVA4'
  844.          CALL ISOVA4(ITYPEL,MATCAS,NMAXCA,NMAXPT,MATARE,2,NMAXAR)
  845. c          write(6,*) 'ISOVA2: MATARE=',(MATARE(iou,1),iou=1,2)
  846. c          write(6,*) '               ',(MATARE(iou,2),iou=1,2)
  847. c          write(6,*) '                ...'
  848. c          write(6,*) 'ISOVA2: MATCAS=',(MATCAS(iou,1),iou=1,16)
  849. c          write(6,*) '               ',(MATCAS(iou,2),iou=1,16)
  850. c          write(6,*) '                ...'
  851. c          write(6,*) '               ',(MATCAS(iou,255),iou=1,16)
  852. c          write(6,*) '               ',(MATCAS(iou,256),iou=1,16)
  853.  
  854. * REM (BP) :
  855. * - on ne fait pas de test avec la tolerance (cf LEGISO autres itypel) et
  856. *   on espere que VAL(IVAL).le.XISO au lieu de VAL(IVAL).lt.XISO suffira.
  857. C ..
  858. * - on ne construit que des triangles TRI3 en se basant sur l algo :
  859. *   "the marching cubes algorithm" (domaine public)
  860.          DO 214 IEL=1,NEL
  861.             ICAS = 1
  862.             DO IVAL=1,8
  863.                VAL(IVAL)=VELCHE(MIN(N1PTEL,IVAL),MIN(N1EL,IEL))
  864. c                if(VAL(IVAL).le.XISO) ICAS=ICAS+(2**(IVAL-1))
  865.             ENDDO
  866. c           on teste les noeuds (attention: numerotation cub8 cast3m
  867. c           differente de celle du marching cubes algorithm)
  868.             if(VAL(1).le.XISO) ICAS=ICAS+1
  869.             if(VAL(4).le.XISO) ICAS=ICAS+2
  870.             if(VAL(3).le.XISO) ICAS=ICAS+4
  871.             if(VAL(2).le.XISO) ICAS=ICAS+8
  872.             if(VAL(5).le.XISO) ICAS=ICAS+16
  873.             if(VAL(8).le.XISO) ICAS=ICAS+32
  874.             if(VAL(7).le.XISO) ICAS=ICAS+64
  875.             if(VAL(6).le.XISO) ICAS=ICAS+128
  876. c             write(6,*) 'ICAS=',ICAS
  877. c             write(6,*) 'MATCAS=',(MATCAS(iou,ICAS),iou=1,16)
  878. *           Y a-t-il une isovaleur ?
  879.             IF (ICAS.gt.1.or.ICAS.lt.256) THEN
  880.                MLENTI=ELEM(3)
  881.                do iar1=1,15
  882.                   NUMAR1=MATCAS(iar1,ICAS)
  883. c                   write(6,*) 'iar1,NUMAR1=',iar1,NUMAR1
  884.                   if(NUMAR1.eq.0) goto 214
  885.                   ipoin1 = MATARE(1,NUMAR1)
  886.                   ipoin2 = MATARE(2,NUMAR1)
  887.                   VAL1=VAL(ipoin1)
  888.                   VAL2=VAL(ipoin2)
  889.                   NUM1=NUM(ipoin1,IEL)
  890.                   NUM2=NUM(ipoin2,IEL)
  891. c                   write(6,*) '  ipoin1 et 2 =',ipoin1,ipoin2,NUM1,NUM2
  892.                   CALL ISOVA3(XISO,VAL1,VAL2,NUM1,NUM2,MLENTI,NEWNOD)
  893.                   IF (IERR.NE.0) RETURN
  894. c                   write(6,*) '  VAL1 et 2=',VAL1,VAL2
  895.                   if(mod(iar1,3).eq.0) LECT(**)=ICOLOR(IEL)
  896.                enddo
  897.             ENDIF
  898.  214     CONTINUE
  899. c          write(6,*) '=> LECT=',(LECT(iou),iou=1,LECT(/1))
  900.          goto 900
  901.  
  902. ******** on termine proprement **************************************
  903.  900     CONTINUE
  904.  
  905.       ENDDO
  906. c     fin de Boucle sur les zones elementaires
  907. C***********************************************************************
  908.  
  909.       RETURN
  910. * Erreur grave
  911.  9999 CONTINUE
  912.       MOTERR(1:8)='ISOVA2  '
  913.       CALL ERREUR(1039)
  914.       RETURN
  915. *
  916. * End of subroutine ISOVA2
  917. *
  918.       END
  919.  
  920.  

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