1 : C PSIP3D SOURCE BP208322 14/01/29 21:15:14 7921 2 : c 3 : SUBROUTINE PSIP3D(ideux,melmai,melfis,melfro,xcrit, 4 : & msoup1,msoup2,msoup3) 5 : 6 : *********************************************************************** 7 : c OPERATEUR : PSIP 8 : 9 : c APPELE PAR PSIPHI DANS LE CAS 3D 10 : c 11 : c CREATION : BP, 14/03/2012 12 : c MODIF : BP, 19/12/2013 plus d'appel au zonage 13 : c 14 : *********************************************************************** 15 : IMPLICIT REAL*8(A-H,O-Z) 16 : IMPLICIT INTEGER(I-N) 17 : -INC CCOPTIO 18 : -INC SMCOORD 19 : -INC SMELEME 20 : -INC SMCHPOI 21 : -INC CCREEL 22 : segment icpr(xcoor(/1)/(idim+1)) 23 : segment xicpr(xcoor(/1)/(idim+1)) 24 : segment isup(npt) 25 : segment xdis(npt) 26 : segment xdmin(npt) 27 : segment xdmin1(npt) 28 : segment rrlim(npt) 29 : segment xdmin2(npt) 30 : segment xdmin3(npt) 31 : segment isupfi(npt) 32 : segment sfis 33 : integer ifis(nbelfi,2) 34 : endsegment 35 : DATA EPSI/1.D-5/ 36 : logical phicor,phi12,phi23,phi31 37 : integer ooperm(7,2) 38 : DATA ooperm/ 2, 4, 6, 1, 3, 5, 7, 39 : $ 4, 1, 5, 2, 6, 3, 7/ 40 : 41 : 42 : *********************************************************************** 43 : * Initialisation, Recup, Activation 44 : *********************************************************************** 45 : if(iimpi.ge.1) then 45 : call gibtem(xkt) 46 : write(ioimp,*)'===== ENTREE DANS PSIP3D ===== | xkt=',xkt 47 : endif 48 : 49 : idim1=idim+1 50 : idbug = 0 51 : 52 : c des chpoints (les msoup sont deja actifs) 53 : if(msoup2.gt.0) mpova2=msoup2.ipoval 54 : if(msoup1.gt.0) mpova1=msoup1.ipoval 55 : if(msoup3.gt.0) mpova3=msoup3.ipoval 56 : c activation des maillages / desactivation 57 : c melmai : fait dans psiphi / ligne 1325 ou 1329 58 : c melfis : fait ligne 64 / ligne 1274 59 : c melfro : fait (si ideux>=1) ligne 78 / ligne 356 60 : ipt3=melmai 61 : c segact,ipt3 62 : npt=ipt3.num(/2) 63 : 64 : c variables parfois testées 65 : imil5=0 66 : 67 : 68 : *********************************************************************** 69 : * Preparation du maillage de la fissure 70 : *********************************************************************** 71 : meleme=melfis 72 : segact meleme 73 : * transfo du front en TRI3 73 : if(itypel.ne.4) call change(melfis,4) 74 : meleme=melfis 75 : segact meleme 76 : c rem: inutile car change laisse melfis actif 77 : nseg=num(/2) 78 : 79 : *------ Cas ou l'on calcule phi ET psi -------------------------------- 80 : if (ideux.ge.1) then 81 : * on rearrange le melfis de facon à prevenir dans ifis(iel,1) quels 82 : * element touchent le front de fissure. Ceux ci sont ordonnées de facon 83 : * que le front de fissure soit du noeud 1 vers le noeud 2 de l'element 84 : ipt2=melfro 85 : segact ipt2 86 : segini icpr 87 : nbelem2 = ipt2.num(/2) 88 : nbnode2 = ipt2.num(/1) 89 : * petite numerotation locale 90 : * icpr(noeud)>0 si le noeud appartient au front de fissure 91 : * =numero du 1er element de melfro ou on a vu ce noeud 92 : * + on en profite pour calculer l'abscisse curviligne du front si ideux.eq.2 93 : * icpr(noeud)<0 si le noeud appartient au reste du contour de la fissure 94 : * =-numero du 1er element de melcon ou on a vu ce noeud 95 : if (ideux.ge.2) then 96 : segini,xicpr 97 : xtau=0.d0 98 : ic=ipt2.num(1,1) 99 : x11= xcoor((ic-1)*idim1 +1) 100 : y11= xcoor((ic-1)*idim1 +2) 101 : z11= xcoor((ic-1)*idim1 +3) 102 : do ia=1,nbelem2 103 : do ib=1,nbnode2 104 : ic=ipt2.num(ib,ia) 105 : if (icpr(ic).eq.0) then 106 : icpr(ic)=ia 107 : x22= xcoor((ic-1)*idim1 +1) 108 : y22= xcoor((ic-1)*idim1 +2) 109 : z22= xcoor((ic-1)*idim1 +3) 110 : dxtau= sqrt( (x22-x11)*(x22-x11)+(y22-y11)*(y22-y11) 111 : $ +(z22-z11)*(z22-z11) ) 112 : xtau=xtau+dxtau 113 : xicpr(ic)=xtau 114 : x11=x22 115 : y11=y22 116 : z11=z22 117 : c if(idebug.ge.1) 118 : c $ write(6,*) 'ia,ib,ic,xicpr(ic)=',ia,ib,ic,xicpr(ic) 119 : endif 120 : enddo 121 : enddo 122 : else 123 : do ia=1,nbelem2 124 : do ib=1,nbnode2 125 : ic=ipt2.num(ib,ia) 126 : if (icpr(ic).eq.0) then 127 : icpr(ic)=ia 128 : endif 129 : enddo 130 : enddo 131 : endif 132 : * bp (23/02/2012) : calcul du contour et fin de remplissage de icpr 132 : CALL ECROBJ('MAILLAGE',melfis) 132 : CALL PRCONT 132 : CALL LIROBJ('MAILLAGE',melcon,1,IRETOU) 133 : IF (IERR.NE.0) RETURN 134 : meleme = melfis 135 : IPT6 = melcon 136 : SEGACT,meleme,IPT6 137 : nbelCON = IPT6.num(/2) 138 : nbnoCON = IPT6.num(/1) 139 : do ia=1,nbelCON 140 : do ib=1,nbnoCON 141 : ic=IPT6.num(ib,ia) 142 : if(icpr(ic).eq.0) icpr(ic)=-ia 143 : enddo 144 : enddo 145 : * nombre de segments appartenant au contour sans le front 146 : nbelem6 = nbelCON - nbelem2 147 : * remplissage effectif de ifis contenu dans le segment sfis 148 : i2firs = ipt2.num(1,1) 149 : i2last = ipt2.num(nbnode2,nbelem2) 150 : ivui2f = 0 151 : ivui2l = 0 152 : c write(6,*) 'i2firs,i2last=',i2firs,i2last 153 : nbelfi=nseg 154 : nbelem=nseg 155 : nbnn=3 156 : nbsous=0 157 : nbref=0 158 : c segini,ipt5=meleme 159 : * bp 20/12/2011 : on propose de remplir ipt5 avec : 160 : c + des elements du front (et du contour) d abord, 161 : c + des elements avec 1 point appartenant au front (et au contour) ensuite, 162 : c + les autres elements enfin. 163 : c bp 06/03/2012 : on propose de remplir ifis avec : 164 : c ifis(iseg,1) =2 si le segment [1-2] du iseg^ieme triangle appartient 165 : c au front de fissure, =3 si en plus c'est l extrémité du front, 166 : c =1 si l un des segments appartient au contour 167 : c =-1 si l'un des noeuds appartient au front 168 : c ifis(iseg,2) =+1+2+4 si les segments [1-2] [2-3] [3-1] du iseg^ieme triangle appartiennent au contour 169 : segini,ipt5 170 : ideb5 = 0 171 : imil5 = nbelem2 172 : ifin5 = nseg+1 173 : segini sfis 174 : 175 : *--------boucle sur les elements de la fissure ---------------------- 176 : do iel=1,nseg 177 : 178 : * io est incrémenté de +1,+2,+4 si node 1,2,3 appartient au front 179 : * ioo est incrémenté de +1,+2,+4 si node 1,2,3 appartient au contour 180 : io=0 181 : ioo=0 182 : * ---boucle sur les noeuds des elements de la fissure ---------- 183 : do ib=1,num(/1) 184 : ia=ib 185 : if(ib.eq.3) ia = 4 186 : ic = num(ib,iel) 187 : if(icpr(ic).gt.0) io=io+ia 188 : if(icpr(ic).lt.0) ioo=ioo+ia 189 : * cas particulier d un segment avec 1 point appartenant au contour et un autre 190 : if(ic.eq.i2firs.or.ic.eq.i2last) ioo=ioo+ia 191 : enddo 192 : * ---fin de boucle sur les noeuds des elements de la fissure --- 193 : 194 : * remplissage de ipt5(i5,...) et ifis(i5,1) 195 : i5=0 196 : iperm=0 197 : * ---cas ou 1 segment appartient au front 198 : if (io.eq.3.or.io.gt.4) then 199 : ideb5=ideb5+1 200 : i5=ideb5 201 : * (eventuelle) permutation des noeuds pour avoir segment 1-2 = front 202 : if (io.eq.3) then 203 : ipt5.num(1,i5)=num(1,iel) 204 : ipt5.num(2,i5)=num(2,iel) 205 : ipt5.num(3,i5)=num(3,iel) 206 : elseif(io.eq.5) then 207 : ipt5.num(1,i5)=num(3,iel) 208 : ipt5.num(2,i5)=num(1,iel) 209 : ipt5.num(3,i5)=num(2,iel) 210 : iperm=1 211 : elseif(io.eq.6) then 212 : ipt5.num(1,i5)=num(2,iel) 213 : ipt5.num(2,i5)=num(3,iel) 214 : ipt5.num(3,i5)=num(1,iel) 215 : iperm=2 216 : elseif(io.eq.7) then 217 : write(ioimp,*) 'ERREUR: ' 218 : $ ,'Plusieurs cotes d un meme triangle du maillage de ' 219 : $ ,'la fissure converti en TRI3 appartiennent au front!' 220 : return 221 : endif 222 : * ifis( ,1)=3 si c'est une extremité 223 : if ((ipt5.num(1,i5)).eq.i2firs) then 224 : ifis(i5,1)=3 225 : ivui2f = ivui2f+1 226 : elseif((ipt5.num(2,i5)).eq.i2last) then 227 : ifis(i5,1)=3 228 : ivui2l = ivui2l+1 229 : else 230 : * ifis( ,1)=2 si c'est un tri3 du front avec des voisins 231 : ifis(i5,1)=2 232 : endif 233 : * ---cas ou 1 (ou 2) segment(s) appartient au contour 234 : elseif(ioo.eq.3.or.ioo.gt.4) then 235 : imil5=imil5+1 236 : i5=imil5 237 : * pas de permutation des noeuds 238 : ipt5.num(1,i5)=num(1,iel) 239 : ipt5.num(2,i5)=num(2,iel) 240 : ipt5.num(3,i5)=num(3,iel) 241 : * ifis( ,1)=1 si un segment appartient au contour 242 : ifis(i5,1)=1 243 : * ---cas ou seulement 1 point appartient au front 244 : elseif(io.eq.1.or.io.eq.2.or.io.eq.4) then 245 : ifin5=ifin5-1 246 : i5=ifin5 247 : * pas de permutation des noeuds 248 : ipt5.num(1,i5)=num(1,iel) 249 : ipt5.num(2,i5)=num(2,iel) 250 : ipt5.num(3,i5)=num(3,iel) 251 : * ifis( ,1)=-1 si 1 seul point du tri3 appartient au front 252 : ifis(i5,1)=-1 253 : * ---autres cas 254 : else 255 : ifin5=ifin5-1 256 : i5=ifin5 257 : * pas de permutation des noeuds 258 : ipt5.num(1,i5)=num(1,iel) 259 : ipt5.num(2,i5)=num(2,iel) 260 : ipt5.num(3,i5)=num(3,iel) 261 : endif 262 : 263 : * remplissage de ifis(i5,2) 264 : * ---uniquement si 1 (ou 2) segment(s) appartient au contour 265 : if(ioo.eq.3.or.ioo.gt.4) then 266 : * attention a l eventuelle permutation effectuee 267 : if(iperm.ne.0) ioo=ooperm(ioo,iperm) 268 : * -ioo=3 => noeuds 1-2 = segment 1 appartient au contour 269 : if (ioo.eq.3) then 270 : ifis(i5,2)=1 271 : * -ioo=5 => noeuds 3-1 = segment 3 appartient au contour 272 : elseif(ioo.eq.5) then 273 : ifis(i5,2)=4 274 : * -ioo=6 => noeuds 2-3 = segment 2 appartient au contour 275 : elseif(ioo.eq.6) then 276 : ifis(i5,2)=2 277 : * -cas ou 3 points appartiennent au contour => travail en + 278 : elseif(ioo.eq.7) then 279 : c write(6,*)'tri3 avec 2 segments dans le contour!',iel,i5 280 : * rem: si les 3 points appartiennent au contour, on n'a pas de segments 281 : * appartenant au front 282 : ii1 = ipt5.num(1,i5) 283 : ii2 = ipt5.num(2,i5) 284 : ii3 = ipt5.num(3,i5) 285 : c il faut trouver quels sont les 2 segments appartenant au contour ! 286 : ia1 = -icpr(ii1) 287 : ia2 = -icpr(ii2) 288 : ia3 = -icpr(ii3) 289 : ia1min = min(ia1,ia2) 290 : ia1min = min(ia1min,ia3) 291 : ia1seg = 0 292 : do 11 ia1=ia1min,nbelCON 293 : ifis1 = 0 294 : do ib=1,nbnoCON 295 : if(ii1.eq.IPT6.num(ib,ia1)) ifis1=ifis1+1 296 : if(ii2.eq.IPT6.num(ib,ia1)) ifis1=ifis1+2 297 : if(ii3.eq.IPT6.num(ib,ia1)) ifis1=ifis1+4 298 : enddo 299 : if(ifis1.le.2) goto 11 300 : if(ifis1.eq.4) goto 11 301 : if(ifis1.eq.3) ifis(i5,2)=ifis(i5,2)+1 302 : if(ifis1.eq.6) ifis(i5,2)=ifis(i5,2)+2 303 : if(ifis1.eq.5) ifis(i5,2)=ifis(i5,2)+4 304 : if(ifis1.ge.7) then 305 : write(ioimp,*) 'tri3 avec 3 segments dans le contour !' 305 : call erreur(26) 306 : endif 307 : if(ifis1.eq.3.or.ifis1.ge.5) ia1seg=ia1seg+1 308 : if(ia1seg.eq.2) goto 12 309 : 11 continue 310 : c on n a pas trouvé les 2 segments appartenant au contour 311 : write(ioimp,*)'on n a pas trouvé les 2 segments', 312 : $ ' du tri3 devant appartenant au contour !' 312 : call erreur(26) 313 : c on a trouvé les 2 segments appartenant au contour 314 : 12 continue 315 : if(ifis(i5,2).lt.3.or.ifis(i5,2).ge.7) then 316 : write(ioimp,*) 'impossible de detecter les 2 segments ', 317 : $ 'appartenant au contour avec les noeuds ',ii1,ii2,ii3 318 : write(ioimp,*) 'ifis(ideb5,2)=',ifis(ideb5,2) 319 : write(ioimp,*) 'IPT6 pour les segments :' 320 : write(ioimp,*) ' ia1',ia1,(IPT6.num(iou,ia1),iou=1,2) 321 : write(ioimp,*) ' ia2',ia2,(IPT6.num(iou,ia2),iou=1,2) 322 : write(ioimp,*) ' ia3',ia3,(IPT6.num(iou,ia3),iou=1,2) 322 : call erreur(26) 323 : return 324 : endif 325 : endif 326 : endif 327 : 328 : 329 : 330 : enddo 331 : * fin de boucle sur les elements de la fissure ----------------- 332 : 333 : 334 : c if ((ifin5-ideb5).ne.1) then 335 : if (ideb5.ne.nbelem2.or.(ifin5-imil5).ne.1) then 336 : write(ioimp,*) 'pb avec le rearrangement de elements de fissure' 337 : write(ioimp,*) ideb5,nbelem2,imil5,ifin5 337 : call erreur(26) 338 : return 339 : endif 340 : * si le front n est pas ferme, il faut avoir detecte ses 2 extremites 341 : * sinon il est probable que la ligne du front melfro ne soit pas 342 : * orientée comme le bord du maillage melfis 343 : if ((ivui2f.ne.1).and.(ivui2l.ne.1)) then 344 : write(ioimp,*) 'ERREUR: Les extremites du front (noeuds', 345 : $ i2firs,i2last,') ont été detectees ',ivui2f, 'et ',ivui2l 346 : $ ,' fois au lieu de 1 fois attendue' 347 : write(ioimp,*) 'Verifiez la coherence avec l orientation' 348 : $ ,'du maillage de la fissure fourni' 349 : $ ,' et INVErser le sens du front le cas échéant.' 350 : c write(6,*) 'front de fissure=' 350 : c CALL ECMAIL(ipt2,JENT2) 351 : c write(6,*) 'surface de fissure=' 351 : c CALL ECMAIL(ipt5,JENT5) 351 : call erreur(26) 352 : return 353 : endif 354 : c segdes,meleme,IPT6 355 : segdes,meleme,IPT6,ipt2 356 : meleme=ipt5 357 : melfis=meleme 358 : segsup icpr 359 : 360 : endif 361 : *------ fin du Cas ou l'on calcule phi ET psi -------------------------- 362 : 363 : * sont actifs: melmai(=struct changée en poi1), ipt5=meleme=melfis, 364 : * sont inactifs: ipt2=melfro 365 : 366 : 367 : 368 : *********************************************************************** 369 : * debut du travail suivant valeur de icle 370 : 371 : if (xcrit.eq.0.d0) then 372 : icle=0 373 : xcrit=1.D10 374 : else 375 : icle=1 376 : * ajout d'un epsilon 1.D-5 de rattrapage 377 : xcrit = 1.00001D0 * xcrit 378 : endif 379 : if(iimpi.ge.1) then 379 : call gibtem(xkt) 380 : write(ioimp,*)'debut zonage si icle=1=',icle,xcrit,'| xkt=',xkt 381 : endif 382 : 383 : melfis=meleme 384 : melpoi=melmai 385 : 386 : *********************************************************************** 387 : * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT 388 : *********************************************************************** 389 : 390 : * pour chaque segment de la fissure (melfis 2eme maillage lu) 391 : * on calcule le champ phi et psi des noeuds autour de l'element 392 : if(iimpi.ge.1) then 392 : call gibtem(xkt) 393 : write(ioimp,*)'debut travail | xkt=',xkt 394 : endif 395 : * 396 : meleme=melfis 397 : c segact,meleme 398 : ipt3=melpoi 399 : c segact,ipt3 400 : npt=ipt3.num(/2) 401 : * attention on initialise isup avec les valeurs inverses,il faudra 402 : * le retourner 403 : segini isup,isupfi,xdis,xdmin,xdmin1,rrlim 404 : if(ideux.ge.1) segini xdmin2,xdmin3 405 : igard=0 406 : igardf=0 407 : 408 : 409 : *===================================================================== 410 : * Boucle sur les elements de la fissure ========================= 411 : 412 : do 400 iseg=1,nseg 413 : 414 : c coordonnes du tri3 415 : ia=num(1,iseg) 416 : x11= xcoor((ia-1)*idim1 +1) 417 : y11= xcoor((ia-1)*idim1 +2) 418 : z11= xcoor((ia-1)*idim1 +3) 419 : ib=num(2,iseg) 420 : x22= xcoor((ib-1)*idim1 +1) 421 : y22= xcoor((ib-1)*idim1 +2) 422 : z22= xcoor((ib-1)*idim1 +3) 423 : ic=num(3,iseg) 424 : x33= xcoor((ic-1)*idim1 +1) 425 : y33= xcoor((ic-1)*idim1 +2) 426 : z33= xcoor((ic-1)*idim1 +3) 427 : if (x11.ge.x22) then 428 : xxmi=min(x22,x33) - xcrit 429 : xxma=max(x11,x33) + xcrit 430 : else 431 : xxmi=min(x11,x33) - xcrit 432 : xxma=max(x22,x33) + xcrit 433 : endif 434 : if (y11.ge.y22) then 435 : yymi=min(y22,y33) - xcrit 436 : yyma=max(y11,y33) + xcrit 437 : else 438 : yymi=min(y11,y33) - xcrit 439 : yyma=max(y22,y33) + xcrit 440 : endif 441 : if (z11.ge.z22) then 442 : zzmi=min(z22,z33) - xcrit 443 : zzma=max(z11,z33) + xcrit 444 : else 445 : zzmi=min(z11,z33) - xcrit 446 : zzma=max(z22,z33) + xcrit 447 : endif 448 : 449 : c vecteurs des aretes 450 : v12x=x22-x11 451 : v12y=y22-y11 452 : v12z=z22-z11 453 : v23x=x33-x22 454 : v23y=y33-y22 455 : v23z=z33-z22 456 : v31x=x11-x33 457 : v31y=y11-y33 458 : v31z=z11-z33 459 : 460 : c vecteur normal 461 : vx=v12y*(z33-z11)-v12z*(y33-y11) 462 : vy=v12z*(x33-x11)-v12x*(z33-z11) 463 : vz=v12x*(y33-y11)-v12y*(x33-x11) 464 : bb=vx*vx + vy*vy+ vz*vz 465 : aa= sqrt(bb) 466 : if (aa.le.XPETIT) then 467 : write(ioimp,*) 'element fissure de longueur nulle ',aa 468 : write(ioimp,*) 'tri3.',iseg,'.',ia,'=',x11,y11,z11 469 : write(ioimp,*) 'tri3.',iseg,'.',ib,'=',x22,y22,z22 470 : write(ioimp,*) 'tri3.',iseg,'.',ic,'=',x33,y33,z33 470 : call erreur(26) 471 : return 472 : endif 473 : sur = bb/2.d0 474 : vx=vx/aa 475 : vy=vy/aa 476 : vz=vz/aa 477 : vxyz = max(abs(vx),abs(vy)) 478 : vxyz = max(vxyz,abs(vz)) 479 : if (vxyz.le.XPETIT) then 480 : write(ioimp,*)'element fissure de longueur nulle',aa,vx,vy,vz 480 : call erreur(26) 481 : return 482 : endif 483 : xsur = 1.D-8*sur 484 : 485 : c vecteurs des aretes (suite) 486 : v12lo= sqrt(v12x*v12x+v12y*v12y+v12z*v12z) 487 : v23lo= sqrt(v23x*v23x+v23y*v23y+v23z*v23z) 488 : v31lo= sqrt(v31x*v31x+v31y*v31y+v31z*v31z) 489 : v12x= v12x/v12lo 490 : v12y= v12y/v12lo 491 : v12z= v12z/v12lo 492 : v23x= v23x/v23lo 493 : v23y= v23y/v23lo 494 : v23z= v23z/v23lo 495 : v31x= v31x/v31lo 496 : v31y= v31y/v31lo 497 : v31z= v31z/v31lo 498 : 499 : if(iimpi.ge.1) 500 : $ write(ioimp,*)'-> triangle ',iseg,'/',nseg,' noeuds ',ia,ib,ic 501 : 502 : 503 : *======================================================================= 504 : * Boucle sur les noeuds ====================================== 505 : do 490 ipt=1,npt 506 : 507 : id=ipt3.num(1,ipt) 508 : dis = xgrand 509 : iphi2 = 0 510 : 511 : *------------ cas ou le noeud structure id = l un des noeuds des tri3 de la fissure 512 : if (ia.eq.id.or.ib.eq.id.or.ic.eq.id) then 513 : 514 : dis=0.D0 515 : iphi2 = 1 516 : 517 : *------------ "autres" cas 518 : else 519 : 520 : * recup des coordonnes du noeud en question 521 : xxp = xcoor((id-1)*idim1 + 1 ) 522 : yyp = xcoor((id-1)*idim1 + 2 ) 523 : zzp = xcoor((id-1)*idim1 + 3 ) 524 : 525 : *bp2013 on saute ce point si trop loin (critere) 526 : if (icle.eq.1) then 527 : if(xxp.lt.xxmi) goto 491 528 : if(xxp.gt.xxma) goto 491 529 : if(yyp.lt.yymi) goto 491 530 : if(yyp.gt.yyma) goto 491 531 : if(zzp.lt.zzmi) goto 491 532 : if(zzp.gt.zzma) goto 491 533 : endif 534 : 535 : d1xp=xxp-x11 536 : d1yp=yyp-y11 537 : d1zp=zzp-z11 538 : d2xp=xxp-x22 539 : d2yp=yyp-y22 540 : d2zp=zzp-z22 541 : d3xp=xxp-x33 542 : d3yp=yyp-y33 543 : d3zp=zzp-z33 544 : * calcul de la distance au triangle 545 : d1= sqrt(d1xp*d1xp + d1yp*d1yp + d1zp*d1zp) 546 : d2= sqrt(d2xp*d2xp + d2yp*d2yp + d2zp*d2zp) 547 : d3= sqrt(d3xp*d3xp + d3yp*d3yp + d3zp*d3zp) 548 : if (d1.lt.d2) then 549 : rr= min(d1,d3) 550 : else 551 : rr= min(d2,d3) 552 : endif 553 : 554 : * test pour eviter plein de calculs => on saute au noeud suivant 555 : if (isup(ipt).eq.0) then 556 : rrlim(ipt) = rr 557 : else 558 : c if(rr.gt.(1.1d0*rrlim(ipt))) goto 491 559 : if(rr.gt.(2.d0*rrlim(ipt))) goto 491 560 : if(rr.lt.rrlim(ipt)) rrlim(ipt)=rr 561 : endif 562 : 563 : * calcul de la distance au plan du triangle 564 : xdd = vx*d1xp + vy*d1yp + vz*d1zp 565 : * determination du pt projetté 566 : xx=xxp-xdd*vx 567 : yy=yyp-xdd*vy 568 : zz=zzp-xdd*vz 569 : * ce pt est-il à l'interieur ou sur un coté du triangle? 570 : d1x=xx-x11 571 : d1y=yy-y11 572 : d1z=zz-z11 573 : d2x=xx-x22 574 : d2y=yy-y22 575 : d2z=zz-z22 576 : d3x=xx-x33 577 : d3y=yy-y33 578 : d3z=zz-z33 579 : pv12x=d1y*d2z-d1z*d2y 580 : pv12y=d1z*d2x-d1x*d2z 581 : pv12z=d1x*d2y-d1y*d2x 582 : pv23x=d2y*d3z-d2z*d3y 583 : pv23y=d2z*d3x-d2x*d3z 584 : pv23z=d2x*d3y-d2y*d3x 585 : pv31x=d3y*d1z-d3z*d1y 586 : pv31y=d3z*d1x-d3x*d1z 587 : pv31z=d3x*d1y-d3y*d1x 588 : ps12v=pv12X*vx+pv12y*vy+pv12z*vz 589 : ps23v=pv23x*vx+pv23y*vy+pv23z*vz 590 : ps31v=pv31x*vx+pv31y*vy+pv31z*vz 591 : 592 : * ce pt est-il à l'interieur ou sur un coté du triangle? 593 : if(ps12v.ge.0..and.ps23v.ge.0..and.ps31v.ge.0.) then 594 : * -> on est a l'interieur 595 : dis = xdd 596 : iphi2 = 1 597 : elseif (abs(ps12v).lt.xsur) then 598 : * -> on est situé sur le segment 12 : est-on a l'interieur? 599 : psc1=v12x*d1x+v12y*d1y+v12z*d1z 600 : psc2=v12x*d2x+v12y*d2y+v12z*d2z 601 : if ( psc1*psc2 .lt.0.d0) then 602 : dis = xdd 603 : iphi2 = 1 604 : endif 605 : elseif(abs(ps23v).lt.xsur)then 606 : * -> on est situé sur le segment 23 : est-on a l'interieur? 607 : psc2=v23x*d2x+v23y*d2y+v23z*d2z 608 : psc3=v23x*d3x+v23y*d3y+v23z*d3z 609 : if ( psc2*psc3 .lt.0.d0) then 610 : dis = xdd 611 : iphi2 = 1 612 : endif 613 : elseif(abs(ps31v).lt.xsur)then 614 : * -> on est situé sur le segment 31 : est-on a l'interieur? 615 : psc3=v31x*d3x+v31y*d3y+v31z*d3z 616 : psc1=v31x*d1x+v31y*d1y+v31z*d1z 617 : if ( psc3*psc1 .lt.0.d0) then 618 : dis = xdd 619 : iphi2 = 1 620 : endif 621 : endif 622 : 623 : * on est a l'exterieur : il faut calculer dis et verifier xdd 624 : if (dis.eq.xgrand) then 625 : iphi2 = 10 626 : 627 : * si on est a l'ext du tri3 et en avant du front+contour on change v 628 : 629 : * -si on a affaire a un triangle du front (ou du contour) : 630 : * est on en avant du front (ou du contour)? -> phicor 631 : phicor=.false. 632 : phi12=.false. 633 : phi23=.false. 634 : phi31=.false. 635 : if (ideux.ge.1) then 636 : if (ifis(iseg,1).ge.1) then 637 : * -Cas ou le segment [1-2] = front ou contour 638 : if(ifis(iseg,1).ge.2.or.ifis(iseg,2).eq.1.or. 639 : $ ifis(iseg,2).eq.3.or.ifis(iseg,2).eq.5) then 640 : c phi12=(ps12v.lt.0.d0) 641 : phi12=(ps12v.lt.(-xsur)) 642 : if(phi12) then 643 : if(ifis(iseg,1).ge.2) then 644 : iphi2=iphi2+5 645 : else 646 : iphi2=iphi2+1 647 : endif 648 : endif 649 : endif 650 : * -Cas ou le segment [2-3] = contour 651 : if(ifis(iseg,2).eq.2.or. 652 : $ ifis(iseg,2).eq.3.or.ifis(iseg,2).eq.6) then 653 : c phi23=(ps23v.lt.0.d0) 654 : phi23=(ps23v.lt.(-xsur)) 655 : if(phi23) iphi2=iphi2+1 656 : endif 657 : * -Cas ou le segment [3-1] = contour 658 : if(ifis(iseg,2).eq.4.or. 659 : $ ifis(iseg,2).eq.5.or.ifis(iseg,2).eq.6) then 660 : c phi31=(ps31v.lt.0.d0) 661 : phi31=(ps31v.lt.(-xsur)) 662 : if(phi31) iphi2=iphi2+1 663 : endif 664 : * -Union des test 665 : phicor=(ideux.ge.1).and.(phi12.or.phi23.or.phi31) 666 : endif 667 : endif 668 : * si phicor, il faudra alors "corriger" phi 669 : * (et prendre xdd au lieu de dis) 670 : 671 : psc1=0.d0 672 : psc2=0.d0 673 : psc3=0.d0 674 : c dis=sign(rr,xdd) 675 : * bp (05/03/2012) : on remplace ci dessus par formule + precise 676 : * --on teste autour du noeud 1-- 677 : if (rr.eq.d1) then 678 : if (ps12v.lt.0.d0.and.ps31v.lt.0d0) then 679 : * on est definitivement + proche du noeud 1 680 : dis = sign(rr,xdd) 681 : else 682 : if (ps12v.lt.0.d0) then 683 : * on est du coté du segment 12 684 : psc1 = v12x*d1x + v12y*d1y + v12z*d1z 685 : elseif (ps31v.lt.0.d0) then 686 : * on est du coté du segment 13 687 : psc1 = -1.d0*(v31x*d1x + v31y*d1y + v31z*d1z) 688 : else 689 : * on ne devrait pas passer par la 690 : psc1=0.d0 691 : endif 692 : if (psc1.le.0d0) then 693 : dis = sign(rr,xdd) 694 : else 695 : dis = sign(sqrt(rr*rr-psc1*psc1),xdd) 696 : endif 697 : endif 698 : * --on teste autour du noeud 2-- 699 : elseif (rr.eq.d2) then 700 : if (ps12v.lt.0.d0.and.ps23v.lt.0d0) then 701 : * on est definitivement + proche du noeud 2 702 : dis = sign(rr,xdd) 703 : else 704 : if (ps12v.lt.0.d0) then 705 : * on est du coté du segment 21 706 : psc2 = -1.d0*(v12x*d2x + v12y*d2y + v12z*d2z) 707 : elseif (ps31v.lt.0.d0) then 708 : * on est du coté du segment 23 709 : psc2 = v23x*d2x + v23y*d2y + v23z*d2z 710 : else 711 : * on ne devrait pas passer par la 712 : psc2=0.d0 713 : endif 714 : if (psc2.le.0d0) then 715 : dis = sign(rr,xdd) 716 : else 717 : dis = sign(sqrt(rr*rr-psc2*psc2),xdd) 718 : endif 719 : endif 720 : * --on teste autour du noeud 3-- 721 : elseif (rr.eq.d3) then 722 : if (ps31v.lt.0.d0.and.ps23v.lt.0d0) then 723 : * on est definitivement + proche du noeud 2 724 : dis = sign(rr,xdd) 725 : else 726 : if (ps31v.lt.0.d0) then 727 : * on est du coté du segment 31 728 : psc3 = v31x*d3x + v31y*d3y + v31z*d3z 729 : elseif (ps31v.lt.0.d0) then 730 : * on est du coté du segment 32 731 : psc3 = -1.d0*(v23x*d3x + v23y*d3y + v23z*d3z) 732 : else 733 : * on ne devrait pas passer par la 734 : psc3=0.d0 735 : endif 736 : if (psc3.le.0d0) then 737 : dis = sign(rr,xdd) 738 : else 739 : dis = sign(sqrt(rr*rr-psc3*psc3),xdd) 740 : endif 741 : endif 742 : endif 743 : 744 : endif 745 : * fin du cas ou on est a l'exterieur 746 : 747 : endif 748 : *------------ fin des "autres" cas (noeud id different de ceux du tri3) 749 : 750 : * -on a calculé (xdd et) dis: s'agit il de phi? 751 : c on prend xdd plutot que dis pour definir Phi dans le cas ou : 752 : c -on a l option DEUX ou TROI 753 : c -le tri3 utilisé (le + proche) appartient au front ou au contour 754 : c -le point se projette en avant du front ou du contour 755 : * iphi2 qu on va eventuellement ranger dans isup(ipt) vaut: 756 : * 1 si on se projette a l int du triangle 757 : * 10 si on se projette a l ext du triangle 758 : * 11 si en + on sort du contour 759 : * 15 si en + on est en avant du front 760 : * 16 si en + on est en avant du front et on sort du contour 761 : 762 : if(id.eq.idbug) then 763 : c if(phicor) then 764 : write(6,*) 'noeud:',id,' (',ipt,') tri3',iseg,':',ia,ib,ic 765 : write(6,*) 'iphi2,isup(ipt),dis,xdd',iphi2,isup(ipt),dis,xdd 766 : write(6,*) 'phi12,phicor=',phi12,phicor 767 : endif 768 : 769 : * -c'est la premiere fois que l'on voit ce point 770 : if (isup(ipt).eq.0) then 771 : if(phicor) then 772 : mpova2.vpocha(ipt,1)=xdd 773 : else 774 : mpova2.vpocha(ipt,1)=dis 775 : endif 776 : xdmin(ipt) =dis 777 : xdmin1(ipt)=xdd 778 : isup(ipt)=iphi2 779 : iphi2=100 780 : igard=igard+1 781 : else 782 : * --on est (franchement) + proche que les precedents essais 783 : if (abs(dis).lt.0.999999*abs(xdmin(ipt))) then 784 : if(phicor) then 785 : mpova2.vpocha(ipt,1)=xdd 786 : else 787 : mpova2.vpocha(ipt,1)=dis 788 : endif 789 : xdmin(ipt) =dis 790 : xdmin1(ipt)=xdd 791 : isup(ipt)=iphi2 792 : iphi2=100 793 : * --on est (quasiment) aussi proche que le meilleur essai 794 : elseif (abs(dis).lt.1.000001*abs(xdmin(ipt))) then 795 : * -on a 2 tri3 du front+contour a egalite 796 : * et l'un des deux a trouvé le point en avant du front : 797 : if ((iseg.le.imil5).and. 798 : $ (phicor.or.isup(ipt).gt.10)) then 799 : c strategie : xdd= max(xdd1 xdd2) 800 : if (abs(xdd).gt.abs(xdmin1(ipt))) then 801 : mpova2.vpocha(ipt,1)=xdd 802 : xdmin(ipt) =dis 803 : xdmin1(ipt)=xdd 804 : c isup(ipt)=iphi2 805 : if(iphi2.gt.isup(ipt)) isup(ipt)=iphi2 806 : iphi2=100 807 : c cas ou on avait pas vu qu il fallait prendre xdd 808 : elseif(isup(ipt).le.10) then 809 : c equivalent à : 810 : c elseif(mpova2.vpocha(ipt,1).ne.xdmin1(ipt)) then 811 : mpova2.vpocha(ipt,1)=xdmin1(ipt) 812 : isup(ipt)=iphi2 813 : iphi2=100 814 : endif 815 : 816 : * -on avait 1 tri3 du front+contour et ce tri3 est a egalite 817 : * et le 1er avait trouvé le point en avant du front : 818 : * on laisse le xdd calculé depuis le 1er tri3 819 : elseif ((iseg.gt.imil5).and. 820 : $ (isup(ipt).gt.10)) then 821 : goto 491 822 : * -quasi-egalite entre 2 concurrents mais signes opposés : 823 : else 824 : ss1=dis*xdmin(ipt) 825 : if (ss1.le.0.d0) then 826 : * on choisit le dis associé au xdd le + marqué 827 : if (abs(xdd).gt.abs(xdmin1(ipt))) then 828 : mpova2.vpocha(ipt,1)=dis 829 : xdmin(ipt) =dis 830 : xdmin1(ipt)=xdd 831 : c isup(ipt)=iphi2 832 : if(iphi2.gt.isup(ipt)) isup(ipt)=iphi2 833 : iphi2=100 834 : endif 835 : endif 836 : endif 837 : endif 838 : * --fin de la distinction + proche / aussi proche 839 : endif 840 : * -fin de la distinction 1ere fois ou pas 841 : 842 : if(id.eq.idbug) then 843 : write(6,*)'iphi2,xdmin,xdmin1,phi',iphi2, 844 : $ xdmin(ipt),xdmin1(ipt),mpova2.vpocha(ipt,1) 845 : endif 846 : if(iimpi.ge.2) then 847 : write(ioimp,*) 'point',id,ipt,'dis,xdd,PHI,isup=', 848 : $ dis,xdd,mpova2.vpocha(ipt,1),isup(ipt) 849 : endif 850 : 851 : 852 : 491 continue 853 : * le tri3 n'est pas un tri3 du front => on saute 854 : *TODO : test a verifier peut etre .eq. et .eq.-1 aussi a ne pas sauter? 855 : if (ideux.lt.1) goto 490 856 : if (ifis(iseg,1).le.1) goto 490 857 : 858 : c ------ Calcul de la 2eme Level Set ------- 859 : xdd1 = 0.D0 860 : 861 : * -cas ou le noeud structure id = l un des noeuds des tri3 de la fissure 862 : if (ia.eq.id.or.ib.eq.id) then 863 : dis1=0.D0 864 : xd2=0.D0 865 : xd3=0.D0 866 : if(ideux.ge.2) xtau=xicpr(id) 867 : * -autres cas 868 : else 869 : * on recupere ce qu on a pas recupéré lors du calcul de phi 870 : if (ic.eq.id) then 871 : xxp = xcoor((id-1)*idim1 + 1 ) 872 : yyp = xcoor((id-1)*idim1 + 2 ) 873 : zzp = xcoor((id-1)*idim1 + 3 ) 874 : endif 875 : * calcul du vecteur vn qui est le sens de propagation 876 : vnx=v12y*vz-v12z*vy 877 : vny=v12z*vx-v12x*vz 878 : vnz=v12x*vy-v12y*vx 879 : * il faut orienter vnx pour que le vecteur soit dans le sens 880 : * de la propagation 881 : qq= vnx*(x33-x11)+vny*(y33-y11)+vnz*(z33-z11) 882 : if (qq.ge.0.D0) then 883 : vnx=-vnx 884 : vny=-vny 885 : vnz=-vnz 886 : endif 887 : * normalement vn est normé on le verifie 888 : qq= sqrt(vnx*vnx + vny*vny + vnz*vnz) 889 : if (abs(qq-1.d0).gt.1.d-2) then 890 : write(ioimp,*) ' pb dans le calcul de vn' 890 : call erreur(26) 891 : return 892 : elseif(abs(qq-1.).gt.1.d-6) then 893 : vnx=vnx/qq 894 : vny=vny/qq 895 : vnz=vnz/qq 896 : endif 897 : * on calcule la distance au plan defini par segment 12 et vecteur vx,vy,vz 898 : * et le point (xx yy zz ) projeté du point sur ce plan 899 : xdd1=(xxp-x11)*vnx+(yyp-y11)*vny+(zzp-z11)*vnz 900 : xx=xxp-xdd1*vnx 901 : yy=yyp-xdd1*vny 902 : zz=zzp-xdd1*vnz 903 : psc1=(xx-x11)*v12x+(yy-y11)*v12y+(zz-z11)*v12z 904 : psc2=(xx-x22)*v12x+(yy-y22)*v12y+(zz-z22)*v12z 905 : if ((psc1*psc2).le.0.d0) then 906 : C le point se projete sur 12 907 : dis1=xdd1 908 : pp1x=x11+(psc1*v12x) 909 : pp1y=y11+(psc1*v12y) 910 : pp1z=z11+(psc1*v12z) 911 : xd2 = 0.D0 912 : xd3 = sqrt((xxp-pp1x)**2+(yyp-pp1y)**2+(zzp-pp1z)**2) 913 : if (ideux.ge.2) then 914 : xtau1=xicpr(ia) 915 : xtau2=xicpr(ib) 916 : stau = psc1 / v12lo 917 : xtau=stau*xtau2+(1.d0-stau)*xtau1 918 : endif 919 : else 920 : if (psc1.ge.0.D0) then 921 : pp2x=xx+((v12lo-psc1)*v12x) 922 : pp2y=yy+((v12lo-psc1)*v12y) 923 : pp2z=zz+((v12lo-psc1)*v12z) 924 : xd2 = abs(v12lo-psc1) 925 : xd3 = sqrt((xxp-x22)**2+(yyp-y22)**2+(zzp-z22)**2) 926 : c if(ideux.ge.2) xtau=xicpr(ib) 927 : if(ideux.ge.2) then 928 : if(ib.eq.i2firs.or.ib.eq.i2last) then 929 : c write(6,*) 'le noeud',id,'detecte comme apres',ib 930 : xtau1=xicpr(ia) 931 : xtau2=xicpr(ib) 932 : stau = psc1 / v12lo 933 : xtau=stau*xtau2+(1.d0-stau)*xtau1 934 : else 935 : xtau=xicpr(ib) 936 : endif 937 : endif 938 : else 939 : pp2x=xx-(psc1*v12x) 940 : pp2y=yy-(psc1*v12y) 941 : pp2z=zz-(psc1*v12z) 942 : xd2 = abs(psc1) 943 : xd3 = sqrt((xxp-x11)**2+(yyp-y11)**2+(zzp-z11)**2) 944 : c if(ideux.ge.2) xtau=xicpr(ia) 945 : if(ideux.ge.2) then 946 : if(ia.eq.i2firs.or.ia.eq.i2last) then 947 : c write(6,*) 'le noeud',id,'detecte comme avant',ia 948 : xtau1=xicpr(ia) 949 : xtau2=xicpr(ib) 950 : stau = psc1 / v12lo 951 : xtau=stau*xtau2+(1.d0-stau)*xtau1 952 : else 953 : xtau=xicpr(ia) 954 : endif 955 : endif 956 : endif 957 : dis1=sqrt((xxp-pp2x)**2+(yyp-pp2y)**2+(zzp-pp2z)**2) 958 : dis1=sign(dis1,xdd1) 959 : endif 960 : 961 : endif 962 : * -fin de la distinction si noeud du tri3 ou pas 963 : 964 : * 1ere fois qu on traite ce point? 965 : if (isupfi(ipt).eq.0) then 966 : * isupfi(ipt)=1 967 : * on utilise isupfi pour definir si la valeur inscrite 968 : * provient du 1er ou dernier seg du front (=2) ou autre(=1)? 969 : isupfi(ipt)=ifis(iseg,1) 970 : igardf=igardf+1 971 : xdmin2(ipt)= xd2 972 : xdmin3(ipt)= xd3 973 : cbp2013 mpova1.vpocha(ipt,1)=dis1 974 : if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau 975 : else 976 : * if(abs(dis1).lt.abs(mpova1.vpocha(ipt,1)))then 977 : * on a trouvé un segment + proche ou = du noeud considéré 978 : * if(abs(dis1).le.abs(mpova1.vpocha(ipt,1)))then 979 : ** if(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))-1.D-8)then 980 : * if(xd2.le.xdmin2(ipt)) then 981 : if (xd3.lt.(abs(xdmin3(ipt))-1.D-8)) then 982 : cbp2013 mpova1.vpocha(ipt,1)=dis1 983 : xdmin2(ipt)=xd2 984 : xdmin3(ipt)=xd3 985 : isupfi(ipt)=ifis(iseg,1) 986 : if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau 987 : endif 988 : * rattrapage dans le cas d'erreur d'arrondi sur dis1 989 : * elseif(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))+1.D-8) 990 : * $ then 991 : * if(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))+1.D-8)then 992 : if (xd3.lt.(abs(xdmin3(ipt))+1.D-8))then 993 : * en cas d'égalité,on n'utilise la distance dis1 que si: 994 : * - le pt projetté appartient au segment (xd2=0) 995 : * ou - le point projeté est + proche du segment 996 : if (xd2.lt.xdmin2(ipt)) then 997 : cbp2013 mpova1.vpocha(ipt,1)=dis1 998 : xdmin2(ipt)=xd2 999 : xdmin3(ipt)=xd3 1000 : isupfi(ipt)=ifis(iseg,1) 1001 : if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau 1002 : endif 1003 : endif 1004 : endif 1005 : 1006 : 490 continue 1007 : * Fin de Boucle sur les noeuds =============================== 1008 : *======================================================================= 1009 : 1010 : 400 continue 1011 : * Fin de Boucle sur les elements de la fissure ========================= 1012 : *======================================================================= 1013 : 1014 : 1015 : *------Petite correction de Psi 1016 : if (ideux.ge.1) then 1017 : * Boucle sur les noeuds =============================== 1018 : do 500 ipt=1,npt 1019 : if (isup(ipt).eq.0) goto 500 1020 : c * pour le cas de fissure debouchante 1021 : c * dans le cas ou c'est un segment debouchant qui a servi au calcul de psi, 1022 : c * il faut prendre la distance xdd2 (et pas dis) qu'on retrouve avec : 1023 : c if ((isupfi(ipt)).eq.2) then 1024 : c dis1 = mpova1.vpocha(ipt,1) 1025 : c xdd1 = sqrt( (dis1**2) - (xdmin2(ipt)**2) ) 1026 : c mpova1.vpocha(ipt,1) = sign(xdd1,dis1) 1027 : c endif 1028 : c bp : essai du 16/03/2012 : psi = distance au front**2-phi**2 1029 : xdd1 = (xdmin3(ipt)**2) - (mpova2.vpocha(ipt,1)**2) 1030 : xdd1 = sqrt(max(xdd1,0.d0)) 1031 : if(ipt3.num(1,ipt).eq.idbug) then 1032 : write(6,*)'ipt,isup,psi_avant,phi,xdd1,= ',ipt,isup(ipt), 1033 : $ xdmin3(ipt),mpova1.vpocha(ipt,1),mpova2.vpocha(ipt,1),xdd1 1034 : endif 1035 : c if (isup(ipt).gt.0) then 1036 : if (isup(ipt).ge.15) then 1037 : *on est a l exterieur => psi >0 1038 : mpova1.vpocha(ipt,1) = xdd1 1039 : else 1040 : mpova1.vpocha(ipt,1) = -1.d0*xdd1 1041 : endif 1042 : 500 continue 1043 : endif 1044 : *------Fin de la Petite correction de Psi 1045 : 1046 : 1047 : *********************************************************************** 1048 : * FIN DU TRAVAIL 1049 : *********************************************************************** 1050 : if(iimpi.ge.1) then 1050 : call gibtem(xkt) 1051 : write(ioimp,*)'fin du travail | xkt=',xkt 1052 : endif 1053 : 1054 : * on inverse isup 1055 : * AVANT : isup(i) =0 veut dire que le point na pas été traité 1056 : * APRES : isup(i) =1 veut dire que le point peut etre enleve (=0 sinon) 1057 : iaju=0 1058 : do ia=1,isup(/1) 1059 : if (isup(ia).eq.0) then 1060 : iaju=iaju+1 1061 : isup(ia)=1 1062 : else 1063 : isup(ia)=0 1064 : endif 1065 : enddo 1066 : if (iaju+igard.ne.isup(/1)) then 1067 : write(ioimp,*) 'psiphi : tout va mal ' 1068 : $ ,iaju,'+',igard,'ne',isup(/1) 1068 : call erreur(26) 1069 : return 1070 : endif 1071 : ipt1=ipt3 1072 : ipt4=melfis 1073 : segdes ipt4 1074 : * iaju est le nombre de points a éliminer 1075 : * isup(i) =1 veut dire que le point peut etre enleve (=0 sinon) 1076 : 1077 : * faut-il ajuster les champs? 1078 : if (iaju.ne.0) then 1079 : * si oui on recree les bons chpoints 1080 : nbelem=ipt1.num(/2)- iaju 1081 : nbnn=1 1082 : nbref=0 1083 : nbsous=0 1084 : segini ipt2 1085 : ipt2.itypel=1 1086 : nc=1 1087 : n=nbelem 1088 : segini mpova5 1089 : *>>>bp&yt: ajout ligne suivante 1090 : if(ideux.ge.1) segini mpova4 1091 : if(ideux.ge.2) segini mpova6 1092 : iel=0 1093 : do i=1,ipt1.num(/2) 1094 : if (isup(i).eq.0) then 1095 : iel=iel+1 1096 : ipt2.num(1,iel)=ipt1.num(1,i) 1097 : mpova5.vpocha(iel,1)=mpova2.vpocha(i,1) 1098 : *>>>bp&yt: ajout ligne suivante 1099 : * on met le meme maillage pour les 2 level set (car toujours utilisées ensemble) 1100 : if(ideux.ge.1) mpova4.vpocha(iel,1)=mpova1.vpocha(i,1) 1101 : if(ideux.ge.2) mpova6.vpocha(iel,1)=mpova3.vpocha(i,1) 1102 : endif 1103 : enddo 1104 : if (ideux.ge.2) then 1105 : c segact msoup3*mod 1106 : msoup3.ipoval=mpova6 1107 : msoup3.igeoc=ipt2 1108 : segdes msoup3,mpova6 1109 : segsup,mpova3 1110 : endif 1111 : if (ideux.ge.1) then 1112 : c segact msoup1*mod 1113 : msoup1.ipoval=mpova4 1114 : msoup1.igeoc=ipt2 1115 : segdes msoup1,mpova4 1116 : segsup,mpova1 1117 : endif 1118 : c segact msoup2*mod 1119 : msoup2.ipoval=mpova5 1120 : msoup2.igeoc=ipt2 1121 : segdes msoup2,mpova5,ipt2 1122 : segsup,mpova2 1123 : segdes ipt1 1124 : else 1125 : if(ideux.ge.2) segdes,msoup3,mpova3 1126 : if(ideux.ge.1) segdes,msoup1,mpova1 1127 : segdes msoup2,mpova2,ipt1 1128 : endif 1129 : if (ideux.ge.1) then 1130 : segsup,icpr,sfis 1131 : if(ideux.ge.2) segsup,xicpr 1132 : endif 1133 : if (icle.eq.1) then 1134 : segsup isup,isupfi,xdis,xdmin,xdmin1,rrlim 1135 : if(ideux.ge.1) segsup xdmin2,xdmin3 1136 : endif 1137 : 1138 : RETURN 1139 : END 1140 : 1141 : 1142 : 1143 : 1144 :
© Cast3M 2003 - All rights reserved.
Disclaimer