Télécharger intgeo.eso

Retour à la liste

Numérotation des lignes :

  1. C INTGEO SOURCE BP208322 16/11/18 21:17:51 9177
  2. C
  3. SUBROUTINE INTGEO(IOB1,IOB2,IOB3,IOB4,ICLE)
  4. C-----------------------------------------------------------------------
  5. C
  6. C INTERSECTION GEOMETRIQUE DE MAILLAGE
  7. C
  8. C-----------------------------------------------------------------------
  9. C
  10. C ENTREE :
  11. C IOB1 MAILLAGE 1 INITIAL
  12. C (IOB2 MAILLAGE 2 INITIAL)
  13. C
  14. C SORTIE :
  15. C IOB1 MAILLAGE 1 MODIFIE
  16. C (IOB2 MAILLAGE 2 MODIFIE)
  17. C IOB3 MAILLAGE 3 = SELF-INTERSECTION DU MAILLAGE 1
  18. C ou INTERSECTION DU MAILLAGE 1 AVEC 2
  19. C (IOB4 MAILLAGE 3 = INTERSECTION DU MAILLAGE 2 AVEC 1)
  20. C
  21. C-----------------------------------------------------------------------
  22. C
  23. C BP208322, 28/08/2012 : Creation
  24. C
  25. C REM : - self-intersection abandonnée car compliquée : on préfère
  26. C chercher à créer directement une bonne surface ...
  27. C
  28. C TO DO : - cas coplanaire à traiter
  29. c - optimiser
  30. C
  31. C-----------------------------------------------------------------------
  32. c
  33. IMPLICIT INTEGER(I-N)
  34. IMPLICIT REAL*8 (A-H,O-Z)
  35. -INC SMCOORD
  36. -INC CCOPTIO
  37. -INC CCGEOME
  38. -INC SMELEME
  39. -INC CCREEL
  40.  
  41. SEGMENT WORKEL
  42. INTEGER IDEJVU(NBELEM)
  43. REAL*8 V(3,NBELEM)
  44. ENDSEGMENT
  45. POINTEUR WORK1.WORKEL,WORK2.WORKEL
  46. INTEGER NODE(2)
  47.  
  48. c NBNMAX = NomBre de Noeuds MAX
  49. c ITA(i) = numero du i^eme noeuds du triangle A courant
  50. c XTA(i,k) = k^ieme coordonnee du i^eme noeuds du triangle A courant
  51. c XTB(i,k) = k^ieme coordonnee du i^eme noeuds du triangle B courant
  52. PARAMETER (NBNMAX=3)
  53. INTEGER ITA(NBNMAX),ITB(NBNMAX)
  54. REAL*8 XTA(NBNMAX,3),XTB(NBNMAX,3)
  55. REAL*8 XSEG(1*2,3)
  56. REAL*8 dA1(3),dA2(3),dA3(3)
  57. REAL*8 dB1(3),dB2(3),dB3(3)
  58. REAL*8 xP2(3),xP3(3),xQ2(3),xQ3(3)
  59. LOGICAL SELF
  60.  
  61. if(iimpi.GE.10) write(ioimp,*)'====== ENTREE dans INTGEO ======'
  62.  
  63. C-----------------------------------------------------------------------
  64. C RECUP + VERIFICATION DES DONNEES D ENTREE
  65.  
  66. c SELF intersection ?
  67. SELF = (IOB2.eq.0)
  68. if(SELF) then
  69. write(ioimp,*) 'Self-intersection non traitée'
  70. MOTERR(1:8)='MAILLAGE'
  71. call erreur(37)
  72. endif
  73. c on ne traite pas aujourd'hui... abandonnée car compliquée !
  74.  
  75. c option NOVE => IVERI
  76. IVERI = ICLE
  77.  
  78. c Critere absolu de proximite entre 2 noeuds
  79. CALL LIRREE(XELIM,0,IELIM)
  80. c Critere absolu de proximite d'un plan (detection du cas coplanaire)
  81. CALL LIRREE(XPLAN,0,IPLAN)
  82. if(IPLAN.eq.0) XPLAN=XZPREC**(0.66)
  83. if(IELIM.ne.0) then
  84. XELIM2=2.d0*XELIM
  85. XPLAN2=0.01d0*XELIM
  86. if(XPLAN2.lt.XPLAN) then
  87. XPLAN=XPLAN2
  88. endif
  89. endif
  90. XPLANn=-1.d0*XPLAN
  91. c si XELIM non fourni on calcule avec la taille du 1er element tq :
  92. c XZPREC << XPLAN << XELIM << taille des elements
  93.  
  94. C Operateur disponible en dimension 3
  95. IF (IDIM.NE.3) THEN
  96. INTERR(1)=IDIM
  97. CALL ERREUR(709)
  98. RETURN
  99. ENDIF
  100. idim1 = idim+1
  101. segact,MCOORD*mod
  102. NBPTS=XCOOR(/1)/idim1
  103.  
  104. C Recup
  105. IPT1 = IOB1
  106. IPT2 = IOB2
  107. IPT3 = 0
  108. IPT4 = 0
  109.  
  110. C +VERIFICATION DE IPT1
  111. C un maillage non complexe
  112. segact,IPT1
  113. NBSOU1=IPT1.LISOUS(/1)
  114. if (NBSOU1.gt.0) then
  115. call erreur(25)
  116. return
  117. endif
  118. C maillage de TRI3 uniquement
  119. ITYPE1=IPT1.ITYPEL
  120. if(ITYPE1.ne.4) then
  121. call erreur(16)
  122. return
  123. endif
  124. NBNN1 = IPT1.NUM(/1)
  125. NBELE1 = IPT1.NUM(/2)
  126.  
  127. C +VERIFICATION DE IPT2
  128. IF(.not.SELF) THEN
  129. C un maillage non complexe
  130. segact,IPT2
  131. NBSOU2=IPT2.LISOUS(/1)
  132. if (NBSOU2.gt.0) then
  133. call erreur(25)
  134. return
  135. endif
  136. C maillage de TRI3 uniquement
  137. ITYPE2=IPT2.ITYPEL
  138. if(ITYPE2.ne.4) then
  139. call erreur(16)
  140. return
  141. endif
  142. NBNN2 = IPT2.NUM(/1)
  143. NBELE2 = IPT2.NUM(/2)
  144. ENDIF
  145.  
  146.  
  147.  
  148. C-----------------------------------------------------------------------
  149. C PREPARATION DES DONNEES DE SORTIE
  150.  
  151. NBSOUS = 0
  152. NBREF = 0
  153.  
  154. C +NOUVEAU IPT1 (dimensionné + grand)
  155. NBNN = NBNN1
  156. NBELEM = NBELE1
  157. segini,MELEME=IPT1
  158. segdes,IPT1
  159. IPT1=MELEME
  160. NBELEM = 4*NBELE1
  161. segadj,IPT1
  162. segini,WORK1
  163.  
  164. C +NOUVEAU IPT2
  165. IF(.not.SELF) THEN
  166. NBNN = NBNN2
  167. NBELEM = NBELE2
  168. segini,MELEME=IPT2
  169. segdes,IPT2
  170. IPT2=MELEME
  171. NBELEM = 4*NBELE2
  172. segadj,IPT2
  173. segini,WORK2
  174. ENDIF
  175.  
  176. C +IPT3 et IPT4 = lignes de seg2 = intersection des tri3
  177. NBNN3 = 2
  178. NBNN = NBNN3
  179. NBELE3 = 0
  180. NBELE4 = 0
  181. NBELEM = 3*NBELE1
  182. segini,IPT3
  183. IPT3.ITYPEL=2
  184. IF(.not.SELF) THEN
  185. segini,IPT4
  186. IPT4.ITYPEL=2
  187. ENDIF
  188.  
  189.  
  190. C-----------------------------------------------------------------------
  191. C POUR LE CAS DE LA SELF-INTERSECTION
  192. IF(SELF) THEN
  193. IPT2=IPT1
  194. NBNN2=NBNN1
  195. NBELE2=NBELE1
  196. WORK2=WORK1
  197. IPT4=IPT3
  198. c attention : on va remplir IPT2 qui ne refere pas au meme element
  199. c que IPT1 mais pas IPT4 (car on remplit deja IPT3)
  200. ENDIF
  201.  
  202. c pour ajuster la taille par bloc
  203. nbloc = max(4,max(NBELE1,NBELE2))
  204.  
  205. C=======================================================================
  206. C BOUCLES POUR LE CALCUL D INTERSECTION
  207.  
  208.  
  209. C==== BOUCLE SUR LES ELEMENTS A ========================================
  210. c on met une grande valeur mais on s'arrete des qu on a fini :
  211. c plus de triangles (initiaux + coupés) a traiter (=NBELE1 atteint)
  212. c IAFIN = 1000*NBELE1
  213. IAFIN = 100*NBELE1+100*NBELE2
  214.  
  215. DO 1000 IA=1,IAFIN
  216.  
  217. C-----------------------------------------------------------------------
  218. C CALCULS RELATIFS AU TRIANGLE A
  219.  
  220. c 1er coup : on prend le triangle "pere" (=avant decoupe)
  221. IBDEB0=0
  222. 1100 CONTINUE
  223. ICOUPA = 0
  224. ICOUPB = 0
  225.  
  226.  
  227. c coordonnees du Triangle A
  228. xAmax=XPETIT
  229. yAmax=XPETIT
  230. zAmax=XPETIT
  231. xAmin=XGRAND
  232. yAmin=XGRAND
  233. zAmin=XGRAND
  234. DO INOD1=1,NBNN1
  235. IPA1 = IPT1.NUM(INOD1,IA)
  236. if(INOD1.gt.1) then
  237. do iprec=1,INOD1-1
  238. if(IPA1.eq.ITA(iprec)) then
  239. c write(ioimp,*) '2 noeuds identiques => on saute'
  240. goto 900
  241. endif
  242. enddo
  243. endif
  244. ITA(INOD1)=IPA1
  245. IREF1= idim1*(IPA1-1)
  246. XTA(INOD1,1)=XCOOR(IREF1+1)
  247. XTA(INOD1,2)=XCOOR(IREF1+2)
  248. XTA(INOD1,3)=XCOOR(IREF1+3)
  249. xAmax = max(xAmax,XTA(INOD1,1))
  250. yAmax = max(yAmax,XTA(INOD1,2))
  251. zAmax = max(zAmax,XTA(INOD1,3))
  252. xAmin = min(xAmin,XTA(INOD1,1))
  253. yAmin = min(yAmin,XTA(INOD1,2))
  254. zAmin = min(zAmin,XTA(INOD1,3))
  255. ENDDO
  256.  
  257. if(iimpi.ge.10) then
  258. write(ioimp,*) '____________________________________________'
  259. write(ioimp,*) 'TRIANGLE A # ',IA,' / ',NBELE1
  260. write(ioimp,FMT="(' noeuds ',I5,I5,I5)")ITA(1),ITA(2),ITA(3)
  261. c write(ioimp,*) ' coordx=',(XTA(iou,1),iou=1,3)
  262. c write(ioimp,*) ' coordy=',(XTA(iou,2),iou=1,3)
  263. c write(ioimp,*) ' coordz=',(XTA(iou,3),iou=1,3)
  264. endif
  265.  
  266. IF(WORK1.IDEJVU(IA).eq.0) THEN
  267.  
  268. c vecteurs des aretes
  269. va12x = XTA(2,1) - XTA(1,1)
  270. va12y = XTA(2,2) - XTA(1,2)
  271. va12z = XTA(2,3) - XTA(1,3)
  272. va13x = XTA(3,1) - XTA(1,1)
  273. va13y = XTA(3,2) - XTA(1,2)
  274. va13z = XTA(3,3) - XTA(1,3)
  275.  
  276. c vecteur normal au Triangle A + aire du triangle
  277. vax = va12y*va13z - va12z*va13y
  278. vay = va12z*va13x - va12x*va13z
  279. vaz = va12x*va13y - va12y*va13x
  280. xlenA = sqrt(vax*vax + vay*vay + vaz*vaz)
  281. if (xlenA.le.XPETIT) then
  282. write(ioimp,*) 'TRI3 A de dimension nulle! ',xlenA,XPETIT
  283. c write(ioimp,*) 'va12=',va12x,va12y,va12z
  284. c write(ioimp,*) 'va13=',va13x,va13y,va13z
  285. c write(ioimp,*) 'va12 PVEC va13 = va=',vax,vay,vaz
  286. if(IVERI.ge.1) then
  287. c write(ioimp,*) 'NOVERIF: ON retourne le maillage en l etat'
  288. c goto 9000
  289. c bp : une autre solution est de "sauter" cet element
  290. write(ioimp,*) 'NOVERIF: l execution continue malgre tout'
  291. goto 900
  292. else
  293. call erreur(26)
  294. return
  295. endif
  296. endif
  297. vax = vax / xlenA
  298. vay = vay / xlenA
  299. vaz = vaz / xlenA
  300. WORK1.V(1,IA) = vax
  301. WORK1.V(2,IA) = vay
  302. WORK1.V(3,IA) = vaz
  303. WORK1.IDEJVU(IA)=1
  304. if(iimpi.ge.12) write(ioimp,*) ' normale vA =',vax,vay,vaz
  305. if(IELIM.eq.0) then
  306. XELIM=1.D-2*xlenA
  307. XELIM2=2.d0*XELIM
  308. IELIM=1
  309. XPLAN2=1.D-2*XELIM
  310. if(XPLAN2.lt.XPLAN) then
  311. XPLAN=XPLAN2
  312. XPLANn=-1.d0*XPLAN
  313. endif
  314. endif
  315.  
  316. ELSE
  317. vax = WORK1.V(1,IA)
  318. vay = WORK1.V(2,IA)
  319. vaz = WORK1.V(3,IA)
  320.  
  321. ENDIF
  322.  
  323. if(IELIM.eq.0) write(ioimp,*)'Valeur par défaut : XELIM=',XELIM
  324. if(IPLAN.eq.0) write(ioimp,*)'Valeur par défaut : XPLAN=',XPLAN
  325. IPLAN=1
  326. if(iimpi.ge.10) write(ioimp,*)'XELIM,XPLAN=',XELIM,XPLAN
  327. if(iimpi.ge.12) write(ioimp,*)'XELIM2,XPLANn=',XELIM2,XPLANn
  328.  
  329. c on definit une boite limitant le domaine de recherche
  330. xAmax = xAmax + (XELIM2)
  331. yAmax = yAmax + (XELIM2)
  332. zAmax = zAmax + (XELIM2)
  333. xAmin = xAmin - (XELIM2)
  334. yAmin = yAmin - (XELIM2)
  335. zAmin = zAmin - (XELIM2)
  336.  
  337.  
  338. C====== BOUCLE SUR LES ELEMENTS B ======================================
  339.  
  340. c dans le cas d'une decoupe, on reprend avec l element IB suivant
  341. IF(IBDEB0.ne.0) THEN
  342. IBDEB=IBDEB0
  343. ELSE
  344. c dans le cas de la self-intersection on teste les triangles d'apres
  345. IF(SELF) THEN
  346. IBDEB=IA+1
  347. ELSE
  348. IBDEB=1
  349. ENDIF
  350. ENDIF
  351. IBFIN=NBELE2
  352.  
  353. DO 2000 IB=IBDEB,IBFIN
  354.  
  355. c if(iimpi.ge.10) write(ioimp,*) '+ TRIANGLE B # ',IB,' / ',NBELE2
  356.  
  357. C-----------------------------------------------------------------------
  358. C CALCUL DE LA DISTANCE AU PLAN TA DES NOEUDS DU Triangle B
  359.  
  360. c coordonnees du triangle B
  361. ibxmax=0
  362. ibymax=0
  363. ibzmax=0
  364. ibxmin=0
  365. ibymin=0
  366. ibzmin=0
  367. DO INOD2=1,NBNN2
  368. IPB2 = IPT2.NUM(INOD2,IB)
  369. if(INOD2.gt.1) then
  370. do iprec=1,INOD2-1
  371. if(IPB2.eq.ITB(iprec)) then
  372. write(ioimp,*) '2 noeuds identiques de TB => on saute'
  373. goto 1900
  374. endif
  375. enddo
  376. endif
  377. ITB(INOD2)=IPB2
  378. IREF2= idim1*(IPB2-1)
  379. XTB(INOD2,1)=XCOOR(IREF2+1)
  380. XTB(INOD2,2)=XCOOR(IREF2+2)
  381. XTB(INOD2,3)=XCOOR(IREF2+3)
  382. if(XTB(INOD2,1).gt.xAmax) ibxmax=ibxmax+1
  383. if(XTB(INOD2,2).gt.yAmax) ibymax=ibymax+1
  384. if(XTB(INOD2,3).gt.zAmax) ibzmax=ibzmax+1
  385. if(XTB(INOD2,1).lt.xAmin) ibxmin=ibxmin+1
  386. if(XTB(INOD2,2).lt.yAmin) ibymin=ibymin+1
  387. if(XTB(INOD2,3).lt.zAmin) ibzmin=ibzmin+1
  388. ENDDO
  389.  
  390. c si on est hors boite selon x y ou z on saute
  391. if(ibxmax.eq.NBNN2.or.ibymax.eq.NBNN2.or.ibzmax.eq.NBNN2)
  392. & goto 1900
  393. if(ibxmin.eq.NBNN2.or.ibymin.eq.NBNN2.or.ibzmin.eq.NBNN2)
  394. & goto 1900
  395.  
  396. c distance au plan TA
  397. do iidim=1,idim
  398. dB1(iidim) = XTB(1,iidim) - XTA(1,iidim)
  399. dB2(iidim) = XTB(2,iidim) - XTA(1,iidim)
  400. dB3(iidim) = XTB(3,iidim) - XTA(1,iidim)
  401. enddo
  402. xddB1 = vax*dB1(1) + vay*dB1(2) + vaz*dB1(3)
  403. xddB2 = vax*dB2(1) + vay*dB2(2) + vaz*dB2(3)
  404. xddB3 = vax*dB3(1) + vay*dB3(2) + vaz*dB3(3)
  405.  
  406. c si tous les points de TB sont en dessous (ou dessus) de TA
  407. c + le cas affleurant
  408. c => il n'y a pas d'intersection : on saute
  409. if(xddB1.lt.XPLANn.and.xddB2.lt.XPLANn.and.xddB3.lt.XPLANn)
  410. & goto 1900
  411. if(xddB1.gt.XPLAN.and.xddB2.gt.XPLAN.and.xddB3.gt.XPLAN)
  412. & goto 1900
  413.  
  414. c combien de points de TB dans le plan de TA ? => iBinA
  415. iBinA = 0
  416. if(abs(xddB1).le.XPLAN) iBinA=iBinA+1
  417. if(abs(xddB2).le.XPLAN) iBinA=iBinA+2
  418. if(abs(xddB3).le.XPLAN) iBinA=iBinA+4
  419.  
  420. c si tous les points de TB sont dans le plan de TA,
  421. c => on doit traiter le cas coplanaire
  422. c Il n'est pas traité pour l'instant car assez rare en pratique en 3D.
  423. if (iBinA.eq.7) then
  424. c L'intersection générée pourrait etre consitutée de polygones [Gander]
  425. c etape 1: calcul des intersection segment de A - segments de B
  426. c 2: recherche des points de A (/B) a l'interieur de TB (/TA)
  427. c 3: ordonner dans le sens trigo les points des étapes 1 et 2
  428. c 4: eliminer les doublons
  429. goto 1900
  430. endif
  431.  
  432. c si 1 seul point de TB est dans le plan de TA,
  433. c et les autres dont du meme coté,
  434. c => on ne traite pas cette intersection ponctuelle : on saute
  435. proB12 = xddB1*xddB2
  436. proB23 = xddB2*xddB3
  437. proB31 = xddB3*xddB1
  438. if(iBinA.eq.1.and.proB23.gt.0.d0) goto 1900
  439. if(iBinA.eq.2.and.proB31.gt.0.d0) goto 1900
  440. if(iBinA.eq.4.and.proB12.gt.0.d0) goto 1900
  441.  
  442. c si 2 points de TB sont dans le plan de TA
  443. c ou si on a une traversée franche => on traite
  444. if(iimpi.ge.10) then
  445. write(ioimp,*) '+ TRIANGLE B # ',IB,' / ',NBELE2
  446. write(ioimp,FMT="('noeuds ',I5,I5,I5)")ITB(1),ITB(2),ITB(3)
  447. c write(ioimp,*) '+ coordx ',(XTB(iou,1),iou=1,3)
  448. c write(ioimp,*) '+ coordy ',(XTB(iou,2),iou=1,3)
  449. c write(ioimp,*) '+ coordz ',(XTB(iou,3),iou=1,3)
  450. if(iimpi.ge.12)
  451. & write(ioimp,*) ' xddB1,2,3=',xddB1,xddB2,xddB3
  452. endif
  453.  
  454. C-----------------------------------------------------------------------
  455. C CALCULS RELATIFS AU TRIANGLE B
  456.  
  457. IF(WORK2.IDEJVU(IB).eq.0) THEN
  458.  
  459. c vecteurs des aretes
  460. vb12x = XTB(2,1) - XTB(1,1)
  461. vb12y = XTB(2,2) - XTB(1,2)
  462. vb12z = XTB(2,3) - XTB(1,3)
  463. vb13x = XTB(3,1) - XTB(1,1)
  464. vb13y = XTB(3,2) - XTB(1,2)
  465. vb13z = XTB(3,3) - XTB(1,3)
  466.  
  467. c vecteur normal au Triangle B + aire du triangle
  468. vbx = vb12y*vb13z - vb12z*vb13y
  469. vby = vb12z*vb13x - vb12x*vb13z
  470. vbz = vb12x*vb13y - vb12y*vb13x
  471. xlenB = sqrt(vbx*vbx + vby*vby + vbz*vbz)
  472. if (xlenB.le.XPETIT) then
  473. write(ioimp,*) 'TRI3 B de dimension nulle! ',xlenB,XPETIT
  474. if(IVERI.ge.1) then
  475. c write(ioimp,*) 'NOVERIF: ON retourne le maillage en l etat'
  476. c goto 9000
  477. c bp : une autre solution est de "sauter" cet element
  478. write(ioimp,*) 'NOVERIF: l execution continue malgre tout'
  479. goto 1900
  480. else
  481. call erreur(26)
  482. return
  483. endif
  484. endif
  485. vbx = vbx / xlenB
  486. vby = vby / xlenB
  487. vbz = vbz / xlenB
  488. WORK2.V(1,IB) = vbx
  489. WORK2.V(2,IB) = vby
  490. WORK2.V(3,IB) = vbz
  491. WORK2.IDEJVU(IB) = 1
  492. if(iimpi.ge.12) write(ioimp,*) ' normale vB =',vbx,vby,vbz
  493.  
  494. ELSE
  495. vbx = WORK2.V(1,IB)
  496. vby = WORK2.V(2,IB)
  497. vbz = WORK2.V(3,IB)
  498.  
  499. ENDIF
  500.  
  501. C-----------------------------------------------------------------------
  502. C CALCUL DE LA DISTANCE AU PLAN TB DES NOEUDS DU Triangle A
  503.  
  504. do iidim=1,idim
  505. dA1(iidim) = XTA(1,iidim) - XTB(1,iidim)
  506. dA2(iidim) = XTA(2,iidim) - XTB(1,iidim)
  507. dA3(iidim) = XTA(3,iidim) - XTB(1,iidim)
  508. enddo
  509. xddA1 = vbx*dA1(1) + vby*dA1(2) + vbz*dA1(3)
  510. xddA2 = vbx*dA2(1) + vby*dA2(2) + vbz*dA2(3)
  511. xddA3 = vbx*dA3(1) + vby*dA3(2) + vbz*dA3(3)
  512. if(iimpi.ge.12) write(ioimp,*) ' xddA1,2,3=',xddA1,xddA2,xddA3
  513.  
  514. c si tous les points de TA sont en dessous (ou dessus) de TB,
  515. c il n'y a pas d'intersection : on quitte
  516. c if(xddA1.lt.xpetA.and.xddA2.lt.xpetA.and.xddA3.lt.xpetA)
  517. c & goto 1900
  518. c if(xddA1.gt.xpetAn.and.xddA2.gt.xpetAn.and.xddA3.gt.xpetAn)
  519. c & goto 1900
  520. c idem, mais on traite desormais le cas affleurant
  521. if(xddA1.lt.XPLANn.and.xddA2.lt.XPLANn.and.xddA3.lt.XPLANn)
  522. & goto 1900
  523. if(xddA1.gt.XPLAN.and.xddA2.gt.XPLAN.and.xddA3.gt.XPLAN)
  524. & goto 1900
  525.  
  526. c combien de points de TA dans le plan de TB ? => iAinB
  527. iAinB = 0
  528. if(abs(xddA1).le.XPLAN) iAinB=iAinB+1
  529. if(abs(xddA2).le.XPLAN) iAinB=iAinB+2
  530. if(abs(xddA3).le.XPLAN) iAinB=iAinB+4
  531.  
  532. c si tous les points de TA sont dans le plan de TB,
  533. c-------- => on doit traiter le cas coplanaire ------------------------
  534. if (iAinB.eq.7.or.iBinA.eq.7) then
  535. write(ioimp,*) 'CAS COPLANAIRE EN COURS DE DEVELOPPEMENT'
  536. goto 1900
  537.  
  538. elseif(SELF) then
  539. c si self et 2 points en commun et pas coplanaire => on saute
  540. ichev = 0
  541. do INOD2=1,NBNN2
  542. IPB2 = IPT2.NUM(INOD2,IB)
  543. if(IPB2.eq.ITA(1)) ichev=ichev+1
  544. if(IPB2.eq.ITA(2)) ichev=ichev+1
  545. if(IPB2.eq.ITA(3)) ichev=ichev+1
  546. enddo
  547. if(ichev.eq.2) goto 1900
  548. endif
  549. c-------- fin du cas coplanaire---------------------------------------
  550.  
  551.  
  552. c si 1 seul point de TA est dans le plan de TB
  553. c et les autres dont du meme coté,
  554. c => on ne traite pas cette intersection ponctuelle : on saute
  555. proA12 = xddA1*xddA2
  556. proA23 = xddA2*xddA3
  557. proA31 = xddA3*xddA1
  558. if(iAinB.eq.1.and.proA23.gt.0.d0) goto 1900
  559. if(iAinB.eq.2.and.proA31.gt.0.d0) goto 1900
  560. if(iAinB.eq.4.and.proA12.gt.0.d0) goto 1900
  561.  
  562. c si 2 points de TA sont dans le plan de TB
  563. c ou si on a une traversée franche => on traite
  564.  
  565. C-----------------------------------------------------------------------
  566. C VECTEUR DIRECTEUR DE LA DROITE INTERSECTANT LES 2 PLANS
  567. dx = vay*vbz - vaz*vby
  568. dy = vaz*vbx - vax*vbz
  569. dz = vax*vby - vay*vbx
  570. if(iimpi.ge.12) write(ioimp,*) ' directrice D=',dx,dy,dz
  571.  
  572.  
  573. C-----------------------------------------------------------------------
  574. C EVENTUELLE PERMUTATION DES NOEUDS DE TA
  575. c afin que xddA1 soit de signe opposé à celui de xddA2 et xddA3
  576. c on fait attention aux cas ou on a des noeuds de A dans plan TB
  577. c noeuds A1 et A2 dans plan TB
  578. if(iAinB.eq.3) then
  579. proA12 = 1.d0
  580. proA23 = -1.d0
  581. c noeuds A1 et A3 dans plan TB
  582. elseif(iAinB.eq.5) then
  583. proA12 = -1.d0
  584. proA23 = -1.d0
  585. c noeuds A2 et A3 dans plan TB
  586. elseif(iAinB.eq.6) then
  587. proA12 = 0.d0
  588. proA23 = 1.d0
  589. else
  590. proA12 = xddA1*xddA2
  591. proA23 = xddA2*xddA3
  592. endif
  593. IF(proA23.ge.0.d0) THEN
  594. c si xddA2 et xddA3 de meme signe, rien a faire
  595. c rem: ge possible (au lieu de gt) car on est sur de traverser franchement
  596. iA1=1
  597. iA2=2
  598. iA3=3
  599. ELSE
  600. c xddA2 et xddA3 de signe oppose : il faut permuter
  601. IF(proA12.le.0.d0) THEN
  602. c xddA2 different des 2 autres : 1<-2 2<-3 3<-1
  603. iA1=2
  604. iA2=3
  605. iA3=1
  606. xddtmp= xddA1
  607. xddA1 = xddA2
  608. xddA2 = xddA3
  609. xddA3 = xddtmp
  610. ELSE
  611. c xddA3 different des 2 autres : 1<-3 2<-1 3<-2
  612. iA1=3
  613. iA2=1
  614. iA3=2
  615. xddtmp= xddA2
  616. xddA2 = xddA1
  617. xddA1 = xddA3
  618. xddA3 = xddtmp
  619. ENDIF
  620. ENDIF
  621. va12x = XTA(iA2,1) - XTA(iA1,1)
  622. va12y = XTA(iA2,2) - XTA(iA1,2)
  623. va12z = XTA(iA2,3) - XTA(iA1,3)
  624. va13x = XTA(iA3,1) - XTA(iA1,1)
  625. va13y = XTA(iA3,2) - XTA(iA1,2)
  626. va13z = XTA(iA3,3) - XTA(iA1,3)
  627. c on recalcule iAinB
  628. iAinB = 0
  629. if(abs(xddA1).le.XPLAN) iAinB=iAinB+1
  630. if(abs(xddA2).le.XPLAN) iAinB=iAinB+2
  631. if(abs(xddA3).le.XPLAN) iAinB=iAinB+4
  632. c ici A1 (=1er noeud du Triangle A) est situé de l'autre coté du plan TB
  633. if(iimpi.ge.12) then
  634. write(ioimp,*) ' noeuds A permutés=',ITA(iA1),ITA(iA2),ITA(iA3)
  635. write(ioimp,*) ' xddA1,2,3=',xddA1,xddA2,xddA3
  636. endif
  637.  
  638. C-----------------------------------------------------------------------
  639. C EVENTUELLE PERMUTATION DES NOEUDS DE TB
  640. c afin que xddB1 soit de signe opposé à celui de xddB2 et xddB3
  641. c on fait attention aux cas ou on a des noeuds de A dans plan TB
  642. if(iBinA.eq.3) then
  643. proB12 = 1.d0
  644. proB23 = -1.d0
  645. elseif(iBinA.eq.5) then
  646. proB12 = -1.d0
  647. proB23 = -1.d0
  648. elseif(iBinA.eq.6) then
  649. proB12 = 0.d0
  650. proB23 = 1.d0
  651. else
  652. proB12 = xddB1*xddB2
  653. proB23 = xddB2*xddB3
  654. endif
  655. IF(proB23.ge.0.d0) THEN
  656. c si xddB2 et xddB3 de meme signe, rien a faire
  657. c rem: ge possible (au lieu de gt) car on est sur de traverser franchement
  658. iB1=1
  659. iB2=2
  660. iB3=3
  661. ELSE
  662. c xddB2 et xddB3 de signe oppose : il faut permuter
  663. IF(proB12.le.0.d0) THEN
  664. c xddB2 different des 2 autres : 1<-2 2<-3 3<-1
  665. iB1=2
  666. iB2=3
  667. iB3=1
  668. xddtmp= xddB1
  669. xddB1 = xddB2
  670. xddB2 = xddB3
  671. xddB3 = xddtmp
  672. ELSE
  673. c xddB3 different des 2 autres : 1<-3 2<-1 3<-2
  674. iB1=3
  675. iB2=1
  676. iB3=2
  677. xddtmp= xddB2
  678. xddB2 = xddB1
  679. xddB1 = xddB3
  680. xddB3 = xddtmp
  681. ENDIF
  682. ENDIF
  683. vb12x = XTB(iB2,1) - XTB(iB1,1)
  684. vb12y = XTB(iB2,2) - XTB(iB1,2)
  685. vb12z = XTB(iB2,3) - XTB(iB1,3)
  686. vb13x = XTB(iB3,1) - XTB(iB1,1)
  687. vb13y = XTB(iB3,2) - XTB(iB1,2)
  688. vb13z = XTB(iB3,3) - XTB(iB1,3)
  689. c on recalcule iBinA
  690. iBinA = 0
  691. if(abs(xddB1).le.XPLAN) iBinA=iBinA+1
  692. if(abs(xddB2).le.XPLAN) iBinA=iBinA+2
  693. if(abs(xddB3).le.XPLAN) iBinA=iBinA+4
  694. c ici B1 (=1er noeud du Triangle B) est situé de l'autre coté du plan TB
  695. if(iimpi.ge.12) then
  696. write(ioimp,*) ' noeuds B permutés=',ITB(iB1),ITB(iB2),ITB(iB3)
  697. write(ioimp,*) ' xddB1,2,3=',xddB1,xddB2,xddB3
  698. endif
  699.  
  700.  
  701. C-----------------------------------------------------------------------
  702. C CALCUL DES POINTS P2 et P3
  703.  
  704. c P2 = intersection arete 1-2 de TA avec plan TB
  705. rap2 = abs(xddA1/(xddA1-xddA2))
  706. xP2(1) = XTA(iA1,1) + rap2*va12x
  707. xP2(2) = XTA(iA1,2) + rap2*va12y
  708. xP2(3) = XTA(iA1,3) + rap2*va12z
  709. c P3 = intersection arete 1-3 de TA avec plan TB
  710. rap3 = abs(xddA1/(xddA1-xddA3))
  711. xP3(1) = XTA(iA1,1) + rap3*va13x
  712. xP3(2) = XTA(iA1,2) + rap3*va13y
  713. xP3(3) = XTA(iA1,3) + rap3*va13z
  714.  
  715. c pour y repondre il faut calculer les points Q2 et Q3
  716. if(iimpi.ge.12) write(ioimp,*) ' P2=',(xP2(iou),iou=1,3)
  717. if(iimpi.ge.12) write(ioimp,*) ' P3=',(xP3(iou),iou=1,3)
  718.  
  719. c le segment [P1-P2] est il inclut dans TB ?
  720. C-----------------------------------------------------------------------
  721. C CALCUL DES POINTS Q2 et Q3
  722.  
  723. c Q2 = intersection arete 1-2 de TB avec plan TA
  724. rap2 = abs(xddB1/(xddB1-xddB2))
  725. xQ2(1) = XTB(iB1,1) + rap2*vb12x
  726. xQ2(2) = XTB(iB1,2) + rap2*vb12y
  727. xQ2(3) = XTB(iB1,3) + rap2*vb12z
  728. c Q3 = intersection arete 1-3 de TB avec plan TB
  729. rap3 = abs(xddB1/(xddB1-xddB3))
  730. xQ3(1) = XTB(iB1,1) + rap3*vb13x
  731. xQ3(2) = XTB(iB1,2) + rap3*vb13y
  732. xQ3(3) = XTB(iB1,3) + rap3*vb13z
  733.  
  734. if(iimpi.ge.12) write(ioimp,*) ' Q2=',(xQ2(iou),iou=1,3)
  735. if(iimpi.ge.12) write(ioimp,*) ' Q3=',(xQ3(iou),iou=1,3)
  736.  
  737.  
  738. C-----------------------------------------------------------------------
  739. C IDENTIFICATION EVENTUELLE DES POINTS P2 P3 Q2 ET Q3
  740. c calcul de la distance de P2 aux noeuds A1 et A2
  741. IP2=0
  742. dP2A2 = sqrt((xP2(1)-XTA(IA2,1))**2
  743. & + (xP2(2)-XTA(IA2,2))**2
  744. & + (xP2(3)-XTA(IA2,3))**2)
  745. if (dP2A2.le.XELIM2) then
  746. IP2 =ITA(IA2)
  747. xP2(1)=XTA(IA2,1)
  748. xP2(2)=XTA(IA2,2)
  749. xP2(3)=XTA(IA2,3)
  750. else
  751. dP2A1 = sqrt((xP2(1)-XTA(IA1,1))**2
  752. & + (xP2(2)-XTA(IA1,2))**2
  753. & + (xP2(3)-XTA(IA1,3))**2)
  754. if (dP2A1.le.XELIM2) then
  755. IP2 =ITA(IA1)
  756. xP2(1)=XTA(IA1,1)
  757. xP2(2)=XTA(IA1,2)
  758. xP2(3)=XTA(IA1,3)
  759. endif
  760. endif
  761. c calcul de la distance de P3 aux noeuds A1 et A3
  762. IP3=0
  763. dP3A3 = sqrt((xP3(1)-XTA(IA3,1))**2
  764. & + (xP3(2)-XTA(IA3,2))**2
  765. & + (xP3(3)-XTA(IA3,3))**2)
  766. if (dP3A3.le.XELIM2) then
  767. IP3 =ITA(IA3)
  768. xP3(1)=XTA(IA3,1)
  769. xP3(2)=XTA(IA3,2)
  770. xP3(3)=XTA(IA3,3)
  771. else
  772. dP3A1 = sqrt((xP3(1)-XTA(IA1,1))**2
  773. & + (xP3(2)-XTA(IA1,2))**2
  774. & + (xP3(3)-XTA(IA1,3))**2)
  775. if (dP3A1.le.XELIM2) then
  776. IP3 =ITA(IA1)
  777. xP3(1)=XTA(IA1,1)
  778. xP3(2)=XTA(IA1,2)
  779. xP3(3)=XTA(IA1,3)
  780. endif
  781. endif
  782.  
  783. c calcul de la distance de Q2 aux noeuds B1 et B2
  784. IQ2=0
  785. dQ2B2 = sqrt((xQ2(1)-XTB(IB2,1))**2
  786. & + (xQ2(2)-XTB(IB2,2))**2
  787. & + (xQ2(3)-XTB(IB2,3))**2)
  788. if (dQ2B2.le.XELIM2) then
  789. IQ2 = ITB(IB2)
  790. xQ2(1)=XTB(IB2,1)
  791. xQ2(2)=XTB(IB2,2)
  792. xQ2(3)=XTB(IB2,3)
  793. else
  794. dQ2B1 = sqrt((xQ2(1)-XTB(IB1,1))**2
  795. & + (xQ2(2)-XTB(IB1,2))**2
  796. & + (xQ2(3)-XTB(IB1,3))**2)
  797. if (dQ2B1.le.XELIM2) then
  798. IQ2 = ITB(IB1)
  799. xQ2(1)=XTB(IB1,1)
  800. xQ2(2)=XTB(IB1,2)
  801. xQ2(3)=XTB(IB1,3)
  802. endif
  803. endif
  804. c calcul de la distance de Q3 aux noeuds B1 et B3
  805. IQ3=0
  806. dQ3B3 = sqrt((xQ3(1)-XTB(IB3,1))**2
  807. & + (xQ3(2)-XTB(IB3,2))**2
  808. & + (xQ3(3)-XTB(IB3,3))**2)
  809. if (dQ3B3.le.XELIM2) then
  810. IQ3 =ITB(IB3)
  811. xQ3(1)=XTB(IB3,1)
  812. xQ3(2)=XTB(IB3,2)
  813. xQ3(3)=XTB(IB3,3)
  814. else
  815. dQ3B1 = sqrt((xQ3(1)-XTB(IB1,1))**2
  816. & + (xQ3(2)-XTB(IB1,2))**2
  817. & + (xQ3(3)-XTB(IB1,3))**2)
  818. if (dQ3B1.le.XELIM2) then
  819. IQ3 =ITB(IB1)
  820. xQ3(1)=XTB(IB1,1)
  821. xQ3(2)=XTB(IB1,2)
  822. xQ3(3)=XTB(IB1,3)
  823. endif
  824. endif
  825. c rem bp : cette identification déplace la position de P2 P3 Q2 et Q3
  826. c vers les A1 A2 etc...
  827.  
  828.  
  829. C-----------------------------------------------------------------------
  830. C ABCISSE s SUR D DES POINTS P2, P3, Q2 et Q3
  831. sP2 = xP2(1)*dx + xP2(2)*dy + xP2(3)*dz
  832. sP3 = xP3(1)*dx + xP3(2)*dy + xP3(3)*dz
  833. sQ2 = xQ2(1)*dx + xQ2(2)*dy + xQ2(3)*dz
  834. sQ3 = xQ3(1)*dx + xQ3(2)*dy + xQ3(3)*dz
  835. if(iimpi.ge.10) then
  836. WRITE(IOIMP,*) 'identif #P2,3=',IP2,IP3,' Q2,3=',IQ2,IQ3
  837. WRITE(IOIMP,*) 'abscisses P2,3=',sP2,sP3,' Q2,3=',sQ2,sQ3
  838. endif
  839.  
  840. c si les segments [P2-P3] et [Q2-Q3] ne se chevauchent pas, on quitte
  841. if(sP2.lt.sP3) then
  842. sminP = sP2
  843. smaxP = sP3
  844. else
  845. sminP = sP3
  846. smaxP = sP2
  847. endif
  848. if(sQ2.lt.sQ3) then
  849. sminQ = sQ2
  850. smaxQ = sQ3
  851. else
  852. sminQ = sQ3
  853. smaxQ = sQ2
  854. endif
  855.  
  856. c si pas de chevauchement (exclusion), on saute
  857. if (sminP.ge.smaxQ.or.smaxP.le.sminQ) goto 1900
  858.  
  859.  
  860. C-----------------------------------------------------------------------
  861. C PREPARATION DES MAILLAGES DE TA et TB + SEGMENT D'INTERSECTION
  862. c write(6,*) 'place dispo dan ipt1=',IPT1.NUM(/2),NBELE1
  863. c write(6,*) 'place dispo dan ipt2=',IPT2.NUM(/2),NBELE2
  864. c Avant tout on vérifie qu'il reste assez de place
  865. if(NBELE1+4.gt.IPT1.NUM(/2)) then
  866. NBNN = NBNN1
  867. NBELEM = NBELE1+nbloc
  868. c segadj,IPT1
  869. segadj,IPT1,WORK1
  870. endif
  871. if(NBELE2+4.gt.IPT2.NUM(/2)) then
  872. NBNN = NBNN2
  873. NBELEM = NBELE2+nbloc
  874. segadj,IPT2,WORK2
  875. endif
  876. if(NBELE3+1.gt.IPT3.NUM(/2)) then
  877. NBNN = NBNN3
  878. NBELEM = NBELE3+nbloc
  879. segadj,IPT3
  880. endif
  881. if(NBELE4+1.gt.IPT4.NUM(/2)) then
  882. NBNN = NBNN3
  883. NBELEM = NBELE4+nbloc
  884. segadj,IPT4
  885. endif
  886.  
  887.  
  888. c si chevauchement exact et 2 points de B dans plan A
  889. c et 2 points de A dans plan B => on saute
  890. if ( (abs(sminP-sminQ).le.XELIM)
  891. & .and.(abs(smaxP-smaxQ).le.XELIM) ) then
  892. c on ajoute les segments a IPT3 et IPT4
  893. if (iAinB.eq.6.and.iBinA.eq.6) then
  894. c +on crée les segments intersection
  895. NBTMP=NBELE3+1
  896. NODE(1)=ITA(IA2)
  897. NODE(2)=ITA(IA3)
  898. call ajouel(NODE,2,2,IPT3,NBTMP)
  899. NBELE3=max(NBELE3,NBTMP)
  900. if(.not.SELF) then
  901. NBTMP=NBELE4+1
  902. NODE(1)=ITB(IB2)
  903. NODE(2)=ITB(IB3)
  904. call ajouel(NODE,2,2,IPT4,NBTMP)
  905. NBELE4=max(NBELE4,NBTMP)
  906. endif
  907. if(iimpi.ge.10) then
  908. WRITE(IOIMP,*) 'chevauchement exact IPT3 elem',NBTMP
  909. WRITE(IOIMP,*) IPT3.NUM(1,NBTMP),IPT3.NUM(2,NBTMP)
  910. endif
  911. goto 1900
  912. endif
  913. endif
  914.  
  915. c arrivé ici, on a detecté une intersection entre TA et TB
  916. if(iimpi.ge.10) then
  917. WRITE(IOIMP,*) 'INTERSECTION DETECTÉE entre ',IA,' et ',IB
  918. WRITE(IOIMP,*) 'P=[',sminP,smaxP,'] Q=[',sminQ,smaxQ,']'
  919. endif
  920.  
  921.  
  922.  
  923.  
  924. C-----------------------------------------------------------------------
  925. C PREPARATION DE LA CONSTRUCTION DES FILS DE TA
  926. C + CONSTRUCTION DU SEGMENT D'INTERSECTION
  927.  
  928.  
  929. c --- Cas [P2-P3] inclus dans [Q2-Q3] => segment P2-P3 ---
  930. IF (sminQ.lt.(sminP+XELIM).and.(smaxP-XELIM).lt.smaxQ) THEN
  931.  
  932. c -on est sur le bord du triangle => on ne decoupe pas A
  933. if(IP2.ne.0.and.IP3.ne.0) then
  934. ICOUPA=0
  935. c -on a 2 points confondus => TA decoupe en 2
  936. elseif(IP2.ne.0) then
  937. call AJOUPO(xP3,IPT3,XELIM,IP3)
  938. ICOUPA=23
  939. elseif(IP3.ne.0) then
  940. call AJOUPO(xP2,IPT3,XELIM,IP2)
  941. ICOUPA=22
  942. c -on est dans le cas general => TA decoupe en 3
  943. else
  944. c on ajoute P2 et P3 a XCOOR
  945. call AJOUPO(xP2,IPT3,XELIM,IP2)
  946. call AJOUPO(xP3,IPT3,XELIM,IP3)
  947. ICOUPA=30
  948. endif
  949. c +on crée le segment intersection
  950. NBTMP=NBELE3+1
  951. NODE(1)=IP2
  952. NODE(2)=IP3
  953. call ajouel(NODE,2,2,IPT3,NBTMP)
  954. NBELE3=max(NBELE3,NBTMP)
  955.  
  956. c --- Cas [P2-P3] exclus de [Q2-Q3] => segment Q2-Q3 ---
  957. ELSEIF(sminQ.gt.(sminP+XELIM).and.(smaxP-XELIM).gt.smaxQ)THEN
  958.  
  959. c -on est sur le bord du triangle => TA decoupe en 3
  960. if(IP2.ne.0.and.IP3.ne.0) then
  961. ICOUPA=31
  962. c -on a 2 points confondus => TA decoupe en 5
  963. else
  964. if(abs(xddA2).ge.abs(xddA3)) then
  965. ICOUPA=52
  966. else
  967. ICOUPA=53
  968. endif
  969. endif
  970. call AJOUPO(xQ2,IPT3,XELIM,IQQ2)
  971. c cas particulier IQ2 = IQ3
  972. if(IQ2.eq.IQ3) then
  973. IQQ3=IQQ2
  974. else
  975. call AJOUPO(xQ3,IPT3,XELIM,IQQ3)
  976. endif
  977. c petite astuce pour garder l orientation
  978. if ( (sP2.lt.sP3.and.sQ2.lt.sQ3)
  979. & .or.(sP2.gt.sP3.and.sQ2.gt.sQ3) ) then
  980. IQQ2=IQQ2
  981. IQQ3=IQQ3
  982. else
  983. IQQtmp=IQQ3
  984. IQQ3=IQQ2
  985. IQQ2=IQQtmp
  986. endif
  987. c +on crée le segment intersection
  988. NBTMP=NBELE3+1
  989. NODE(1)=IQQ2
  990. NODE(2)=IQQ3
  991. call ajouel(NODE,2,2,IPT3,NBTMP)
  992. NBELE3=max(NBELE3,NBTMP)
  993. c endif
  994.  
  995. c --- Cas chevauchement selon [Q2 P3] ---
  996. ELSEIF (((sP2+XELIM).le.sQ2.and.sQ2.le.(sP3-XELIM)
  997. & .and.sQ2.le.sQ3)
  998. & .OR.((sP2-XELIM).ge.sQ2.and.sQ2.ge.(sP3+XELIM)
  999. & .and.sQ2.ge.sQ3)) THEN
  1000.  
  1001. c -on a a minima 2 points confondus
  1002. if(IP3.ne.0) then
  1003. c -on est sur le bord du triangle => TA decoupe en 2
  1004. if(IP2.ne.0) then
  1005. ICOUPA=21
  1006. c -on a 2 points confondus => TA decoupe en 3
  1007. else
  1008. ICOUPA=39
  1009. endif
  1010. call AJOUPO(xQ2,IPT3,XELIM,IQQ)
  1011. c -TA decoupe en 4
  1012. else
  1013. call AJOUPO(xP3,IPT3,XELIM,IP3)
  1014. call AJOUPO(xQ2,IPT3,XELIM,IQQ)
  1015. ICOUPA=43
  1016. endif
  1017. c +on crée le segment intersection
  1018. NBTMP=NBELE3+1
  1019. NODE(1)=IQQ
  1020. NODE(2)=IP3
  1021. call ajouel(NODE,2,2,IPT3,NBTMP)
  1022. NBELE3=max(NBELE3,NBTMP)
  1023. c endif
  1024.  
  1025. c --- Cas chevauchement selon [Q3 P3] ---
  1026. ELSEIF (((sP2+XELIM).le.sQ3.and.sQ3.le.(sP3-XELIM)
  1027. & .and.sQ3.le.sQ2)
  1028. & .OR.((sP2-XELIM).ge.sQ3.and.sQ3.ge.(sP3+XELIM)
  1029. & .and.sQ3.ge.sQ2)) THEN
  1030.  
  1031. c -on a a minima 2 points confondus
  1032. if(IP3.ne.0) then
  1033. c -on est sur le bord du triangle => TA decoupe en 2
  1034. if(IP2.ne.0) then
  1035. ICOUPA=21
  1036. c -on a 2 points confondus => TA decoupe en 3
  1037. else
  1038. ICOUPA=39
  1039. endif
  1040. call AJOUPO(xQ3,IPT3,XELIM,IQQ)
  1041. c -TA decoupe en 4
  1042. else
  1043. call AJOUPO(xP3,IPT3,XELIM,IP3)
  1044. call AJOUPO(xQ3,IPT3,XELIM,IQQ)
  1045. ICOUPA=43
  1046. endif
  1047. c +on crée le segment intersection
  1048. NBTMP=NBELE3+1
  1049. NODE(1)=IQQ
  1050. NODE(2)=IP3
  1051. call ajouel(NODE,2,2,IPT3,NBTMP)
  1052. NBELE3=max(NBELE3,NBTMP)
  1053.  
  1054. c --- Cas chevauchement selon [Q2 P2] ---
  1055. ELSEIF (((sP3+XELIM).le.sQ2.and.sQ2.le.(sP2-XELIM)
  1056. & .and.sQ2.le.sQ3)
  1057. & .OR.((sP3-XELIM).ge.sQ2.and.sQ2.ge.(sP2+XELIM)
  1058. & .and.sQ2.ge.sQ3)) THEN
  1059.  
  1060. c -on a a minima 2 points confondus
  1061. if(IP2.ne.0) then
  1062. c -on est sur le bord du triangle => TA decoupe en 2
  1063. if(IP3.ne.0) then
  1064. ICOUPA=21
  1065. c -on a 2 points confondus => TA decoupe en 3
  1066. else
  1067. ICOUPA=39
  1068. endif
  1069. call AJOUPO(xQ2,IPT3,XELIM,IQQ)
  1070. c -TA decoupe en 4
  1071. else
  1072. call AJOUPO(xP2,IPT3,XELIM,IP2)
  1073. call AJOUPO(xQ2,IPT3,XELIM,IQQ)
  1074. ICOUPA=42
  1075. endif
  1076. c +on crée le segment intersection
  1077. NBTMP=NBELE3+1
  1078. NODE(1)=IQQ
  1079. NODE(2)=IP2
  1080. call ajouel(NODE,2,2,IPT3,NBTMP)
  1081. NBELE3=max(NBELE3,NBTMP)
  1082. c endif
  1083.  
  1084. c --- Cas chevauchement selon [Q3 P2] ---
  1085. ELSEIF (((sP3+XELIM).le.sQ3.and.sQ3.le.(sP2-XELIM)
  1086. & .and.sQ3.le.sQ2)
  1087. & .OR.((sP3-XELIM).ge.sQ3.and.sQ3.ge.(sP2+XELIM)
  1088. & .and.sQ3.ge.sQ2)) THEN
  1089.  
  1090. c -on a a minima 2 points confondus
  1091. if(IP2.ne.0) then
  1092. c -on est sur le bord du triangle => TA decoupe en 2
  1093. if(IP3.ne.0) then
  1094. ICOUPA=21
  1095. c -on a 2 points confondus => TA decoupe en 3
  1096. else
  1097. ICOUPA=39
  1098. endif
  1099. call AJOUPO(xQ3,IPT3,XELIM,IQQ)
  1100. c -TA decoupe en 4
  1101. else
  1102. call AJOUPO(xP2,IPT3,XELIM,IP2)
  1103. call AJOUPO(xQ3,IPT3,XELIM,IQQ)
  1104. ICOUPA=42
  1105. endif
  1106. c +on crée le segment intersection
  1107. NBTMP=NBELE3+1
  1108. NODE(1)=IP2
  1109. NODE(2)=IQQ
  1110. call ajouel(NODE,2,2,IPT3,NBTMP)
  1111. NBELE3=max(NBELE3,NBTMP)
  1112. c endif
  1113.  
  1114. ENDIF
  1115. if(iimpi.ge.10) then
  1116. write(ioimp,fmt="('IPT3(1et2',I5,')=',I5,I5)")
  1117. & NBELE3,IPT3.NUM(1,NBELE3),IPT3.NUM(2,NBELE3)
  1118. endif
  1119.  
  1120.  
  1121. C-----------------------------------------------------------------------
  1122. C PREPARATION DE LA CONSTRUCTION DES FILS DE TB
  1123. C + CONSTRUCTION DU SEGMENT D'INTERSECTION
  1124.  
  1125. c --- Cas [Q2-Q3] inclus dans [P2-P3] => segment Q2-Q3 ---
  1126. IF (sminP.lt.(sminQ+XELIM).and.(smaxQ-XELIM).lt.smaxP) THEN
  1127.  
  1128. c -on est sur le bord du triangle => on ne decoupe pas B
  1129. if(IQ2.ne.0.and.IQ3.ne.0) then
  1130. ICOUPB=0
  1131. c -on a 2 points confondus => TB decoupe en 2
  1132. elseif(IQ2.ne.0) then
  1133. call AJOUPO(xQ3,IPT4,XELIM,IQ3)
  1134. ICOUPB=23
  1135. elseif(IQ3.ne.0) then
  1136. call AJOUPO(xQ2,IPT4,XELIM,IQ2)
  1137. IQ3 = ITB(IB3)
  1138. ICOUPB=22
  1139. c -on est dans le cas general => TB decoupe en 3
  1140. else
  1141. c on ajoute Q2 et Q3 a XCOOR
  1142. call AJOUPO(xQ2,IPT4,XELIM,IQ2)
  1143. call AJOUPO(xQ3,IPT4,XELIM,IQ3)
  1144. ICOUPB=30
  1145. endif
  1146. c +on crée le segment intersection
  1147. if(.not.SELF) then
  1148. NBTMP=NBELE4+1
  1149. NODE(1)=IQ2
  1150. NODE(2)=IQ3
  1151. call ajouel(NODE,2,2,IPT4,NBTMP)
  1152. NBELE4=max(NBELE4,NBTMP)
  1153. endif
  1154.  
  1155. c --- Cas [Q2-Q3] exclus de [P2-P3] => segment P2-P3 ---
  1156. ELSEIF(sminP.gt.(sminQ+XELIM).and.(smaxQ-XELIM).gt.smaxP)THEN
  1157.  
  1158. c -on est sur le bord du triangle => TB decoupe en 3
  1159. if(IQ2.ne.0.and.IQ3.ne.0) then
  1160. ICOUPB=31
  1161. c -on a 2 points confondus => TB decoupe en 5
  1162. else
  1163. if(abs(xddB2).ge.abs(xddB3)) then
  1164. ICOUPB=52
  1165. else
  1166. ICOUPB=53
  1167. endif
  1168. endif
  1169. call AJOUPO(xP2,IPT4,XELIM,IPP2)
  1170. c cas particulier IP2 = IP3
  1171. if(IP2.eq.IP3) then
  1172. IPP3=IPP2
  1173. else
  1174. call AJOUPO(xP3,IPT4,XELIM,IPP3)
  1175. endif
  1176. c petite astuce pour garder l orientation
  1177. if ( (sQ2.lt.sQ3.and.sP2.lt.sP3)
  1178. & .or.(sQ2.gt.sQ3.and.sP2.gt.sP3) ) then
  1179. IPP2=IPP2
  1180. IPP3=IPP3
  1181. else
  1182. IPPtmp=IPP3
  1183. IPP3=IPP2
  1184. IPP2=IPPtmp
  1185. endif
  1186. c +on crée le segment intersection
  1187. if(.not.SELF) then
  1188. NBTMP=NBELE4+1
  1189. NODE(1)=IPP2
  1190. NODE(2)=IPP3
  1191. call ajouel(NODE,2,2,IPT4,NBTMP)
  1192. NBELE4=max(NBELE4,NBTMP)
  1193. endif
  1194.  
  1195. c --- Cas chevauchement selon [P2 Q3] ---
  1196. ELSEIF (((sQ2+XELIM).le.sP2.and.sP2.le.(sQ3-XELIM)
  1197. & .and.sP2.le.sP3)
  1198. & .OR.((sQ2-XELIM).ge.sP2.and.sP2.ge.(sQ3+XELIM)
  1199. & .and.sP2.ge.sP3)) THEN
  1200.  
  1201. c -on a a minima 2 points confondus
  1202. if(IQ3.ne.0) then
  1203. IQ3 = ITB(IB3)
  1204. c -on est sur le bord du triangle => TB decoupe en 2
  1205. if(IQ2.ne.0) then
  1206. ICOUPB=21
  1207. c -on a 2 points confondus => TB decoupe en 3
  1208. else
  1209. ICOUPB=39
  1210. endif
  1211. call AJOUPO(xP2,IPT4,XELIM,IPP)
  1212. c -TB decoupe en 4
  1213. else
  1214. call AJOUPO(xQ3,IPT4,XELIM,IQ3)
  1215. call AJOUPO(xP2,IPT4,XELIM,IPP)
  1216. ICOUPB=43
  1217. endif
  1218. c +on crée le segment intersection
  1219. if(.not.SELF) then
  1220. NBTMP=NBELE4+1
  1221. NODE(1)=IPP
  1222. NODE(2)=IQ3
  1223. call ajouel(NODE,2,2,IPT4,NBTMP)
  1224. NBELE4=max(NBELE4,NBTMP)
  1225. endif
  1226.  
  1227. c --- Cas chevauchement selon [P3 Q3] ---
  1228. ELSEIF (((sQ2+XELIM).le.sP3.and.sP3.le.(sQ3-XELIM)
  1229. & .and.sP3.le.sP2)
  1230. & .OR.((sQ2-XELIM).ge.sP3.and.sP3.ge.(sQ3+XELIM)
  1231. & .and.sP3.ge.sP2)) THEN
  1232.  
  1233. c -on a a minima 2 points confondus
  1234. if(IQ3.ne.0) then
  1235. IQ3 = ITB(IB3)
  1236. c -on est sur le bord du triangle => TB decoupe en 2
  1237. if(IQ2.ne.0) then
  1238. ICOUPB=21
  1239. c -on a 2 points confondus => TB decoupe en 3
  1240. else
  1241. ICOUPB=39
  1242. endif
  1243. call AJOUPO(xP3,IPT4,XELIM,IPP)
  1244. c -TB decoupe en 4
  1245. else
  1246. call AJOUPO(xQ3,IPT4,XELIM,IQ3)
  1247. call AJOUPO(xP3,IPT4,XELIM,IPP)
  1248. ICOUPB=43
  1249. endif
  1250. c +on crée le segment intersection
  1251. if(.not.SELF) then
  1252. NBTMP=NBELE4+1
  1253. NODE(1)=IPP
  1254. NODE(2)=IQ3
  1255. call ajouel(NODE,2,2,IPT4,NBTMP)
  1256. NBELE4=max(NBELE4,NBTMP)
  1257. endif
  1258.  
  1259. c --- Cas chevauchement selon [P2 Q2] ---
  1260. ELSEIF (((sQ3+XELIM).le.sP2.and.sP2.le.(sQ2-XELIM)
  1261. & .and.sP2.le.sP3)
  1262. & .OR.((sQ3-XELIM).ge.sP2.and.sP2.ge.(sQ2+XELIM)
  1263. & .and.sP2.ge.sP3)) THEN
  1264.  
  1265. c -on a a minima 2 points confondus
  1266. if(IQ2.ne.0) then
  1267. c -on est sur le bord du triangle => TB decoupe en 2
  1268. if(IQ3.ne.0) then
  1269. ICOUPB=21
  1270. c -on a 2 points confondus => TB decoupe en 3
  1271. else
  1272. ICOUPB=39
  1273. endif
  1274. call AJOUPO(xP2,IPT4,XELIM,IPP)
  1275. c -TB decoupe en 4
  1276. else
  1277. call AJOUPO(xQ2,IPT4,XELIM,IQ2)
  1278. call AJOUPO(xP2,IPT4,XELIM,IPP)
  1279. ICOUPB=42
  1280. endif
  1281. c +on crée le segment intersection
  1282. if(.not.SELF) then
  1283. NBTMP=NBELE4+1
  1284. NODE(1)=IQ2
  1285. NODE(2)=IPP
  1286. call ajouel(NODE,2,2,IPT4,NBTMP)
  1287. NBELE4=max(NBELE4,NBTMP)
  1288. endif
  1289.  
  1290. c --- Cas chevauchement selon [P3 Q2] ---
  1291. ELSEIF (((sQ3+XELIM).le.sP3.and.sP3.le.(sQ2-XELIM)
  1292. & .and.sP3.le.sP2)
  1293. & .OR.((sQ3-XELIM).ge.sP3.and.sP3.ge.(sQ2+XELIM)
  1294. & .and.sP3.ge.sP2)) THEN
  1295.  
  1296. c -on a a minima 2 points confondus
  1297. if(IQ2.ne.0) then
  1298. c -on est sur le bord du triangle => TB decoupe en 2
  1299. if(IQ3.ne.0) then
  1300. ICOUPB=21
  1301. c -on a 2 points confondus => TB decoupe en 3
  1302. else
  1303. ICOUPB=39
  1304. endif
  1305. call AJOUPO(xP3,IPT4,XELIM,IPP)
  1306. c -TB decoupe en 4
  1307. else
  1308. call AJOUPO(xQ2,IPT4,XELIM,IQ2)
  1309. call AJOUPO(xP3,IPT4,XELIM,IPP)
  1310. ICOUPB=42
  1311. endif
  1312. c +on crée le segment intersection
  1313. if(.not.SELF) then
  1314. NBTMP=NBELE4+1
  1315. NODE(1)=IQ2
  1316. NODE(2)=IPP
  1317. call ajouel(NODE,2,2,IPT4,NBTMP)
  1318. NBELE4=max(NBELE4,NBTMP)
  1319. endif
  1320.  
  1321. ENDIF
  1322.  
  1323.  
  1324. if(ICOUPA.eq.0.and.ICOUPB.eq.0) goto 1900
  1325. C-----------------------------------------------------------------------
  1326. C CREATION DU SEGMENT INTERSECTION (fait + haut)
  1327. c CREATION DES TRIANGLES DECOUPES ISSUS DE TA
  1328. IF (ICOUPA.EQ.21) THEN
  1329. c on decoupe TA en 2 triangles (rem: ITA(IA3)=IP3)
  1330. c -1er sous triangle de TA
  1331. IPT1.NUM(1,IA)=ITA(IA1)
  1332. IPT1.NUM(2,IA)=ITA(IA2)
  1333. IPT1.NUM(3,IA)=IQQ
  1334. c -2eme sous triangle de TA
  1335. IPT1.NUM(1, NBELE1+1)=IQQ
  1336. IPT1.NUM(2, NBELE1+1)=ITA(IA3)
  1337. IPT1.NUM(3, NBELE1+1)=ITA(IA1)
  1338. IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1339. NBELE1 = NBELE1+1
  1340. ELSEIF(ICOUPA.EQ.22) THEN
  1341. c on decoupe TA en 2 triangles (rem: ITA(IA3)=IP3)
  1342. c -1er sous triangle de TA
  1343. IPT1.NUM(1,IA)=ITA(IA1)
  1344. IPT1.NUM(2,IA)=IP2
  1345. IPT1.NUM(3,IA)=ITA(IA3)
  1346. c -2eme sous triangle de TA
  1347. IPT1.NUM(1, NBELE1+1)=IP2
  1348. IPT1.NUM(2, NBELE1+1)=ITA(IA2)
  1349. IPT1.NUM(3, NBELE1+1)=ITA(IA3)
  1350. IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1351. NBELE1 = NBELE1+1
  1352. ELSEIF(ICOUPA.EQ.23) THEN
  1353. c on decoupe TA en 2 triangles (rem: ITA(IA2)=IP2)
  1354. c -1er sous triangle de TA
  1355. IPT1.NUM(1,IA)=ITA(IA1)
  1356. IPT1.NUM(2,IA)=ITA(IA2)
  1357. IPT1.NUM(3,IA)=IP3
  1358. c -2eme sous triangle de TA
  1359. IPT1.NUM(1, NBELE1+1)=ITA(IA2)
  1360. IPT1.NUM(2, NBELE1+1)=ITA(IA3)
  1361. IPT1.NUM(3, NBELE1+1)=IP3
  1362. IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1363. NBELE1 = NBELE1+1
  1364. ELSEIF(ICOUPA.EQ.30) THEN
  1365. c on decoupe TA en 3 triangles
  1366. c -1er sous triangle de TA
  1367. IPT1.NUM(1,IA)=ITA(IA1)
  1368. IPT1.NUM(2,IA)=IP2
  1369. IPT1.NUM(3,IA)=IP3
  1370. c -2eme sous triangle de TA
  1371. IPT1.NUM(1, NBELE1+1)=IP2
  1372. IPT1.NUM(2, NBELE1+1)=ITA(IA2)
  1373. IPT1.NUM(3, NBELE1+1)=ITA(IA3)
  1374. IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1375. c -3eme sous triangle de TA
  1376. IPT1.NUM(1, NBELE1+2)=ITA(IA3)
  1377. IPT1.NUM(2, NBELE1+2)=IP3
  1378. IPT1.NUM(3, NBELE1+2)=IP2
  1379. IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1380. NBELE1 = NBELE1+2
  1381. ELSEIF(ICOUPA.EQ.31) THEN
  1382. c on decoupe TA en 3 triangles
  1383. c -1er sous triangle de TA
  1384. IPT1.NUM(1,IA)=ITA(IA1)
  1385. IPT1.NUM(2,IA)=ITA(IA2)
  1386. IPT1.NUM(3,IA)=IQQ2
  1387. c -2eme sous triangle de TA
  1388. IPT1.NUM(1, NBELE1+1)=ITA(IA1)
  1389. IPT1.NUM(2, NBELE1+1)=IQQ2
  1390. IPT1.NUM(3, NBELE1+1)=IQQ3
  1391. IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1392. c -3eme sous triangle de TA
  1393. IPT1.NUM(1, NBELE1+2)=ITA(IA1)
  1394. IPT1.NUM(2, NBELE1+2)=IQQ3
  1395. IPT1.NUM(3, NBELE1+2)=ITA(IA3)
  1396. IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1397. NBELE1 = NBELE1+2
  1398. ELSEIF(ICOUPA.EQ.39) THEN
  1399. c on decoupe TA en 3 triangles
  1400. c -1er sous triangle de TA
  1401. IPT1.NUM(1,IA)=ITA(IA1)
  1402. IPT1.NUM(2,IA)=ITA(IA2)
  1403. IPT1.NUM(3,IA)=IQQ
  1404. c -2eme sous triangle de TA
  1405. IPT1.NUM(1, NBELE1+1)=ITA(IA2)
  1406. IPT1.NUM(2, NBELE1+1)=IQQ
  1407. IPT1.NUM(3, NBELE1+1)=ITA(IA3)
  1408. IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1409. c -3eme sous triangle de TA
  1410. IPT1.NUM(1, NBELE1+2)=ITA(IA3)
  1411. IPT1.NUM(2, NBELE1+2)=IQQ
  1412. IPT1.NUM(3, NBELE1+2)=ITA(IA1)
  1413. IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1414. NBELE1 = NBELE1+2
  1415. ELSEIF(ICOUPA.EQ.42) THEN
  1416. c on decoupe TA en 4 triangles (arete 1-2 coupee)
  1417. c -1er sous triangle de TA
  1418. IPT1.NUM(1,IA)=ITA(IA1)
  1419. IPT1.NUM(2,IA)=IP2
  1420. IPT1.NUM(3,IA)=IQQ
  1421. c -2eme sous triangle de TA
  1422. IPT1.NUM(1, NBELE1+1)=IP2
  1423. IPT1.NUM(2, NBELE1+1)=ITA(IA2)
  1424. IPT1.NUM(3, NBELE1+1)=IQQ
  1425. IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1426. c -3eme sous triangle de TA
  1427. IPT1.NUM(1, NBELE1+2)=ITA(IA2)
  1428. IPT1.NUM(2, NBELE1+2)=ITA(IA3)
  1429. IPT1.NUM(3, NBELE1+2)=IQQ
  1430. IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1431. c -4eme sous triangle de TA
  1432. IPT1.NUM(1, NBELE1+3)=ITA(IA3)
  1433. IPT1.NUM(2, NBELE1+3)=ITA(IA1)
  1434. IPT1.NUM(3, NBELE1+3)=IQQ
  1435. IPT1.ICOLOR(NBELE1+3)=IPT1.ICOLOR(IA)
  1436. NBELE1 = NBELE1+3
  1437. ELSEIF(ICOUPA.EQ.43) THEN
  1438. c on decoupe TA en 4 triangles (arete 1-3 coupee)
  1439. c -1er sous triangle de TA
  1440. IPT1.NUM(1,IA)=ITA(IA1)
  1441. IPT1.NUM(2,IA)=ITA(IA2)
  1442. IPT1.NUM(3,IA)=IQQ
  1443. c -2eme sous triangle de TA
  1444. IPT1.NUM(1, NBELE1+1)=ITA(IA2)
  1445. IPT1.NUM(2, NBELE1+1)=ITA(IA3)
  1446. IPT1.NUM(3, NBELE1+1)=IQQ
  1447. IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1448. c -3eme sous triangle de TA
  1449. IPT1.NUM(1, NBELE1+2)=ITA(IA3)
  1450. IPT1.NUM(2, NBELE1+2)=IP3
  1451. IPT1.NUM(3, NBELE1+2)=IQQ
  1452. IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1453. c -4eme sous triangle de TA
  1454. IPT1.NUM(1, NBELE1+3)=IP3
  1455. IPT1.NUM(2, NBELE1+3)=ITA(IA1)
  1456. IPT1.NUM(3, NBELE1+3)=IQQ
  1457. IPT1.ICOLOR(NBELE1+3)=IPT1.ICOLOR(IA)
  1458. NBELE1 = NBELE1+3
  1459. ELSEIF(ICOUPA.EQ.52) THEN
  1460. c on decoupe TA en 5 triangles
  1461. c -1er sous triangle de TA
  1462. IPT1.NUM(1,IA)=ITA(IA1)
  1463. IPT1.NUM(2,IA)=IQQ2
  1464. IPT1.NUM(3,IA)=IQQ3
  1465. c -2eme sous triangle de TA
  1466. IPT1.NUM(1,NBELE1+1)=ITA(IA1)
  1467. IPT1.NUM(2,NBELE1+1)=ITA(IA2)
  1468. IPT1.NUM(3,NBELE1+1)=IQQ2
  1469. IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1470. c -3eme sous triangle de TA
  1471. IPT1.NUM(1,NBELE1+2)=ITA(IA2)
  1472. IPT1.NUM(2,NBELE1+2)=IQQ3
  1473. IPT1.NUM(3,NBELE1+2)=IQQ2
  1474. IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1475. c -4eme sous triangle de TA
  1476. IPT1.NUM(1,NBELE1+3)=ITA(IA2)
  1477. IPT1.NUM(2,NBELE1+3)=ITA(IA3)
  1478. IPT1.NUM(3,NBELE1+3)=IQQ3
  1479. IPT1.ICOLOR(NBELE1+3)=IPT1.ICOLOR(IA)
  1480. c -5eme sous triangle de TA
  1481. IPT1.NUM(1,NBELE1+4)=ITA(IA3)
  1482. IPT1.NUM(2,NBELE1+4)=ITA(IA1)
  1483. IPT1.NUM(3,NBELE1+4)=IQQ3
  1484. IPT1.ICOLOR(NBELE1+4)=IPT1.ICOLOR(IA)
  1485. NBELE1 = NBELE1+4
  1486. ELSEIF(ICOUPA.EQ.53) THEN
  1487. c on decoupe TA en 5 triangles
  1488. c -1er sous triangle de TA
  1489. IPT1.NUM(1,IA)=ITA(IA1)
  1490. IPT1.NUM(2,IA)=IQQ2
  1491. IPT1.NUM(3,IA)=IQQ3
  1492. c -2eme sous triangle de TA
  1493. IPT1.NUM(1,NBELE1+1)=ITA(IA1)
  1494. IPT1.NUM(2,NBELE1+1)=ITA(IA2)
  1495. IPT1.NUM(3,NBELE1+1)=IQQ2
  1496. IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1497. c -3eme sous triangle de TA
  1498. IPT1.NUM(1,NBELE1+2)=ITA(IA2)
  1499. IPT1.NUM(2,NBELE1+2)=ITA(IA3)
  1500. IPT1.NUM(3,NBELE1+2)=IQQ2
  1501. IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1502. c -4eme sous triangle de TA
  1503. IPT1.NUM(1,NBELE1+3)=ITA(IA3)
  1504. IPT1.NUM(2,NBELE1+3)=IQQ3
  1505. IPT1.NUM(3,NBELE1+3)=IQQ2
  1506. IPT1.ICOLOR(NBELE1+3)=IPT1.ICOLOR(IA)
  1507. c -5eme sous triangle de TA
  1508. IPT1.NUM(1,NBELE1+4)=ITA(IA3)
  1509. IPT1.NUM(2,NBELE1+4)=ITA(IA1)
  1510. IPT1.NUM(3,NBELE1+4)=IQQ3
  1511. IPT1.ICOLOR(NBELE1+4)=IPT1.ICOLOR(IA)
  1512. NBELE1 = NBELE1+4
  1513. ENDIF
  1514.  
  1515. c CREATION DES TRIANGLES DECOUPES ISSUS DE TB
  1516. if(SELF) NBELE2=NBELE1
  1517. IF (ICOUPB.EQ.21) THEN
  1518. c on decoupe TB en 2 triangles
  1519. c -1er sous triangle de TB
  1520. IPT2.NUM(1,IB)=ITB(IB1)
  1521. IPT2.NUM(2,IB)=ITB(IB2)
  1522. IPT2.NUM(3,IB)=IPP
  1523. c -2eme sous triangle de TB
  1524. IPT2.NUM(1, NBELE2+1)=IPP
  1525. IPT2.NUM(2, NBELE2+1)=ITB(IB3)
  1526. IPT2.NUM(3, NBELE2+1)=ITB(IB1)
  1527. IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1528. NBELE2 = NBELE2+1
  1529. ELSEIF(ICOUPB.EQ.22) THEN
  1530. c on decoupe TB en 2 triangles (rem: ITB(IB3)=IQ3)
  1531. c -1er sous triangle de TB
  1532. IPT2.NUM(1,IB)=ITB(IB1)
  1533. IPT2.NUM(2,IB)=IQ2
  1534. IPT2.NUM(3,IB)=ITB(IB3)
  1535. c -2eme sous triangle de TB
  1536. IPT2.NUM(1, NBELE2+1)=IQ2
  1537. IPT2.NUM(2, NBELE2+1)=ITB(IB2)
  1538. IPT2.NUM(3, NBELE2+1)=ITB(IB3)
  1539. IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1540. NBELE2 = NBELE2+1
  1541. ELSEIF(ICOUPB.EQ.23) THEN
  1542. c on decoupe TB en 2 triangles (rem: ITB(IB2)=IQ2)
  1543. c -1er sous triangle de TB
  1544. IPT2.NUM(1,IB)=ITB(IB1)
  1545. IPT2.NUM(2,IB)=ITB(IB2)
  1546. IPT2.NUM(3,IB)=IQ3
  1547. c -2eme sous triangle de TB
  1548. IPT2.NUM(1, NBELE2+1)=ITB(IB2)
  1549. IPT2.NUM(2, NBELE2+1)=ITB(IB3)
  1550. IPT2.NUM(3, NBELE2+1)=IQ3
  1551. IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1552. NBELE2 = NBELE2+1
  1553. ELSEIF(ICOUPB.EQ.30) THEN
  1554. c on decoupe TB en 3 triangles
  1555. c -1er sous triangle de TB
  1556. IPT2.NUM(1,IB)=ITB(IB1)
  1557. IPT2.NUM(2,IB)=IQ2
  1558. IPT2.NUM(3,IB)=IQ3
  1559. c -2eme sous triangle de TB
  1560. IPT2.NUM(1,NBELE2+1)=IQ2
  1561. IPT2.NUM(2,NBELE2+1)=ITB(IB2)
  1562. IPT2.NUM(3,NBELE2+1)=ITB(IB3)
  1563. IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1564. c -3eme sous triangle de TB
  1565. IPT2.NUM(1,NBELE2+2)=IQ2
  1566. IPT2.NUM(2,NBELE2+2)=ITB(IB3)
  1567. IPT2.NUM(3,NBELE2+2)=IQ3
  1568. IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1569. NBELE2=NBELE2+2
  1570. ELSEIF(ICOUPB.EQ.31) THEN
  1571. c on decoupe TB en 3 triangles
  1572. c -1er sous triangle de TB
  1573. IPT2.NUM(1,IB)=ITB(IB1)
  1574. IPT2.NUM(2,IB)=ITB(IB2)
  1575. IPT2.NUM(3,IB)=IPP2
  1576. c -2eme sous triangle de TB
  1577. IPT2.NUM(1, NBELE2+1)=ITB(IB1)
  1578. IPT2.NUM(2, NBELE2+1)=IPP2
  1579. IPT2.NUM(3, NBELE2+1)=IPP3
  1580. IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1581. c -3eme sous triangle de TB
  1582. IPT2.NUM(1, NBELE2+2)=ITB(IB1)
  1583. IPT2.NUM(2, NBELE2+2)=IPP3
  1584. IPT2.NUM(3, NBELE2+2)=ITB(IB3)
  1585. IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1586. NBELE2=NBELE2+2
  1587. ELSEIF(ICOUPB.EQ.39) THEN
  1588. c on decoupe TB en 3 triangles
  1589. c -1er sous triangle de TB
  1590. IPT2.NUM(1,IB)=ITB(IB1)
  1591. IPT2.NUM(2,IB)=ITB(IB2)
  1592. IPT2.NUM(3,IB)=IPP
  1593. c -2eme sous triangle de TB
  1594. IPT2.NUM(1, NBELE2+1)=ITB(IB2)
  1595. IPT2.NUM(2, NBELE2+1)=IPP
  1596. IPT2.NUM(3, NBELE2+1)=ITB(IB3)
  1597. IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1598. c -3eme sous triangle de TB
  1599. IPT2.NUM(1, NBELE2+2)=ITB(IB3)
  1600. IPT2.NUM(2, NBELE2+2)=IPP
  1601. IPT2.NUM(3, NBELE2+2)=ITB(IB1)
  1602. IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1603. NBELE2=NBELE2+2
  1604. ELSEIF(ICOUPB.EQ.42) THEN
  1605. c on decoupe TB en 4 triangles (arete 1-2 coupee)
  1606. c -1er sous triangle de TB
  1607. IPT2.NUM(1,IB)=ITB(IB1)
  1608. IPT2.NUM(2,IB)=IQ2
  1609. IPT2.NUM(3,IB)=IPP
  1610. c -2eme sous triangle de TB
  1611. IPT2.NUM(1, NBELE2+1)=IQ2
  1612. IPT2.NUM(2, NBELE2+1)=ITB(IB2)
  1613. IPT2.NUM(3, NBELE2+1)=IPP
  1614. IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1615. c -3eme sous triangle de TB
  1616. IPT2.NUM(1, NBELE2+2)=ITB(IB2)
  1617. IPT2.NUM(2, NBELE2+2)=ITB(IB3)
  1618. IPT2.NUM(3, NBELE2+2)=IPP
  1619. IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1620. c -4eme sous triangle de TB
  1621. IPT2.NUM(1, NBELE2+3)=ITB(IB3)
  1622. IPT2.NUM(2, NBELE2+3)=ITB(IB1)
  1623. IPT2.NUM(3, NBELE2+3)=IPP
  1624. IPT2.ICOLOR(NBELE2+3)=IPT2.ICOLOR(IB)
  1625. NBELE2=NBELE2+3
  1626. ELSEIF(ICOUPB.EQ.43) THEN
  1627. c on decoupe TB en 4 triangles (arete 1-3 coupee)
  1628. c -1er sous triangle de TB
  1629. IPT2.NUM(1,IB)=ITB(IB1)
  1630. IPT2.NUM(2,IB)=ITB(IB2)
  1631. IPT2.NUM(3,IB)=IPP
  1632. c -2eme sous triangle de TB
  1633. IPT2.NUM(1, NBELE2+1)=ITB(IB2)
  1634. IPT2.NUM(2, NBELE2+1)=ITB(IB3)
  1635. IPT2.NUM(3, NBELE2+1)=IPP
  1636. IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1637. c -3eme sous triangle de TB
  1638. IPT2.NUM(1, NBELE2+2)=ITB(IB3)
  1639. IPT2.NUM(2, NBELE2+2)=IQ3
  1640. IPT2.NUM(3, NBELE2+2)=IPP
  1641. IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1642. c -4eme sous triangle de TB
  1643. IPT2.NUM(1, NBELE2+3)=IQ3
  1644. IPT2.NUM(2, NBELE2+3)=ITB(IB1)
  1645. IPT2.NUM(3, NBELE2+3)=IPP
  1646. IPT2.ICOLOR(NBELE2+3)=IPT2.ICOLOR(IB)
  1647. NBELE2=NBELE2+3
  1648. ELSEIF(ICOUPB.EQ.52) THEN
  1649. c on decoupe TB en 5 triangles
  1650. c -1er sous triangle de TB
  1651. IPT2.NUM(1,IB)=ITB(IB1)
  1652. IPT2.NUM(2,IB)=IPP2
  1653. IPT2.NUM(3,IB)=IPP3
  1654. c -2eme sous triangle de TB
  1655. IPT2.NUM(1,NBELE2+1)=ITB(IB1)
  1656. IPT2.NUM(2,NBELE2+1)=ITB(IB2)
  1657. IPT2.NUM(3,NBELE2+1)=IPP2
  1658. IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1659. c -3eme sous triangle de TB
  1660. IPT2.NUM(1,NBELE2+2)=ITB(IB2)
  1661. IPT2.NUM(2,NBELE2+2)=IPP3
  1662. IPT2.NUM(3,NBELE2+2)=IPP2
  1663. IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1664. c -4eme sous triangle de TB
  1665. IPT2.NUM(1,NBELE2+3)=ITB(IB2)
  1666. IPT2.NUM(2,NBELE2+3)=ITB(IB3)
  1667. IPT2.NUM(3,NBELE2+3)=IPP3
  1668. IPT2.ICOLOR(NBELE2+3)=IPT2.ICOLOR(IB)
  1669. c -5eme sous triangle de TB
  1670. IPT2.NUM(1,NBELE2+4)=ITB(IB3)
  1671. IPT2.NUM(2,NBELE2+4)=ITB(IB1)
  1672. IPT2.NUM(3,NBELE2+4)=IPP3
  1673. IPT2.ICOLOR(NBELE2+4)=IPT2.ICOLOR(IB)
  1674. NBELE2=NBELE2+4
  1675. ELSEIF(ICOUPB.EQ.53) THEN
  1676. c on decoupe TB en 5 triangles
  1677. c -1er sous triangle de TB
  1678. IPT2.NUM(1,IB)=ITB(IB1)
  1679. IPT2.NUM(2,IB)=IPP2
  1680. IPT2.NUM(3,IB)=IPP3
  1681. c -2eme sous triangle de TB
  1682. IPT2.NUM(1,NBELE2+1)=ITB(IB1)
  1683. IPT2.NUM(2,NBELE2+1)=ITB(IB2)
  1684. IPT2.NUM(3,NBELE2+1)=IPP2
  1685. IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1686. c -3eme sous triangle de TB
  1687. IPT2.NUM(1,NBELE2+2)=ITB(IB2)
  1688. IPT2.NUM(2,NBELE2+2)=ITB(IB3)
  1689. IPT2.NUM(3,NBELE2+2)=IPP2
  1690. IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1691. c -4eme sous triangle de TB
  1692. IPT2.NUM(1,NBELE2+3)=ITB(IB3)
  1693. IPT2.NUM(2,NBELE2+3)=IPP3
  1694. IPT2.NUM(3,NBELE2+3)=IPP2
  1695. IPT2.ICOLOR(NBELE2+3)=IPT2.ICOLOR(IB)
  1696. c -5eme sous triangle de TB
  1697. IPT2.NUM(1,NBELE2+4)=ITB(IB3)
  1698. IPT2.NUM(2,NBELE2+4)=ITB(IB1)
  1699. IPT2.NUM(3,NBELE2+4)=IPP3
  1700. IPT2.ICOLOR(NBELE2+4)=IPT2.ICOLOR(IB)
  1701. NBELE2=NBELE2+4
  1702. ENDIF
  1703. if(SELF) NBELE1=NBELE2
  1704.  
  1705. if(iimpi.ge.10) then
  1706. write(ioimp,*) '------------------------------->',ICOUPA
  1707. write(ioimp,*) '------------------------------->',ICOUPB
  1708. write(ioimp,*) '-----------------------------------------'
  1709. endif
  1710.  
  1711. c Tant qu'il y a une decoupe, on reprend le traitement relatif
  1712. c au 1er element "fils" du triangle decoupe ...
  1713. IBDEB0 = IB+1
  1714. if(IBDEB0.le.NBELE2) GOTO 1100
  1715.  
  1716.  
  1717. c 1900 : label ou l on arrive s il n y a pas eu de decoupe
  1718. 1900 CONTINUE
  1719. IBDEB0=0
  1720. if(IA.ge.NBELE1.and.IB.ge.NBELE2) goto 9000
  1721.  
  1722.  
  1723. 2000 CONTINUE
  1724. C====== fin de BOUCLE SUR LES ELEMENTS B ===============================
  1725.  
  1726.  
  1727. 900 CONTINUE
  1728. if(IA.ge.NBELE1) goto 9000
  1729.  
  1730. 1000 CONTINUE
  1731. C==== fin de BOUCLE SUR LES ELEMENTS A =================================
  1732.  
  1733.  
  1734.  
  1735. 9000 CONTINUE
  1736. C-----------------------------------------------------------------------
  1737. C FIN NORMALE
  1738.  
  1739. c on supprime les segments de travaux locaux
  1740. segsup,WORK1
  1741.  
  1742. c on ajuste IPT1
  1743. NBNN =NBNN1
  1744. NBELEM=NBELE1
  1745. segadj,IPT1
  1746. c on ajuste IPT3
  1747. NBNN = 2
  1748. NBELEM = NBELE3
  1749. segadj,IPT3
  1750. c on desactive
  1751. segdes,IPT1,IPT3
  1752.  
  1753. IF(.not.SELF) THEN
  1754. segsup,WORK2
  1755. c on ajuste IPT2
  1756. NBNN =NBNN2
  1757. NBELEM=NBELE2
  1758. segadj,IPT2
  1759. c on ajuste IPT4
  1760. NBNN = 2
  1761. NBELEM = NBELE4
  1762. segadj,IPT4
  1763. c on desactive
  1764. segdes,IPT2,IPT4
  1765. ENDIF
  1766.  
  1767.  
  1768. C ON REGENERE et on ne garde que les triangles (ITYPEL=4)
  1769. CALL REGENE(IPT1)
  1770. IPT1ok = IPT1
  1771. segact,IPT1
  1772. NBSOU1=IPT1.LISOUS(/1)
  1773. if(NBSOU1.gt.0) then
  1774. do isou1=1,NBSOU1
  1775. MELEME=IPT1.LISOUS(isou1)
  1776. segact,MELEME
  1777. if(ITYPEL.eq.4) then
  1778. IPT1ok=MELEME
  1779. segdes,MELEME
  1780. segsup,IPT1
  1781. goto 9010
  1782. else
  1783. segsup,MELEME
  1784. endif
  1785. segdes,MELEME
  1786. enddo
  1787. endif
  1788. 9010 CONTINUE
  1789.  
  1790. IF(SELF) GOTO 9020
  1791. C ON REGENERE et on ne garde que les triangles (ITYPEL=4)
  1792. CALL REGENE(IPT2)
  1793. IPT2ok = IPT2
  1794. segact,IPT2
  1795. NBSOU1=IPT2.LISOUS(/1)
  1796. if(NBSOU1.gt.0) then
  1797. do isou1=1,NBSOU1
  1798. MELEME=IPT2.LISOUS(isou1)
  1799. segact,MELEME
  1800. if(ITYPEL.eq.4) then
  1801. IPT2ok=MELEME
  1802. segdes,MELEME
  1803. segsup,IPT2
  1804. goto 9020
  1805. else
  1806. segsup,MELEME
  1807. endif
  1808. segdes,MELEME
  1809. enddo
  1810. endif
  1811. 9020 CONTINUE
  1812.  
  1813.  
  1814. c on renvoie
  1815. IOB1 = IPT1ok
  1816. IF(.not.SELF) IOB2 = IPT2ok
  1817. IOB3 = IPT3
  1818. IF(.not.SELF) IOB4 = IPT4
  1819.  
  1820.  
  1821. END
  1822.  
  1823.  
  1824.  
  1825.  
  1826.  
  1827.  

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