Télécharger isova2.eso

Retour à la liste

Numérotation des lignes :

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

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