Télécharger delaun.eso

Retour à la liste

Numérotation des lignes :

  1. C DELAUN SOURCE BP208322 16/11/18 21:16:15 9177
  2. SUBROUTINE DELAUN(MPOVAL,IPT1,XBOI,IBOI,IPT2,LNBOIT,ITRC1,ILEA1)
  3. C----------------------------------------------------------------------C
  4. C Triangulation d'un maillage de points (POI1), 1D, 2D ou 3D, C
  5. C suivant l'agorithme de DELAUNAY donne dans P.-L. George, C
  6. C Delaunay triangulation and meshing, Hermes, 1998, p. 41. C
  7. C C
  8. C IPT1 = pointeur sur le maillage de POI1 en entree ; C
  9. C suppose ACTIF en ENTREE et en SORTIE C
  10. C XBOI = rapport entre le "diametre" du nuage de points (IPT1) et C
  11. C la taille de la boite de triangulation ; C
  12. C IBOI = indicateur permettant d'indiquer si on veut "garder" la C
  13. C boite de triangulation : C
  14. C - IBOI = 0 : on ne la garde pas (defaut) ; C
  15. C - IBOI = 1 : on la garde ; C
  16. C C
  17. C IPT2 = pointeur sur le maillage triangule en sortie C
  18. C LNBOIT = liste des noeuds formant la boite de triangulation C
  19. C ITRC1 = table contenant les centres des cercles circonscrits C
  20. C ILEA1 = table des d'adjacence de la triangulation C
  21. C C
  22. C----------------------------------------------------------------------C
  23. C Version V0 : centre et rayon cercles circons. calcules chaque fois C
  24. C Version V1 : centre et rayon cercles circons. stockes dans SEGMENT C
  25. C MCIRCONS + elements adjacents stockes dans SEGMENT C
  26. C ILEA1 C
  27. C----------------------------------------------------------------------C
  28. IMPLICIT INTEGER(I-N)
  29. IMPLICIT REAL*8 (A-H,O-Z)
  30. c
  31. -INC SMCHPOI
  32. -INC CCOPTIO
  33. -INC SMELEME
  34. -INC SMCOORD
  35. -INC CCGEOME
  36. -INC SMLENTI
  37. -INC SMTABLE
  38. PARAMETER (ICMIN=1000)
  39. PARAMETER (ZERO=1D-13)
  40. DIMENSION LNBOIT(8)
  41. DIMENSION XCMIMA(2,3),X0(4,3),XCMIL(3),XPI1(3),XCBJ(3),XG0(3)
  42. LOGICAL BMOD,BHORS,BFF
  43. C
  44. POINTEUR LCAVI.MLENTI,LBOUL.MLENTI,LBASE.MLENTI,LREF1.MLENTI
  45. POINTEUR LEADJ.MLENTI,LNVADJ.MLENTI,LSUP.MLENTI,LSCP.MLENTI
  46. C
  47. SEGMENT,MCIRCONS
  48. REAL*8 TRC(NBE1,4)
  49. ENDSEGMENT
  50. POINTEUR ITRC1.MCIRCONS
  51. C
  52. SEGMENT,MADJACEN
  53. INTEGER LEAC(NBL1,IDIM+1,2)
  54. ENDSEGMENT
  55. POINTEUR ILEA1.MADJACEN
  56. C
  57. SEGMENT,MLIAIS
  58. INTEGER LIAS(NB1,NB1)
  59. ENDSEGMENT
  60. POINTEUR ITL1.MLIAIS
  61. C
  62. SEGMENT RAYONOE
  63. REAL*8 RAYN(NBNER)
  64. ENDSEGMENT
  65. POINTEUR IRN1.RAYONOE
  66. c
  67. C----------------------------------------------------------------------C
  68. C PARTIE 1 C
  69. C Construction de la boite englobant tous les points C
  70. C----------------------------------------------------------------------C
  71. C
  72. C---- 1.1 Je determine la taille du maillage
  73. SEGACT,MCOORD
  74. IDIMP1=IDIM+1
  75. IPT2=IPT1
  76. IF ((IBOI.EQ.1).AND.(IPT1.LISREF(/1).NE.0)) THEN
  77. IPT2=IPT1.LISREF(1)
  78. SEGACT,IPT2
  79. ENDIF
  80. IREF1=(IPT1.NUM(1,1)-1)*IDIMP1
  81. DDMAX=0.
  82. DO I=1,IDIM
  83. XCMIMA(1,I)=XCOOR(IREF1+I)
  84. XCMIMA(2,I)=XCOOR(IREF1+I)
  85. DO J=2,IPT1.NUM(/2)
  86. IREFJ=(IPT1.NUM(1,J)-1)*IDIMP1
  87. XCNJ1=XCOOR(IREFJ+I)
  88. IF (XCNJ1.LT.XCMIMA(1,I)) XCMIMA(1,I)=XCNJ1
  89. IF (XCNJ1.GT.XCMIMA(2,I)) XCMIMA(2,I)=XCNJ1
  90. ENDDO
  91. IF ((IBOI.EQ.1).AND.(IPT1.LISREF(/1).NE.0)) THEN
  92. DO J=1,IPT2.NUM(/2)
  93. IREFJ=(IPT2.NUM(1,J)-1)*IDIMP1
  94. XCNJ1=XCOOR(IREFJ+I)
  95. IF (XCNJ1.LT.XCMIMA(1,I)) XCMIMA(1,I)=XCNJ1
  96. IF (XCNJ1.GT.XCMIMA(2,I)) XCMIMA(2,I)=XCNJ1
  97. ENDDO
  98. ENDIF
  99. DDMAXI=0.5*(XCMIMA(2,I)-XCMIMA(1,I))
  100. IF (DDMAXI.GT.DDMAX) DDMAX=DDMAXI
  101. ENDDO
  102. C
  103. C---- 1.2 On construit la boite
  104. C 1.2.1 Points coins de la boite
  105. DBOIT=XBOI*DDMAX
  106. NBPTS0=XCOOR(/1)/IDIMP1
  107. NBPTS=NBPTS0+(2**IDIM)
  108. SEGADJ,MCOORD
  109. IF (IDIM.EQ.1) THEN
  110. XCMIL(1)=0.5*(XCMIMA(1,1)+XCMIMA(2,1))
  111. IREF1=(NBPTS-2)*IDIMP1
  112. IREF2=(NBPTS-1)*IDIMP1
  113. XCOOR(IREF1+1)=XCMIL(1)-DBOIT
  114. XCOOR(IREF2+1)=XCMIL(1)+DBOIT
  115. XCOOR(IREF1+2)=0.
  116. XCOOR(IREF2+2)=0.
  117. ELSEIF (IDIM.EQ.2) THEN
  118. DO I=1,IDIM
  119. XCMIL(I)=0.5*(XCMIMA(1,I)+XCMIMA(2,I))
  120. ENDDO
  121. IREF1=(NBPTS-4)*IDIMP1
  122. IREF2=(NBPTS-3)*IDIMP1
  123. IREF3=(NBPTS-2)*IDIMP1
  124. IREF4=(NBPTS-1)*IDIMP1
  125. XCOOR(IREF1+1)=XCMIL(1)-DBOIT
  126. XCOOR(IREF2+1)=XCMIL(1)+DBOIT
  127. XCOOR(IREF3+1)=XCMIL(1)+DBOIT
  128. XCOOR(IREF4+1)=XCMIL(1)-DBOIT
  129. XCOOR(IREF1+2)=XCMIL(2)-DBOIT
  130. XCOOR(IREF2+2)=XCMIL(2)-DBOIT
  131. XCOOR(IREF3+2)=XCMIL(2)+DBOIT
  132. XCOOR(IREF4+2)=XCMIL(2)+DBOIT
  133. XCOOR(IREF1+3)=0.
  134. XCOOR(IREF2+3)=0.
  135. XCOOR(IREF3+3)=0.
  136. XCOOR(IREF4+3)=0.
  137. ELSEIF (IDIM.EQ.3) THEN
  138. DO I=1,IDIM
  139. XCMIL(I)=0.5*(XCMIMA(1,I)+XCMIMA(2,I))
  140. ENDDO
  141. IREF1=(NBPTS-8)*IDIMP1
  142. IREF2=(NBPTS-7)*IDIMP1
  143. IREF3=(NBPTS-6)*IDIMP1
  144. IREF4=(NBPTS-5)*IDIMP1
  145. IREF5=(NBPTS-4)*IDIMP1
  146. IREF6=(NBPTS-3)*IDIMP1
  147. IREF7=(NBPTS-2)*IDIMP1
  148. IREF8=(NBPTS-1)*IDIMP1
  149. XCOOR(IREF1+1)=XCMIL(1)-DBOIT
  150. XCOOR(IREF2+1)=XCMIL(1)+DBOIT
  151. XCOOR(IREF3+1)=XCMIL(1)+DBOIT
  152. XCOOR(IREF4+1)=XCMIL(1)-DBOIT
  153. XCOOR(IREF5+1)=XCMIL(1)-DBOIT
  154. XCOOR(IREF6+1)=XCMIL(1)+DBOIT
  155. XCOOR(IREF7+1)=XCMIL(1)+DBOIT
  156. XCOOR(IREF8+1)=XCMIL(1)-DBOIT
  157. XCOOR(IREF1+2)=XCMIL(2)-DBOIT
  158. XCOOR(IREF2+2)=XCMIL(2)-DBOIT
  159. XCOOR(IREF3+2)=XCMIL(2)+DBOIT
  160. XCOOR(IREF4+2)=XCMIL(2)+DBOIT
  161. XCOOR(IREF5+2)=XCMIL(2)-DBOIT
  162. XCOOR(IREF6+2)=XCMIL(2)-DBOIT
  163. XCOOR(IREF7+2)=XCMIL(2)+DBOIT
  164. XCOOR(IREF8+2)=XCMIL(2)+DBOIT
  165. XCOOR(IREF1+3)=XCMIL(3)-DBOIT
  166. XCOOR(IREF2+3)=XCMIL(3)-DBOIT
  167. XCOOR(IREF3+3)=XCMIL(3)-DBOIT
  168. XCOOR(IREF4+3)=XCMIL(3)-DBOIT
  169. XCOOR(IREF5+3)=XCMIL(3)+DBOIT
  170. XCOOR(IREF6+3)=XCMIL(3)+DBOIT
  171. XCOOR(IREF7+3)=XCMIL(3)+DBOIT
  172. XCOOR(IREF8+3)=XCMIL(3)+DBOIT
  173. XCOOR(IREF1+4)=0.
  174. XCOOR(IREF2+4)=0.
  175. XCOOR(IREF3+4)=0.
  176. XCOOR(IREF4+4)=0.
  177. XCOOR(IREF5+4)=0.
  178. XCOOR(IREF6+4)=0.
  179. XCOOR(IREF7+4)=0.
  180. XCOOR(IREF8+4)=0.
  181. ENDIF
  182. C
  183. C 1.2.2 Maillage de la boite
  184. NBSOUS=0
  185. NBREF=0
  186. NNOE1=NBPTS-(2**IDIM)+1
  187. IF (IDIM.EQ.1) THEN
  188. C Creation des elements SEG2 en definissant les noeuds
  189. NBNN=2
  190. NBELEM=IPT1.NUM(/2)+1
  191. SEGINI,IPT2
  192. IPT2.ITYPEL=2
  193. NBELE2=1
  194. IPT2.NUM(1,1)=NNOE1
  195. IPT2.NUM(2,1)=NNOE1+1
  196. C Liste des numeros des noeuds de la boite
  197. NNBOIT=2
  198. LNBOIT(1)=NNOE1
  199. LNBOIT(2)=NNOE1+1
  200. C Initialisation du segment ILEA1 contenant la table de
  201. C connectivite des elements de IPT2
  202. NBL1=NBELEM
  203. SEGINI,ILEA1
  204. ILEA1.LEAC(1,1,1)= 0
  205. ILEA1.LEAC(1,1,2)= 0
  206. ILEA1.LEAC(1,2,1)= 0
  207. ILEA1.LEAC(1,2,2)= 0
  208. ELSEIF (IDIM.EQ.2) THEN
  209. C Creation des elements TRI3 en definissant les noeuds
  210. NBNN=3
  211. NBELEM=IPT1.NUM(/2)*3
  212. SEGINI,IPT2
  213. IPT2.ITYPEL=4
  214. NBELE2=2
  215. IPT2.NUM(1,1)=NNOE1
  216. IPT2.NUM(2,1)=NNOE1+1
  217. IPT2.NUM(3,1)=NNOE1+3
  218. IPT2.NUM(1,2)=NNOE1+1
  219. IPT2.NUM(2,2)=NNOE1+2
  220. IPT2.NUM(3,2)=NNOE1+3
  221. C Liste des numeros des noeuds de la boite
  222. NNBOIT=4
  223. LNBOIT(1)=NNOE1
  224. LNBOIT(2)=NNOE1+1
  225. LNBOIT(3)=NNOE1+2
  226. LNBOIT(4)=NNOE1+3
  227. C Initialisation du segment ILEA1 contenant la table de
  228. C connectivite des elements de IPT2
  229. NBL1=NBELEM
  230. SEGINI,ILEA1
  231. ILEA1.LEAC(1,1,1)= 2
  232. ILEA1.LEAC(1,1,2)= 2
  233. ILEA1.LEAC(1,2,1)= 0
  234. ILEA1.LEAC(1,2,2)= 0
  235. ILEA1.LEAC(1,3,1)= 0
  236. ILEA1.LEAC(1,3,2)= 0
  237. ILEA1.LEAC(2,1,1)= 0
  238. ILEA1.LEAC(2,1,2)= 0
  239. ILEA1.LEAC(2,2,1)= 1
  240. ILEA1.LEAC(2,2,2)= 1
  241. ILEA1.LEAC(2,3,1)= 0
  242. ILEA1.LEAC(2,3,2)= 0
  243. ELSEIF (IDIM.EQ.3) THEN
  244. C Creation des elements TET4 en definissant les noeuds
  245. NBNN=4
  246. NBELEM=IPT1.NUM(/2)*7
  247. SEGINI,IPT2
  248. IPT2.ITYPEL=23
  249. NBELE2=5
  250. IPT2.NUM(1,1)=NNOE1
  251. IPT2.NUM(2,1)=NNOE1+1
  252. IPT2.NUM(3,1)=NNOE1+3
  253. IPT2.NUM(4,1)=NNOE1+4
  254. IPT2.NUM(1,2)=NNOE1+1
  255. IPT2.NUM(2,2)=NNOE1+2
  256. IPT2.NUM(3,2)=NNOE1+3
  257. IPT2.NUM(4,2)=NNOE1+6
  258. IPT2.NUM(1,3)=NNOE1+1
  259. IPT2.NUM(2,3)=NNOE1+4
  260. IPT2.NUM(3,3)=NNOE1+5
  261. IPT2.NUM(4,3)=NNOE1+6
  262. IPT2.NUM(1,4)=NNOE1+3
  263. IPT2.NUM(2,4)=NNOE1+4
  264. IPT2.NUM(3,4)=NNOE1+6
  265. IPT2.NUM(4,4)=NNOE1+7
  266. IPT2.NUM(1,5)=NNOE1+1
  267. IPT2.NUM(2,5)=NNOE1+3
  268. IPT2.NUM(3,5)=NNOE1+4
  269. IPT2.NUM(4,5)=NNOE1+6
  270. C Liste des numeros des noeuds de la boite
  271. NNBOIT=8
  272. LNBOIT(1)=NNOE1
  273. LNBOIT(2)=NNOE1+1
  274. LNBOIT(3)=NNOE1+2
  275. LNBOIT(4)=NNOE1+3
  276. LNBOIT(5)=NNOE1+4
  277. LNBOIT(6)=NNOE1+5
  278. LNBOIT(7)=NNOE1+6
  279. LNBOIT(8)=NNOE1+7
  280. C Initialisation du segment ILEA1 contenant la table de
  281. C connectivite des elements de IPT2
  282. NBL1=NBELEM
  283. SEGINI,ILEA1
  284. ILEA1.LEAC(1,1,1)= 5
  285. ILEA1.LEAC(1,1,2)= 4
  286. ILEA1.LEAC(1,2,1)= 0
  287. ILEA1.LEAC(1,2,2)= 0
  288. ILEA1.LEAC(1,3,1)= 0
  289. ILEA1.LEAC(1,3,2)= 0
  290. ILEA1.LEAC(1,4,1)= 0
  291. ILEA1.LEAC(1,4,2)= 0
  292. ILEA1.LEAC(2,1,1)= 0
  293. ILEA1.LEAC(2,1,2)= 0
  294. ILEA1.LEAC(2,2,1)= 5
  295. ILEA1.LEAC(2,2,2)= 3
  296. ILEA1.LEAC(2,3,1)= 0
  297. ILEA1.LEAC(2,3,2)= 0
  298. ILEA1.LEAC(2,4,1)= 0
  299. ILEA1.LEAC(2,4,2)= 0
  300. ILEA1.LEAC(3,1,1)= 0
  301. ILEA1.LEAC(3,1,2)= 0
  302. ILEA1.LEAC(3,2,1)= 0
  303. ILEA1.LEAC(3,2,2)= 0
  304. ILEA1.LEAC(3,3,1)= 5
  305. ILEA1.LEAC(3,3,2)= 2
  306. ILEA1.LEAC(3,4,1)= 0
  307. ILEA1.LEAC(3,4,2)= 0
  308. ILEA1.LEAC(4,1,1)= 0
  309. ILEA1.LEAC(4,1,2)= 0
  310. ILEA1.LEAC(4,2,1)= 0
  311. ILEA1.LEAC(4,2,2)= 0
  312. ILEA1.LEAC(4,3,1)= 0
  313. ILEA1.LEAC(4,3,2)= 0
  314. ILEA1.LEAC(4,4,1)= 5
  315. ILEA1.LEAC(4,4,2)= 1
  316. ILEA1.LEAC(5,1,1)= 4
  317. ILEA1.LEAC(5,1,2)= 4
  318. ILEA1.LEAC(5,2,1)= 3
  319. ILEA1.LEAC(5,2,2)= 3
  320. ILEA1.LEAC(5,3,1)= 2
  321. ILEA1.LEAC(5,3,2)= 2
  322. ILEA1.LEAC(5,4,1)= 1
  323. ILEA1.LEAC(5,4,2)= 1
  324. ENDIF
  325. C
  326. C---- 1.3 Initialisation du segment contenant les rayons et centres
  327. C de la triangulation :
  328. IF (IDIM.NE.1) THEN
  329. NBE1=IPT2.NUM(/2)
  330. SEGINI,ITRC1
  331. NBNER = LNBOIT(NNBOIT)
  332. SEGINI,IRN1
  333. C
  334. if (MPOVAL .ne. 0) then
  335. DO I=1,IPT1.NUM(/2)
  336. IRN1.RAYN(IPT1.NUM(1,I))= MPOVAL.VPOCHA(I,1)
  337. END DO
  338. endif
  339. C
  340. DO I0=1,NBELE2
  341. C BCIRCO renvoie le centre et le rayon de l element considere
  342. C
  343. CALL BCIRCO(IRN1,IPT2,I0,XCBJ,XRBJ)
  344. C
  345. ITRC1.TRC(I0,1)=XRBJ
  346. DO J0=1,IDIM
  347. ITRC1.TRC(I0,J0+1)=XCBJ(J0)
  348. ENDDO
  349. ENDDO
  350. ENDIF
  351. C
  352. C
  353. C
  354. C
  355. C----------------------------------------------------------------------C
  356. C PARTIE 2 C
  357. C Triangulation incrementale C
  358. C----------------------------------------------------------------------C
  359. C
  360. JG=1
  361. SEGINI,LSUP
  362. C On boucle sur tous les points a trianguler
  363. NBSOUS=0
  364. NBREF=0
  365. NBNN=IPT2.NUM(/1)
  366. DO I=1,IPT1.NUM(/2)
  367. C WRITE(6,*)
  368. C WRITE(6,*) 'Iteration I=',I,' /',IPT1.NUM(/2)
  369. C
  370. C------- 2.1 Construction de la base
  371. C
  372. C Je tire un point a ajouter a la triangulation
  373. NPI1=IPT1.NUM(1,I)
  374. C XPI1 : vecteur contenant les coordonnees du point
  375. DO J=1,IDIM
  376. XPI1(J)=XCOOR((NPI1-1)*IDIMP1+J)
  377. ENDDO
  378. C
  379. C 2.1.1 Recherche d'un element de la base (NBASE1)
  380. C Je demarre avec le dernier elem. de IPT2 (NEL0) et cherche si
  381. C NPI1 est dedans. Si oui, c'est gagne, sinon :
  382. C je determine le voisin de NEL0 oriente vers NPI1
  383. C => nouveau NEL0, puis je recommence jusqu'a y arriver
  384. NEL0=NBELE2
  385. NBASE=1
  386. C Cas particulier de la dimension 1
  387. IF (IDIM.EQ.1) THEN
  388. DO J=1,NBELE2
  389. C Coordonnees des noeuds A et B de l'element NEL0
  390. C A est a gauche et B est a droite
  391. NREFA=IPT2.NUM(1,NEL0)
  392. NREFB=IPT2.NUM(2,NEL0)
  393. XA=XCOOR((NREFA-1)*IDIMP1+1)
  394. XB=XCOOR((NREFB-1)*IDIMP1+1)
  395. C Cas ou NPI1 est confondu avec A ou B
  396. X1=ABS(XPI1(1))
  397. IF (X1.LT.ZERO) X1=1.
  398. DTEST1=ABS(XPI1(1)-XA)/X1
  399. DTEST2=ABS(XPI1(1)-XB)/X1
  400. IF ((DTEST1.LT.ZERO).OR.(DTEST2.LT.ZERO)) THEN
  401. NBNN=IPT2.NUM(/1)
  402. NBELEM=NBELE2
  403. NBSOUS=0
  404. NBREF=0
  405. SEGADJ,IPT2
  406. C Erreur de triangulation : vérifiez qu'il n'y a pas de
  407. C points confondus
  408. CALL ERREUR(846)
  409. GOTO 2501
  410. C Cas ou le NPT1 est entre A et B, on a trouve la base
  411. ELSEIF ((XPI1(1).GT.XA).AND.(XPI1(1).LT.XB)) THEN
  412. NBASE1=NEL0
  413. GOTO 2381
  414. C Si NPI1 est a gauche de A, on va chercher le voisin
  415. ELSEIF (XPI1(1).LT.XA) THEN
  416. NEL0=ILEA1.LEAC(NEL0,2,1)
  417. C Si NPI1 est a droite de B, on va chercher le voisin
  418. ELSEIF (XPI1(1).GT.XB) THEN
  419. NEL0=ILEA1.LEAC(NEL0,1,1)
  420. ENDIF
  421. ENDDO
  422. C Si l'on passe ici, c'est que l'on a pas trouve NBASE1 !
  423. NBNN=IPT2.NUM(/1)
  424. NBELEM=NBELE2
  425. NBSOUS=0
  426. NBREF=0
  427. SEGADJ,IPT2
  428. C Erreur de triangulation : vérifiez qu'il n'y a pas de
  429. C points confondus
  430. CALL ERREUR(846)
  431. GOTO 2501
  432. ENDIF
  433. C Cas generaux en dimension 2 et 3
  434. JG=1
  435. SEGINI,LSCP
  436. NEL00=NEL0
  437. NEL000=NEL0
  438. DO J=1,NBELE2
  439. NCP=0
  440. C Coordonnees des noeuds de l'element NEL0
  441. DO K=1,IDIMP1
  442. NREF0=IPT2.NUM(K,NEL0)
  443. DO KK=1,IDIM
  444. X0(K,KK)=XCOOR((NREF0-1)*IDIMP1+KK)
  445. ENDDO
  446. ENDDO
  447. C Coordonnees du barycentre G de NEL0
  448. DO KK=1,IDIM
  449. XSOM=0.
  450. DO K=1,IDIMP1
  451. XSOM=XSOM+X0(K,KK)
  452. ENDDO
  453. XG0(KK)=XSOM/(IDIMP1)
  454. ENDDO
  455. IK1=0
  456. IK2=0
  457. IK3=0
  458. C En dimension 2
  459. IF (IDIM.EQ.2) THEN
  460. C On effectue le meme test sur chaque cote
  461. DO K=1,3
  462. IF (K.EQ.1) THEN
  463. NA=2
  464. NB=3
  465. ELSEIF (K.EQ.2) THEN
  466. NA=1
  467. NB=3
  468. ELSEIF (K.EQ.3) THEN
  469. NA=1
  470. NB=2
  471. ENDIF
  472. C Normale au cote AB du triangle
  473. XN=X0(NA,2)-X0(NB,2)
  474. YN=X0(NB,1)-X0(NA,1)
  475. XXN=SQRT((XN**2)+(YN**2))
  476. XN=XN/XXN
  477. YN=YN/XXN
  478. C On teste si NPI1 est du meme cote que G par rapport a
  479. C la droite (AB)
  480. PS1=(XN*(XPI1(1)-X0(NA,1)))+(YN*(XPI1(2)-X0(NA,2)))
  481. XX0=SQRT(((XPI1(1)-X0(NA,1))**2)+
  482. & ((XPI1(2)-X0(NA,2))**2))
  483. PS11=PS1/XX0
  484. IF (ABS(PS11).LT.ZERO) THEN
  485. NSCP=K
  486. NCP=1
  487. NBASE=2
  488. GOTO 2111
  489. ENDIF
  490. PS2=(XN*(XG0(1)-X0(NA,1)))+(YN*(XG0(2)-X0(NA,2)))
  491. IF ((SIGN(1D0,PS1)).EQ.(SIGN(1D0,PS2))) THEN
  492. GOTO 2111
  493. ELSE
  494. NADJ0=K
  495. GOTO 2112
  496. ENDIF
  497. 2111 CONTINUE
  498. IF (K.EQ.1) IK1=1
  499. IF (K.EQ.2) IK2=1
  500. IF (K.EQ.3) IK3=1
  501. ENDDO
  502. IF ((IK1.EQ.1).AND.(IK2.EQ.1).AND.(IK3.EQ.1)) THEN
  503. C NPI1 est dans NEL0, on a trouve un element de la base
  504. NBASE1=NEL0
  505. GOTO 2115
  506. ENDIF
  507. 2112 CONTINUE
  508. C NPI1 et G ne sont pas du meme cote par rapport au cote
  509. C vu par le noeud NADJ0 de NEL0 => on passe a l'element
  510. C voisin NEL1
  511. IF (J.GE.3) NEL000=NEL00
  512. IF (J.GE.2) NEL00=NEL0
  513. NEL0=ILEA1.LEAC(NEL0,NADJ0,1)
  514. IF (NEL0.EQ.NEL00) THEN
  515. C PRINT*,'Erreur dans la table des elements adjacents'
  516. C PRINT*,'L''element',NEL0,'est voisin de lui meme !'
  517. CALL ERREUR(846)
  518. GOTO 2501
  519. ENDIF
  520. IF (NEL0.EQ.NEL000) THEN
  521. C PRINT*,'Erreur dans la construction du maillage'
  522. C PRINT*,'Les elements',NEL0,NEL00,'se chevauchent !'
  523. CALL ERREUR(846)
  524. GOTO 2501
  525. ENDIF
  526. ENDIF
  527. C En dimension 3
  528. IF (IDIM.EQ.3) THEN
  529. IK4=0
  530. C On effectue le test sur chaque face
  531. DO K=1,4
  532. IF (K.EQ.1) THEN
  533. NA=2
  534. NB=3
  535. NC=4
  536. ELSEIF (K.EQ.2) THEN
  537. NA=1
  538. NB=3
  539. NC=4
  540. ELSEIF (K.EQ.3) THEN
  541. NA=1
  542. NB=2
  543. NC=4
  544. ELSEIF (K.EQ.4) THEN
  545. NA=1
  546. NB=2
  547. NC=3
  548. ENDIF
  549. C Normale au plan ABC du triangle (n = AB^AC)
  550. XN=((X0(NB,2)-X0(NA,2))*(X0(NC,3)-X0(NA,3))) -
  551. & ((X0(NB,3)-X0(NA,3))*(X0(NC,2)-X0(NA,2)))
  552. YN=6#40;4/3pan>(X0(NB,3)-X0(NA,3))*(X0(NC,1)-X0(NA,1))) -
  553. & ((X0(NB,1)-X0(NA,1))*(X0(NC,3)-X0(NA,3)))
  554. ZN=((X0(NB,1)-X0(NA,1))*(X0(NC,2)-X0(NA,2))) -
  555. & ((X0(NB,2)-X0(NA,2))*(X0(NC,1)-X0(NA,1)))
  556. XXN=SQRT((XN**2)+(YN**2)+(ZN**2))
  557. XN=XN/XXN
  558. YN=YN/XXN
  559. ZN=ZN/XXN
  560. C On teste si NPI1 est du meme cote que G par rapport au
  561. C plan ABC
  562. PS1=(XN*(XPI1(1)-X0(NA,1))) +
  563. & (YN*(XPI1(2)-X0(NA,2))) +
  564. & (ZN*(XPI1(3)-X0(NA,3)))
  565. XX0=SQRT(((XPI1(1)-X0(NA,1))**2) +
  566. & ((XPI1(2)-X0(NA,2))**2)+
  567. & ((XPI1(3)-X0(NA,3))**2))
  568. PS11=PS1/XX0
  569. IF (ABS(PS11).LT.ZERO) THEN
  570. NCP=NCP+1
  571. JG=NCP
  572. SEGADJ,LSCP
  573. LSCP.LECT(NCP)=K
  574. GOTO 2113
  575. ENDIF
  576. PS2=(XN*(XG0(1)-X0(NA,1))) +
  577. & (YN*(XG0(2)-X0(NA,2))) +
  578. & (ZN*(XG0(3)-X0(NA,3)))
  579. IF ((SIGN(1D0,PS1)).EQ.(SIGN(1D0,PS2))) THEN
  580. GOTO 2113
  581. ELSE
  582. GOTO 2114
  583. ENDIF
  584. 2113 CONTINUE
  585. IF (K.EQ.1) IK1=1
  586. IF (K.EQ.2) IK2=1
  587. IF (K.EQ.3) IK3=1
  588. IF (K.EQ.4) IK4=1
  589. ENDDO
  590. IF ((IK1.EQ.1).AND.(IK2.EQ.1).AND.
  591. & (IK3.EQ.1).AND.(IK4.EQ.1)) THEN
  592. C NPI1 est bien dans NEL0, on a donc trouve la base
  593. NBASE1=NEL0
  594. GOTO 2115
  595. ENDIF
  596. 2114 CONTINUE
  597. C NPI1 et G ne sont pas du meme cote par rapport au plan
  598. C vu par le noeud NADJ0 de NEL0 => on passe a l'element
  599. C adjacent a NEL0 vu par le noeud K
  600. IF (J.GE.3) NEL000=NEL00
  601. IF (J.GE.2) NEL00=NEL0
  602. NEL0=ILEA1.LEAC(NEL0,K,1)
  603. IF (NEL0.EQ.NEL00) THEN
  604. C PRINT*,'Erreur dans la table des elements adjacents'
  605. C PRINT*,'L''element',NEL0,'est voisin de lui meme !'
  606. CALL ERREUR(846)
  607. GOTO 2501
  608. ENDIF
  609. IF (NEL0.EQ.NEL000) THEN
  610. C PRINT*,'Erreur dans la construction du maillage'
  611. C PRINT*,'Les elements',NEL0,NEL00,'se chevauchent !'
  612. CALL ERREUR(846)
  613. GOTO 2501
  614. ENDIF
  615. ENDIF
  616. ENDDO
  617. C Si l'on passe ici, c'est que l'on a pas trouve NBASE1 !
  618. NBNN=IPT2.NUM(/1)
  619. NBELEM=NBELE2
  620. NBSOUS=0
  621. NBREF=0
  622. SEGADJ,IPT2
  623. C Erreur de triangulation : vérifiez qu'il n'y a pas de points
  624. C confondus
  625. CALL ERREUR(846)
  626. GOTO 2501
  627. 2115 CONTINUE
  628. C Test si NPI1 n'est pas confondu avec un des noeuds de NBASE1
  629. DO J=1,IDIMP1
  630. NREFJ=IPT2.NUM(J,NBASE1)
  631. NC=0
  632. DO JJ=1,IDIM
  633. X1=ABS(XPI1(JJ))
  634. IF (X1.LT.ZERO) X1=1.
  635. DTEST=ABS(XPI1(JJ)-XCOOR((NREFJ-1)*IDIMP1+JJ))/X1
  636. IF (DTEST.LT.ZERO) NC=NC+1
  637. ENDDO
  638. IF (NC.EQ.IDIM) THEN
  639. NBNN=IPT2.NUM(/1)
  640. NBELEM=NBELE2
  641. NBSOUS=0
  642. NBREF=0
  643. SEGADJ,IPT2
  644. C Erreur de triangulation : vérifiez qu'il n'y a pas de
  645. C points confondus
  646. CALL ERREUR(846)
  647. GOTO 2501
  648. ENDIF
  649. ENDDO
  650. C
  651. C 2.1.2 Ajout eventuel d'autres elements a la base
  652. C LBASE : liste des numeros des elements de IPT2 appartenant a
  653. C la base, taille initialisee a NBASE
  654. JG=NBASE
  655. SEGINI,LBASE
  656. LBASE.LECT(1)=NBASE1
  657. C Ajout des elments supplementaires a LBASE si necessaire
  658. IF (NCP.NE.0) THEN
  659. C En dimension 2, un seul adjacent est ajoute a LBASE, celui
  660. C partageant le cote
  661. IF (IDIM.EQ.2) THEN
  662. LBASE.LECT(2)=ILEA1.LEAC(NBASE1,NSCP,1)
  663. C En dimension 3 plusieurs cas sont possibles
  664. ELSEIF (IDIM.EQ.3) THEN
  665. C Cas ou NPI1 est sur un des 4 plans
  666. IF (NCP.EQ.1) THEN
  667. C On recupere l'adjacent partageant le plan
  668. NBASE=2
  669. JG=NBASE
  670. SEGADJ,LBASE
  671. LBASE.LECT(2)=ILEA1.LEAC(NBASE1,LSCP.LECT(1),1)
  672. C Cas ou NPI1 est sur une des aretes
  673. ELSEIF (NCP.EQ.2) THEN
  674. C On recupere tous les elements partageant l'arete
  675. NA=1
  676. NB=2
  677. NC=LSCP.LECT(1)
  678. ND=LSCP.LECT(2)
  679. NREFC=IPT2.NUM(NC,NBASE1)
  680. NREFD=IPT2.NUM(ND,NBASE1)
  681. DO L=1,4
  682. IF ((NA.EQ.NC).OR.(NA.EQ.ND).OR.(NA.EQ.NB)) NA=NA+1
  683. IF ((NB.EQ.NC).OR.(NB.EQ.ND).OR.(NB.EQ.NA)) NB=NB+1
  684. ENDDO
  685. NREFA=IPT2.NUM(NA,NBASE1)
  686. NREFB=IPT2.NUM(NB,NBASE1)
  687. N1=NBASE1
  688. DO L=1,100
  689. N2=ILEA1.LEAC(N1,NC,1)
  690. ND=ILEA1.LEAC(N1,NC,2)
  691. IF (N2.EQ.NBASE1) GOTO 2121
  692. NBASE=NBASE+1
  693. JG=NBASE
  694. SEGADJ,LBASE
  695. LBASE.LECT(NBASE)=N2
  696. DO M=1,4
  697. IF (IPT2.NUM(M,N2).EQ.NREFA) NA=M
  698. IF (IPT2.NUM(M,N2).EQ.NREFB) NB=M
  699. IF (IPT2.NUM(M,N2).EQ.NREFD) NC=M
  700. ENDDO
  701. NREFD=IPT2.NUM(ND,N2)
  702. N1=N2
  703. ENDDO
  704. 2121 CONTINUE
  705. C Cas ou NPI1 est confondu avec un des sommets
  706. ELSEIF (NCP.EQ.3) THEN
  707. NBNN=IPT2.NUM(/1)
  708. NBELEM=NBELE2
  709. NBSOUS=0
  710. NBREF=0
  711. SEGADJ,IPT2
  712. C Erreur de triangulation : vérifiez qu'il n'y a pas de
  713. C points confondus
  714. CALL ERREUR(846)
  715. GOTO 2501
  716. ENDIF
  717. ENDIF
  718. ENDIF
  719. C
  720. C------- 2.2 Construction de la cavite
  721. C
  722. C 2.2.1 Construction de la cavite (MELEME IPT3)
  723. C par adjacence a partir de l'element NBASE1 de la base
  724. C
  725. C LCAVI : liste des numeros des elements de IPT2 appartenant a
  726. C la cavite, taille initialisee a NBELE2
  727. JG=NBELE2
  728. SEGINI,LCAVI
  729. C Le premier element de la cavite est NBASE1
  730. LCAVI.LECT(1)=NBASE1
  731. C NCAVI : nombre d'elements dans la cavite (toujours a jour !)
  732. NCAVI=1
  733. C Initialisation de la cavite IPT3 a NBELE2
  734. NBELEM=NBELE2
  735. NBNN=IPT2.NUM(/1)
  736. SEGINI,IPT3
  737. IPT3.ITYPEL=IPT2.ITYPEL
  738. DO K=1,IPT3.NUM(/1)
  739. IPT3.NUM(K,NCAVI)=IPT2.NUM(K,NBASE1)
  740. ENDDO
  741. C Parcours des adjacents de maniere recursive jusqu'a stabilite
  742. ITRA=1
  743. DO J=1,NBELE2
  744. C Element a traiter de la liste LCAVI
  745. NEJ=LCAVI.LECT(ITRA)
  746. IF ((NEJ.EQ.0).OR.(J.EQ.NBELE2)) GOTO 2213
  747. C On prend les elements adjacents de NEJ
  748. DO K=1,IDIMP1
  749. NEAJ=ILEA1.LEAC(NEJ,K,1)
  750. C On verifie qu'il y a bien un element adjacent
  751. IF (NEAJ.EQ.0) GOTO 2212
  752. C On verifie que NEAJ n'est pas deja dans LCAVI
  753. DO L=1,LCAVI.LECT(/1)
  754. IF (LCAVI.LECT(L).EQ.0) GOTO 2211
  755. IF (LCAVI.LECT(L).EQ.NEAJ) GOTO 2212
  756. ENDDO
  757. 2211 CONTINUE
  758. C XRBJ : rayon de la boule circonscrite a l'elem. NEAJ
  759. C XCBJ(i) : ieme coordonnee de son centre
  760. XRBJ=ITRC1.TRC(NEAJ,1)
  761. DO L=1,IDIM
  762. XCBJ(L)=ITRC1.TRC(NEAJ,L+1)
  763. ENDDO
  764. C DIJ1 : distance de NPI1 au centre de la boule de
  765. C l'element considere
  766. DIJ1=0.
  767. DO L=1,IDIM
  768. DIJ1=DIJ1+((XPI1(L)-XCBJ(L))**2)
  769. ENDDO
  770. C
  771. IF (MPOVAL .NE. 0)
  772. & DIJ1 = DIJ1 - ((IRN1.RAYN(IPT1.NUM(1,I)))**2)
  773. C
  774. DIJ1=(DIJ1/(XRBJ**2))-1.
  775. C Si NEAJ est un nouvel element de la cavite,
  776. IF (DIJ1.LT.ZERO) THEN
  777. C alors on l'ajoute a LCAVI et on incremente NCAVI
  778. NCAVI=NCAVI+1
  779. LCAVI.LECT(NCAVI)=NEAJ
  780. C et on remplis IPT3 avec ce nouvel element
  781. DO L=1,IPT3.NUM(/1)
  782. IPT3.NUM(L,NCAVI)=IPT2.NUM(L,NEAJ)
  783. ENDDO
  784. ENDIF
  785. 2212 CONTINUE
  786. ENDDO
  787. ITRA=ITRA+1
  788. ENDDO
  789. 2213 CONTINUE
  790. C Ajustement de LCAVI et IPT3 a la bonne taille si besoin
  791. IF (LCAVI.LECT(/1).GT.NCAVI) THEN
  792. JG=NCAVI
  793. SEGADJ,LCAVI
  794. NBNN=IPT3.NUM(/1)
  795. NBELEM=NCAVI
  796. SEGADJ,IPT3
  797. ENDIF
  798. NBELE3=IPT3.NUM(/2)
  799. C
  800. C 2.2.2 Verification : la cavite doit etre connexe, chaque
  801. C element doit avoir au moins un element adjacent, c'est a dire
  802. C partageant un cote/face. De plus, un cote/face ne peut etre
  803. C commun qu'a deux elements au plus.
  804. C On determine egalement le contour/enveloppe, IPT4, de la cavite
  805. NNBELE2=NBELE2
  806. 2221 CONTINUE
  807. NBELE2=NNBELE2
  808. NNCON=0
  809. JG=IDIM
  810. SEGINI,LREF1
  811. NBNN=IPT3.NUM(/1)-1
  812. NBELEM=NBELE3*3
  813. SEGINI,IPT4
  814. IF (IDIM.EQ.2) IPT4.ITYPEL=2
  815. IF (IDIM.EQ.3) IPT4.ITYPEL=4
  816. NBELE4=0
  817. C On parcoure les elements de la cavite
  818. DO I1=1,IPT3.NUM(/2)
  819. C Nombre d'adjacents a l'element I1
  820. C un element est adjacent s'il partage un cote/face
  821. NVI1=0
  822. C On parcoure les 3/4 cotes/faces de l'element I1
  823. DO I2=1,IPT3.NUM(/1)
  824. C Nombre d'adjacents a la face I2 (1 ou 0)
  825. NVI2=0
  826. IF (IDIM.EQ.2) THEN
  827. IF (I2.EQ.1) THEN
  828. LREF1.LECT(1)=IPT3.NUM(2,I1)
  829. LREF1.LECT(2)=IPT3.NUM(3,I1)
  830. ELSEIF (I2.EQ.2) THEN
  831. LREF1.LECT(1)=IPT3.NUM(1,I1)
  832. LREF1.LECT(2)=IPT3.NUM(3,I1)
  833. ELSEIF (I2.EQ.3) THEN
  834. LREF1.LECT(1)=IPT3.NUM(1,I1)
  835. LREF1.LECT(2)=IPT3.NUM(2,I1)
  836. ENDIF
  837. ELSEIF (IDIM.EQ.3) THEN
  838. IF (I2.EQ.1) THEN
  839. LREF1.LECT(1)=IPT3.NUM(2,I1)
  840. LREF1.LECT(2)=IPT3.NUM(3,I1)
  841. LREF1.LECT(3)=IPT3.NUM(4,I1)
  842. ELSEIF (I2.EQ.2) THEN
  843. LREF1.LECT(1)=IPT3.NUM(1,I1)
  844. LREF1.LECT(2)=IPT3.NUM(3,I1)
  845. LREF1.LECT(3)=IPT3.NUM(4,I1)
  846. ELSEIF (I2.EQ.3) THEN
  847. LREF1.LECT(1)=IPT3.NUM(1,I1)
  848. LREF1.LECT(2)=IPT3.NUM(2,I1)
  849. LREF1.LECT(3)=IPT3.NUM(4,I1)
  850. ELSEIF (I2.EQ.4) THEN
  851. LREF1.LECT(1)=IPT3.NUM(1,I1)
  852. LREF1.LECT(2)=IPT3.NUM(2,I1)
  853. LREF1.LECT(3)=IPT3.NUM(3,I1)
  854. ENDIF
  855. ENDIF
  856. DO I11=1,IPT3.NUM(/2)
  857. NNCI11=0
  858. IF (I11.EQ.I1) GOTO 2222
  859. DO I22=1,IPT3.NUM(/1)
  860. NREFX=IPT3.NUM(I22,I11)
  861. DO K=1,IDIM
  862. IF (NREFX.EQ.LREF1.LECT(K)) THEN
  863. NNCI11=NNCI11+1
  864. ENDIF
  865. ENDDO
  866. ENDDO
  867. IF (NNCI11.EQ.IDIM) THEN
  868. NVI2=NVI2+1
  869. NVI1=NVI1+1
  870. ENDIF
  871. 2222 CONTINUE
  872. ENDDO
  873. C Cas ou I1 a plus de 1 voisin sur la face I2
  874. IF (NVI2.GT.1) THEN
  875. NBNN=IPT2.NUM(/1)
  876. NBELEM=NBELE2
  877. NBSOUS=0
  878. NBREF=0
  879. SEGADJ,IPT2
  880. C Erreur de triangulation : vérifiez qu'il n'y a pas de
  881. C points confondus
  882. CALL ERREUR(846)
  883. GOTO 2501
  884. ENDIF
  885. C Si la face I2 n'est que sur l'element I1, alors on
  886. C l'ajoute au maillage du contour/enceloppe de la cavite
  887. IF (NVI2.EQ.0) THEN
  888. NBELE4=NBELE4+1
  889. IF (IPT4.NUM(/2).LT.NBELE4) THEN
  890. NBNN=IPT4.NUM(/1)
  891. NBELEM=IPT4.NUM(/2)+100
  892. SEGADJ,IPT4
  893. ENDIF
  894. DO K=1,IPT4.NUM(/1)
  895. IF (IDIM.EQ.2) THEN
  896. IF (I2.EQ.1) THEN
  897. IPT4.NUM(1,NBELE4)=IPT3.NUM(2,I1)
  898. IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1)
  899. ELSEIF (I2.EQ.2) THEN
  900. IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1)
  901. IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1)
  902. ELSEIF (I2.EQ.3) THEN
  903. IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1)
  904. IPT4.NUM(2,NBELE4)=IPT3.NUM(2,I1)
  905. ENDIF
  906. ELSEIF (IDIM.EQ.3) THEN
  907. IF (I2.EQ.1) THEN
  908. IPT4.NUM(1,NBELE4)=IPT3.NUM(2,I1)
  909. IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1)
  910. IPT4.NUM(3,NBELE4)=IPT3.NUM(4,I1)
  911. ELSEIF (I2.EQ.2) THEN
  912. IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1)
  913. IPT4.NUM(2,NBELE4)=IPT3.NUM(3,I1)
  914. IPT4.NUM(3,NBELE4)=IPT3.NUM(4,I1)
  915. ELSEIF (I2.EQ.3) THEN
  916. IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1)
  917. IPT4.NUM(2,NBELE4)=IPT3.NUM(2,I1)
  918. IPT4.NUM(3,NBELE4)=IPT3.NUM(4,I1)
  919. ELSEIF (I2.EQ.4) THEN
  920. IPT4.NUM(1,NBELE4)=IPT3.NUM(1,I1)
  921. IPT4.NUM(2,NBELE4)=IPT3.NUM(2,I1)
  922. IPT4.NUM(3,NBELE4)=IPT3.NUM(3,I1)
  923. ENDIF
  924. ENDIF
  925. ENDDO
  926. ENDIF
  927. ENDDO
  928. C Cas ou I1 n'a aucun voisin => cavite non connexe, on retire
  929. C l'element de la cavite
  930. IF (NVI1.EQ.0) THEN
  931. IF (NBELE3.NE.1) THEN
  932. NNCON=1
  933. NESUP=LCAVI.LECT(I1)
  934. ENDIF
  935. ENDIF
  936. ENDDO
  937. NBNN=IPT4.NUM(/1)
  938. NBELEM=NBELE4
  939. SEGADJ,IPT4
  940. IF (NNCON.EQ.1) GOTO 2341
  941. C
  942. C 2.2.3 Verification : la cavite et son contour/enveloppe doivent
  943. C avoir le meme nombre de noeuds
  944. CALL ECROBJ('MAILLAGE',IPT3)
  945. IF (IERR.NE.0) RETURN
  946. CALL NBNO
  947. SEGACT,IPT3
  948. CALL LIRENT(NBNO3,1,IRETOU)
  949. CALL ECROBJ('MAILLAGE',IPT4)
  950. IF (IERR.NE.0) RETURN
  951. CALL NBNO
  952. SEGACT,IPT4
  953. CALL LIRENT(NBNO4,1,IRETOU)
  954. IF (NBNO3.NE.NBNO4) THEN
  955. NBNN=IPT2.NUM(/1)
  956. NBELEM=NBELE2
  957. NBSOUS=0
  958. NBREF=0
  959. SEGADJ,IPT2
  960. C Erreur de triangulation : vérifiez qu'il n'y a pas de
  961. C points confondus
  962. CALL ERREUR(846)
  963. GOTO 2501
  964. ENDIF
  965. C
  966. C------- 2.3 Construction de la boule
  967. C
  968. C 2.3.1 Construction de la liste des elements de la boule LBOUL
  969. C La boule est formee par les elements construits en joignant
  970. C NPI1 aux elements du contour de la cavite.
  971. NBOUL=NBELE4
  972. JG=NBOUL
  973. SEGINI,LBOUL
  974. C Nouveau nombre d'elements de la triangulation
  975. NEW2=NBELE2-NCAVI+NBOUL
  976. C On regarde si le nombre d'elements dans la cavite est inferieur
  977. C au nombre d'elements de son contour/enveloppe car il faudra un
  978. C traitement special lors de la creation de la boule
  979. N43=NBELE4-NBELE3
  980. C Cas general ou l'on va ajouter de nouveaux elements a IPT2
  981. IF (N43.GE.0) THEN
  982. DO J=1,NBOUL
  983. C Les premiers elements de LBOUL sont ceux de LCAVI
  984. IF (J.LE.NCAVI) THEN
  985. LBOUL.LECT(J)=LCAVI.LECT(J)
  986. C et je complete LBOUL avec les nouveaux numeros d'elements
  987. ELSE
  988. LBOUL.LECT(J)=NBELE2+1
  989. NBELE2=NBELE2+1
  990. ENDIF
  991. ENDDO
  992. C J'ajuste IPT2, ITRC1 et ILEA1 si besoin
  993. IF (NEW2.GT.IPT2.NUM(/2)) THEN
  994. NBNN=IPT2.NUM(/1)
  995. NBELEM=NBELE2+ICMIN
  996. SEGADJ,IPT2
  997. NBE1=NBELEM
  998. SEGADJ,ITRC1
  999. NBL1=NBELEM
  1000. SEGADJ,ILEA1
  1001. ENDIF
  1002. C Cas particulier ou l'on va enlever des elements a IPT2
  1003. ELSE
  1004. NSUP=NCAVI-NBOUL
  1005. JG=NSUP
  1006. SEGINI,LSUP
  1007. C LBOUL est fait des NBOUL premiers elements de LCAVI
  1008. DO J=1,NCAVI
  1009. NELJ=LCAVI.LECT(J)
  1010. IF (J.LE.NBOUL) THEN
  1011. LBOUL.LECT(J)=NELJ
  1012. ELSE
  1013. LSUP.LECT(J-NBOUL)=NELJ
  1014. C Mise a zero des numeros des noeuds de l'element NELJ
  1015. DO K=1,IPT2.NUM(/1)
  1016. IPT2.NUM(K,NELJ)=0
  1017. ENDDO
  1018. ENDIF
  1019. ENDDO
  1020. ENDIF
  1021. NBELE2=NEW2
  1022. C
  1023. C 2.3.2 Je recherche les elements exterieurs a la cavite
  1024. C partageant une face de l'enveloppe IPT4, pour cela on regarde
  1025. C les elements de IPT3 et leurs adjacents grace a la table de
  1026. C connectivite
  1027. JG=NBELE4
  1028. SEGINI,LEADJ
  1029. SEGINI,LNVADJ
  1030. JG=IDIM
  1031. SEGINI,LREF1
  1032. C On parcoure les elements II de IPT4 du contour/enveloppe
  1033. DO II=1,NBELE4
  1034. C Numeros des noeuds JJ de l'element II
  1035. DO JJ=1,IDIM
  1036. LREF1.LECT(JJ)=IPT4.NUM(JJ,II)
  1037. ENDDO
  1038. C On parcoure les elements J de IPT3 de la cavite
  1039. DO J=1,NBELE3
  1040. NCOM=0
  1041. NSOM=0
  1042. DO K=1,IDIMP1
  1043. NREFX=IPT3.NUM(K,J)
  1044. C Test si le noeud K de l'element J est commun a II
  1045. DO L=1,IDIM
  1046. IF (NREFX.EQ.LREF1.LECT(L)) THEN
  1047. NCOM=NCOM+1
  1048. NSOM=NSOM+K
  1049. ENDIF
  1050. ENDDO
  1051. C Si l'on a trouve l'element appuye sur la face II
  1052. IF (NCOM.EQ.IDIM) THEN
  1053. NELK2=LCAVI.LECT(J)
  1054. GOTO 2321
  1055. ENDIF
  1056. ENDDO
  1057. ENDDO
  1058. 2321 CONTINUE
  1059. C Determination du numero du noeud voyeur dans IPT2
  1060. IF (IDIM.EQ.2) NNOK2=6-NSOM
  1061. IF (IDIM.EQ.3) NNOK2=10-NSOM
  1062. LEADJ.LECT(II)=ILEA1.LEAC(NELK2,NNOK2,1)
  1063. LNVADJ.LECT(II)=ILEA1.LEAC(NELK2,NNOK2,2)
  1064. ENDDO
  1065. C
  1066. C 2.3.3 Verification : le contour/enveloppe de la cavite doit
  1067. C etre etoile par rapport au point NPI1
  1068. NESUP=0
  1069. IF (IDIM.EQ.3) THEN
  1070. DO J=1,NBELE4
  1071. NA=IPT4.NUM(1,J)
  1072. NB=IPT4.NUM(2,J)
  1073. NC=IPT4.NUM(3,J)
  1074. XA=XCOOR((NA-1)*(IDIM+1)+1)
  1075. YA=XCOOR((NA-1)*(IDIM+1)+2)
  1076. ZA=XCOOR((NA-1)*(IDIM+1)+3)
  1077. XB=XCOOR((NB-1)*(IDIM+1)+1)
  1078. YB=XCOOR((NB-1)*(IDIM+1)+2)
  1079. ZB=XCOOR((NB-1)*(IDIM+1)+3)
  1080. XC=XCOOR((NC-1)*(IDIM+1)+1)
  1081. YC=XCOOR((NC-1)*(IDIM+1)+2)
  1082. ZC=XCOOR((NC-1)*(IDIM+1)+3)
  1083. C Normale au plan ABC : n = AB^AC
  1084. XN=((YB-YA)*(ZC-ZA))-((ZB-ZA)*(YC-YA))
  1085. YN=((ZB-ZA)*(XC-XA))-((XB-XA)*(ZC-ZA))
  1086. ZN=((XB-XA)*(YC-YA))-((YB-YA)*(XC-XA))
  1087. XXN=SQRT((XN**2)+(YN**2)+(ZN**2))
  1088. XN=XN/XXN
  1089. YN=YN/XXN
  1090. ZN=ZN/XXN
  1091. C Test si NPI1 est dans le plan ABC
  1092. PS1=(XN*(XPI1(1)-XA))+(YN*(XPI1(2)-YA))+(ZN*(XPI1(3)-ZA))
  1093. XX0=SQRT(((XPI1(1)-XA)**2)+((XPI1(2)-YA)**2)+
  1094. & ((XPI1(3)-ZA)**2))
  1095. PS11=PS1/XX0
  1096. IF (ABS(PS11).LT.ZERO) THEN
  1097. NESUP=ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),1)
  1098. GOTO 2341
  1099. ENDIF
  1100. C Point D voyant la face ABC mais situe hors de la cavite
  1101. IF (LEADJ.LECT(J).EQ.0) GOTO 2331
  1102. ND=IPT2.NUM(LNVADJ.LECT(J),LEADJ.LECT(J))
  1103. XD=XCOOR((ND-1)*(IDIM+1)+1)
  1104. YD=XCOOR((ND-1)*(IDIM+1)+2)
  1105. ZD=XCOOR((ND-1)*(IDIM+1)+3)
  1106. C Test si NPI1 est du meme cote que D par rapport au plan
  1107. C ABC
  1108. PS2=(XN*(XD-XA)) + (YN*(YD-YA)) + (ZN*(ZD-ZA))
  1109. XX0=SQRT(((XD-XA)**2)+((YD-YA)**2)+((ZD-ZA)**2))
  1110. PS22=PS2/XXN
  1111. IF (ABS(PS22).LT.ZERO) PS2=0.
  1112. IF ((SIGN(1D0,PS1)).EQ.(SIGN(1D0,PS2))) THEN
  1113. NESUP=ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),1)
  1114. GOTO 2341
  1115. ELSE
  1116. GOTO 2331
  1117. ENDIF
  1118. 2331 CONTINUE
  1119. ENDDO
  1120. ENDIF
  1121. C
  1122. C 2.3.4 Suppression d'un element (NESUP) de la cavite
  1123. 2341 CONTINUE
  1124. IF (NESUP.NE.0) THEN
  1125. C Test si NESUP est dans la base, pas touche a la base !
  1126. DO K=1,NBASE
  1127. IF (LBASE.LECT(K).EQ.NESUP) THEN
  1128. NBNN=IPT2.NUM(/1)
  1129. NBELEM=NBELE2
  1130. NBSOUS=0
  1131. NBREF=0
  1132. SEGADJ,IPT2
  1133. C Erreur de triangulation : vérifiez qu'il n'y a pas de
  1134. C points confondus
  1135. CALL ERREUR(846)
  1136. GOTO 2501
  1137. ENDIF
  1138. ENDDO
  1139. NBNN=IPT3.NUM(/1)
  1140. NBELEM=NBELE3-1
  1141. SEGINI,IPT5
  1142. IPT5.ITYPEL=IPT3.ITYPEL
  1143. C Suppression de l'element NESUP le la liste LCAVI
  1144. NS=0
  1145. IF (LCAVI.LECT(NCAVI).EQ.NESUP) THEN
  1146. NCAVI=NCAVI-1
  1147. NS=1
  1148. DO JJ=1,NCAVI
  1149. DO KK=1,IPT3.NUM(/1)
  1150. IPT5.NUM(KK,JJ)=IPT3.NUM(KK,JJ)
  1151. ENDDO
  1152. ENDDO
  1153. ELSE
  1154. DO K=1,NCAVI-1
  1155. IF (LCAVI.LECT(K).EQ.NESUP) THEN
  1156. NS=NS+1
  1157. NCAVI=NCAVI-1
  1158. ENDIF
  1159. LCAVI.LECT(K)=LCAVI.LECT(K+NS)
  1160. DO KK=1,IPT3.NUM(/1)
  1161. IPT5.NUM(KK,K)=IPT3.NUM(KK,K+NS)
  1162. ENDDO
  1163. ENDDO
  1164. ENDIF
  1165. C Ajustement de la liste LCAVI
  1166. JG=NCAVI
  1167. SEGADJ,LCAVI
  1168. SEGSUP,IPT3
  1169. IPT3=IPT5
  1170. NBELE3=NCAVI
  1171. C On repart calculer le contour/enveloppe de la cavite
  1172. GOTO 2221
  1173. ENDIF
  1174. C
  1175. C 2.3.5 Je parcours les elements du contour/enveloppe pour
  1176. C construire le maillage de la boule
  1177. DO J=1,NBELE4
  1178. JEL=LBOUL.LECT(J)
  1179. C Attribution des premiers noeuds des nouveux elements
  1180. C de IPT2 : on choisi ceux du contour IPT4
  1181. DO K=1,IPT4.NUM(/1)
  1182. IPT2.NUM(K,JEL)=IPT4.NUM(K,J)
  1183. ENDDO
  1184. C Attribution du dernier noeud des nouveaux elements de IPT2 :
  1185. C il s'agit de NPI1 au entre de la cavite
  1186. IPT2.NUM(IDIMP1,JEL)=NPI1
  1187. C Mise a jour de la table de connectivite (partie 1 sur 2)
  1188. C uniquement sur les derniers noeuds (IDIMP1) des elements de
  1189. C la boule (adjacents extexrieurs a la boule).
  1190. C Pour cela, on utilise LEADJ et LNVADJ calcules juste avant
  1191. ILEA1.LEAC(JEL,IDIMP1,1)=LEADJ.LECT(J)
  1192. ILEA1.LEAC(JEL,IDIMP1,2)=LNVADJ.LECT(J)
  1193. IF (LEADJ.LECT(J).NE.0) THEN
  1194. ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),1)=JEL
  1195. ILEA1.LEAC(LEADJ.LECT(J),LNVADJ.LECT(J),2)=IDIMP1
  1196. ENDIF
  1197. C et je rajoute au segment ITRC1 ce qu'il faut :
  1198. C
  1199. CALL BCIRCO(IRN1,IPT2,JEL,XCBJ,XRBJ)
  1200. C
  1201. ITRC1.TRC(JEL,1)=XRBJ
  1202. DO K=1,IDIM
  1203. ITRC1.TRC(JEL,K+1)=XCBJ(K)
  1204. ENDDO
  1205. ENDDO
  1206. C
  1207. C 2.3.6 Mise a jour de la table de connectivite (partie 2 sur 2)
  1208. C On parcoure les elements de la boule et on determine toutes les
  1209. C elements adjacents internes a la boule
  1210. JG=IDIM-1
  1211. SEGINI,LREF1
  1212. DO K=1,NBOUL
  1213. NELK=LBOUL.LECT(K)
  1214. C On recherche l'adjacent de NELK vu par le noeud M
  1215. DO M=1,IDIM
  1216. C Numeros des noeuds de l'element NELK a observer
  1217. IF (IDIM.EQ.2) THEN
  1218. IF (M.EQ.1) THEN
  1219. LREF1.LECT(1)=IPT2.NUM(2,NELK)
  1220. ELSEIF (M.EQ.2) THEN
  1221. LREF1.LECT(1)=IPT2.NUM(1,NELK)
  1222. ENDIF
  1223. ENDIF
  1224. IF (IDIM.EQ.3) THEN
  1225. IF (M.EQ.1) THEN
  1226. LREF1.LECT(1)=IPT2.NUM(2,NELK)
  1227. LREF1.LECT(2)=IPT2.NUM(3,NELK)
  1228. ELSEIF (M.EQ.2) THEN
  1229. LREF1.LECT(1)=IPT2.NUM(1,NELK)
  1230. LREF1.LECT(2)=IPT2.NUM(3,NELK)
  1231. ELSEIF (M.EQ.3) THEN
  1232. LREF1.LECT(1)=IPT2.NUM(1,NELK)
  1233. LREF1.LECT(2)=IPT2.NUM(2,NELK)
  1234. ENDIF
  1235. ENDIF
  1236. C On parcoure les autres elements de la boule
  1237. NNOKK=0
  1238. DO KK=1,NBOUL
  1239. NELKK=LBOUL.LECT(KK)
  1240. IF (NELKK.EQ.NELK) GOTO 2361
  1241. NCOM=1
  1242. NSOMKK=IDIMP1
  1243. C On parcoure les 2/3 premiers noeuds LL de NELKK
  1244. DO LL=1,IDIM
  1245. NREFX=IPT2.NUM(LL,NELKK)
  1246. C Test si le noeud LL de NELKK est commun a NELK
  1247. DO L=1,IDIM-1
  1248. IF (NREFX.EQ.LREF1.LECT(L)) THEN
  1249. NCOM=NCOM+1
  1250. NSOMKK=NSOMKK+LL
  1251. ENDIF
  1252. ENDDO
  1253. C Test si NELKK est bien un element adjacent a NELK
  1254. IF (NCOM.EQ.IDIM) GOTO 2362
  1255. ENDDO
  1256. 2361 CONTINUE
  1257. ENDDO
  1258. 2362 CONTINUE
  1259. C Determination du numero du noeud de voyeur de K dans IPT2
  1260. IF (IDIM.EQ.2) NNOKK=6-NSOMKK
  1261. IF (IDIM.EQ.3) NNOKK=10-NSOMKK
  1262. ILEA1.LEAC(NELK,M,1)=NELKK
  1263. ILEA1.LEAC(NELK,M,2)=NNOKK
  1264. ENDDO
  1265. ENDDO
  1266. C
  1267. C 2.3.7 Renumerotation et nettoyage des elements de IPT2 dans le
  1268. C cas particulier ou le nombre d'elements a diminue
  1269. IF (N43.LT.0) THEN
  1270. DO J=1,IPT2.NUM(/2)
  1271. C Determination du nombre d'elements a decaler
  1272. NDEC=0
  1273. DO K=1,NSUP
  1274. C Si l'element examine n'est plus dans la triangulation
  1275. C alors on ne fait rien
  1276. IF (J.EQ.LSUP.LECT(K)) GOTO 2371
  1277. IF (J.GE.LSUP.LECT(K)) NDEC=NDEC+1
  1278. ENDDO
  1279. C Nouveau numero que portera l'element J
  1280. JJ=J-NDEC
  1281. C Renumerotation dans IPT2
  1282. DO K=1,IPT2.NUM(/1)
  1283. IPT2.NUM(K,JJ)=IPT2.NUM(K,J)
  1284. ENDDO
  1285. C Renumerotation dans ILEA1
  1286. DO K=1,IDIMP1
  1287. IOLDVJK=ILEA1.LEAC(J,K,1)
  1288. C Determination du nombre d'elements a decaler
  1289. NDEC2=0
  1290. DO KK=1,NSUP
  1291. IF (IOLDVJK.EQ.LSUP.LECT(KK)) THEN
  1292. NBNN=IPT2.NUM(/1)
  1293. NBELEM=NBELE2
  1294. NBSOUS=0
  1295. NBREF=0
  1296. SEGADJ,IPT2
  1297. C Erreur de triangulation : vérifiez qu'il n'y a
  1298. C pas de points confondus
  1299. CALL ERREUR(846)
  1300. GOTO 2501
  1301. ENDIF
  1302. IF (IOLDVJK.GE.LSUP.LECT(KK)) NDEC2=NDEC2+1
  1303. ENDDO
  1304. NEWVJK=IOLDVJK-NDEC2
  1305. ILEA1.LEAC(JJ,K,1)=NEWVJK
  1306. ILEA1.LEAC(JJ,K,2)=ILEA1.LEAC(J,K,2)
  1307. ENDDO
  1308. C Renumerotation dans ITCR1
  1309. ITRC1.TRC(JJ,1)=ITRC1.TRC(J,1)
  1310. DO K=1,IDIM
  1311. ITRC1.TRC(JJ,K+1)=ITRC1.TRC(J,K+1)
  1312. ENDDO
  1313. 2371 CONTINUE
  1314. ENDDO
  1315. NBNN=IPT2.NUM(/1)
  1316. NBELEM=NBELE2
  1317. SEGADJ,IPT2
  1318. NBL1=NBELE2
  1319. SEGADJ,ILEA1
  1320. NBE1=NBELE2
  1321. SEGADJ,ITRC1
  1322. ENDIF
  1323. SEGSUP,IPT3,IPT4
  1324. C
  1325. C 2.3.8 Cas particulier de la dimension 1, on ajoute le point
  1326. C NPI1 au maillage
  1327. 2381 IF (IDIM.EQ.1) THEN
  1328. C Nouveau nombre d'element dans IPT2
  1329. NEW2=NBELE2+1
  1330. C J'ajuste IPT2, ITRC1 et ILEA1 si besoin
  1331. IF (NEW2.GT.IPT2.NUM(/2)) THEN
  1332. NBNN=IPT2.NUM(/1)
  1333. NBELEM=NBELE2+ICMIN
  1334. SEGADJ,IPT2
  1335. NBE1=NBELEM
  1336. SEGADJ,ITRC1
  1337. NBL1=NBELEM
  1338. SEGADJ,ILEA1
  1339. ENDIF
  1340. NBELE2=NEW2
  1341. C Modification de IPT2 pour ajouter un element
  1342. IPT2.NUM(1,NBASE1)=NREFA
  1343. IPT2.NUM(2,NBASE1)=NPI1
  1344. IPT2.NUM(1,NBELE2)=NPI1
  1345. IPT2.NUM(2,NBELE2)=NREFB
  1346. C Mise a jour de la table de connectivite ILEA1
  1347. NED=ILEA1.LEAC(NBASE1,1,1)
  1348. ILEA1.LEAC(NBASE1,1,1)=NBELE2
  1349. ILEA1.LEAC(NBASE1,1,2)=2
  1350. ILEA1.LEAC(NBELE2,2,1)=NBASE1
  1351. ILEA1.LEAC(NBELE2,2,2)=1
  1352. IF (NED.EQ.0) THEN
  1353. ILEA1.LEAC(NBELE2,1,1)=0
  1354. ILEA1.LEAC(NBELE2,1,2)=0
  1355. ELSE
  1356. ILEA1.LEAC(NBELE2,1,1)=NED
  1357. ILEA1.LEAC(NBELE2,1,2)=2
  1358. ILEA1.LEAC(NED,2,1)=NBELE2
  1359. ILEA1.LEAC(NED,2,2)=1
  1360. ENDIF
  1361. ENDIF
  1362. ENDDO
  1363. C Ajustement final de IPT2, ITRC1 et ILEA1
  1364. IF (IDIM.NE.1) THEN
  1365. NBNN=IPT2.NUM(/1)
  1366. NBELEM=NBELE2
  1367. SEGADJ,IPT2
  1368. NBE1=NBELE2
  1369. SEGADJ,ITRC1
  1370. NBL1=NBELE2
  1371. SEGADJ,ILEA1
  1372. ENDIF
  1373. C
  1374. C
  1375. C------- 2.4 Retour du maillage dans le cas ou l'on garde la boite
  1376. C Si IBOI egal a 1, on renvoie IPT2 (maillage avec boite)
  1377. 2501 CONTINUE
  1378. IF (IERR.NE.0) GOTO 999
  1379. IF (IBOI.EQ.1) THEN
  1380. IF (IDIM.NE.1) THEN
  1381. SEGSUP,LBASE,LCAVI,LBOUL,LREF1,LEADJ,LNVADJ
  1382. SEGSUP,LSUP,LSCP
  1383. SEGDES,IPT2
  1384. ENDIF
  1385. GOTO 999
  1386. ENDIF
  1387. C
  1388. C
  1389. C
  1390. C----------------------------------------------------------------------C
  1391. C PARTIE 3 C
  1392. C Elimination des elements appuyes sur les coins de la boite C
  1393. C----------------------------------------------------------------------C
  1394. C
  1395. C------- 3.1 Suppression des elements touchant les noeuds de la boite
  1396. C Je parcours les noeuds de la boite en cherchant les elements du
  1397. C maillage IPT2 qui en contiennent. Ceux qui n'en contiennent pas
  1398. C sont recopies dans IPT3, qui sera le maillage de sortie :
  1399. NBSOUS=0
  1400. NBREF=0
  1401. NBNN=IPT2.NUM(/1)
  1402. NBELEM=NBELE2
  1403. SEGINI,IPT3
  1404. IPT3.ITYPEL=IPT2.ITYPEL
  1405. NBELE3=0
  1406. DO I=1,NBELE2
  1407. DO J=1,IPT2.NUM(/1)
  1408. NPI2=IPT2.NUM(J,I)
  1409. DO K=1,NNBOIT
  1410. IF (NPI2.EQ.LNBOIT(K)) GOTO 3101
  1411. ENDDO
  1412. ENDDO
  1413. NBELE3=NBELE3+1
  1414. DO J=1,IPT2.NUM(/1)
  1415. IPT3.NUM(J,NBELE3)=IPT2.NUM(J,I)
  1416. ENDDO
  1417. 3101 CONTINUE
  1418. ENDDO
  1419. NBNN=IPT3.NUM(/1)
  1420. NBELEM=NBELE3
  1421. SEGADJ,IPT3
  1422. C
  1423. C------- 3.2 Cas ou IPT3 ne contient pas d'elements.
  1424. IF (IPT3.NUM(/2).EQ.0) THEN
  1425. C En dimension 2, cela veut dire que les points de IPT1 sont
  1426. C alignes. On va alors rendre le maillage des lignes ne touchant
  1427. C pas les noeuds de la boite
  1428. IF (IDIM.EQ.2) GOTO 3203
  1429. C En dimension 3, cela veut dire que les points sont coplanaires.
  1430. C On refait la meme operation pour rendre le maillage faces ne
  1431. C touchant pas les noeuds de la boite
  1432. JG=IDIM
  1433. SEGINI,LREF1
  1434. JG=IPT2.NUM(/2)
  1435. SEGINI,LSUP
  1436. NBSOUS=0
  1437. NBREF=0
  1438. NBNN=IDIM
  1439. NBELEM=NBELE2
  1440. SEGINI,IPT3
  1441. C On va rendre un maillage de TRI3
  1442. IPT3.ITYPEL=4
  1443. NBELE3=0
  1444. C Boucle sur les elements de la triangulation avec la boite
  1445. DO I=1,NBELE2
  1446. NNS=0
  1447. NSOM=0
  1448. NEL=0
  1449. C On regarde si le noeud J de l'element I fait partie de ceux
  1450. C de la boite
  1451. DO J=1,IPT2.NUM(/1)
  1452. NPI2=IPT2.NUM(J,I)
  1453. DO K=1,NNBOIT
  1454. IF (NPI2.EQ.LNBOIT(K)) THEN
  1455. NNS=NNS+1
  1456. NSOM=J
  1457. NEL=ILEA1.LEAC(I,J,1)
  1458. ENDIF
  1459. ENDDO
  1460. ENDDO
  1461. C Si l'element I n'a qu'un noeud faisant partie de ceux de la
  1462. C boite, alors on cree l'element appuyes sur ses autres noeuds
  1463. IF (NNS.EQ.1) THEN
  1464. DO K=1,LSUP.LECT(/1)
  1465. IF (LSUP.LECT(K).EQ.0) GOTO 3201
  1466. IF (I.EQ.LSUP.LECT(K)) THEN
  1467. GOTO 3202
  1468. ENDIF
  1469. ENDDO
  1470. 3201 CONTINUE
  1471. NBELE3=NBELE3+1
  1472. LSUP.LECT(NBELE3)=NEL
  1473. IF (NSOM.EQ.1) THEN
  1474. LREF1.LECT(1)=2
  1475. LREF1.LECT(2)=3
  1476. LREF1.LECT(3)=4
  1477. ELSEIF (NSOM.EQ.2) THEN
  1478. LREF1.LECT(1)=1
  1479. LREF1.LECT(2)=3
  1480. LREF1.LECT(3)=4
  1481. ELSEIF (NSOM.EQ.3) THEN
  1482. LREF1.LECT(1)=1
  1483. LREF1.LECT(2)=2
  1484. LREF1.LECT(3)=4
  1485. ELSEIF (NSOM.EQ.4) THEN
  1486. LREF1.LECT(1)=1
  1487. LREF1.LECT(2)=2
  1488. LREF1.LECT(3)=3
  1489. ENDIF
  1490. DO J=1,IPT3.NUM(/1)
  1491. IPT3.NUM(J,NBELE3)=IPT2.NUM(LREF1.LECT(J),I)
  1492. ENDDO
  1493. ENDIF
  1494. 3202 CONTINUE
  1495. ENDDO
  1496. NBNN=IPT3.NUM(/1)
  1497. NBELEM=NBELE3
  1498. SEGADJ,IPT3
  1499. ENDIF
  1500. C
  1501. C Cas ou les points de IPT1 sont alignes. On change le maillage en
  1502. C lignes et ne garde que les lignes ne touchant pas les noeuds de la
  1503. C boite
  1504. 3203 CONTINUE
  1505. IF (IPT3.NUM(/2).EQ.0) THEN
  1506. CALL ECROBJ('MAILLAGE',IPT2)
  1507. CALL CHANLG
  1508. CALL LIROBJ('MAILLAGE',IPT4,1,IRETOU)
  1509. IF (IERR.NE.0) RETURN
  1510. SEGACT,IPT4
  1511. NBELE4=IPT4.NUM(/2)
  1512. NBSOUS=0
  1513. NBREF=0
  1514. NBNN=IPT4.NUM(/1)
  1515. NBELEM=NBELE4
  1516. SEGINI,IPT3
  1517. C On va rendre un maillage de SEG2
  1518. IPT3.ITYPEL=IPT4.ITYPEL
  1519. NBELE3=0
  1520. C Boucle sur les lignes de la triangulation avec la boite
  1521. DO I=1,NBELE4
  1522. NNS=0
  1523. C On regarde si les noeuds l'element I font partie de ceux de
  1524. C la boite
  1525. DO J=1,IPT4.NUM(/1)
  1526. NPI4=IPT4.NUM(J,I)
  1527. DO K=1,NNBOIT
  1528. IF (NPI4.EQ.LNBOIT(K)) THEN
  1529. NNS=NNS+1
  1530. ENDIF
  1531. ENDDO
  1532. ENDDO
  1533. C Si l'element I n'a aucun noeud faisant partie de ceux de la
  1534. C boite, alors on cree l'element appuyes sur ses noeuds
  1535. IF (NNS.EQ.0) THEN
  1536. NBELE3=NBELE3+1
  1537. IPT3.NUM(1,NBELE3)=IPT4.NUM(1,I)
  1538. IPT3.NUM(2,NBELE3)=IPT4.NUM(2,I)
  1539. ENDIF
  1540. ENDDO
  1541. NBNN=IPT3.NUM(/1)
  1542. NBELEM=NBELE3
  1543. SEGADJ,IPT3
  1544. SEGDES,IPT4
  1545. ENDIF
  1546. C
  1547. C Un peu de menage avant de quitter
  1548. IF (IDIM.EQ.1) THEN
  1549. SEGSUP,LSUP
  1550. SEGDES,IPT2
  1551. ELSE
  1552. SEGDES,IPT3
  1553. SEGSUP,IPT2,LBASE,LCAVI,LBOUL,LREF1,LEADJ,LNVADJ
  1554. SEGSUP,LSUP,LSCP
  1555. ENDIF
  1556. IPT2=IPT3
  1557. NBPTS=NBPTS0
  1558. SEGADJ,MCOORD
  1559. C
  1560. 999 RETURN
  1561. END
  1562.  
  1563.  
  1564.  
  1565.  
  1566.  
  1567.  
  1568.  
  1569.  
  1570.  

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