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