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