1 : C PSIP2D SOURCE BP208322 13/12/19 21:15:05 7888 2 : c 3 : SUBROUTINE PSIP2D(ideux,melmai,melfis,ip1,ip2,xcrit, 4 : & msoup1,msoup2) 5 : *********************************************************************** 6 : c OPERATEUR : PSIP 7 : c 8 : c APPELE PAR PSIPHI DANS LE CAS 2D 9 : c 10 : c CREATION : BP 14/03/2012 11 : c MODIF : BP 18/12/2013 : plus d'appel a zonag2, mais travail sur critere direct 12 : c 13 : *********************************************************************** 14 : IMPLICIT REAL*8(A-H,O-Z) 15 : IMPLICIT INTEGER(I-N) 16 : -INC CCOPTIO 17 : -INC SMCOORD 18 : -INC SMELEME 19 : -INC SMCHPOI 20 : -INC CCREEL 21 : segment ttrav 22 : integer ilig(inoe) 23 : endsegment 24 : segment isup(npt) 25 : c segment xdis(npt) 26 : c segment xdmin(npt) 27 : c segment xdmin1(npt) 28 : c segment rrlim(npt) 29 : c segment xdmin2(npt) 30 : c segment xdmin3(npt) 31 : c segment xrr2b(npt) 32 : c segment xpsc1b(npt) 33 : c segment isupfi(npt) 34 : c SEGMENT ISEG3 35 : c INTEGER NIZO(NZO+1) 36 : c ENDSEGMENT 37 : c c SEGMENT ISEG4 38 : c c INTEGER NUMZO(NZO) 39 : c c ENDSEGMENT 40 : c SEGMENT ISEG5 41 : c INTEGER NNMEL(ILON),IDEJ(NZO) 42 : c ENDSEGMENT 43 : DATA EPSI/1.D-5/ 44 : c REAL*8 XZONE(3,3) 45 : c INTEGER NZONE(3) 46 : idim1=idim+1 47 : 48 : 49 : *********************************************************************** 50 : * Initialisation, Recup, Activation 51 : *********************************************************************** 52 : idim1=idim+1 53 : idbug = 0 54 : if(msoup2.ne.0) mpova2=msoup2.ipoval 55 : if(msoup1.ne.0) mpova1=msoup1.ipoval 56 : c activation des maillages / desactivation 57 : c melmai : fait dans psiphi / ligne 1325 ou 1329 58 : c melfis : fait ligne 56 / ligne 72 59 : 60 : *********************************************************************** 61 : * Preparation du maillage de la fissure 62 : *********************************************************************** 63 : meleme=melfis 64 : segact meleme 65 : 66 : if (lisous(/1).ne.0) then 66 : call erreur(25) 67 : return 68 : endif 69 : if (itypel.ne.2.and.itypel.ne.3) then 70 : * seg2 et seg3 seulement! 70 : call erreur(16) 71 : return 72 : endif 73 : * on ordonne la ligne depuis la pointe de fissure 73 : call ligmai(meleme,ttrav,0) 74 : if(ierr.ne.0) return 75 : segact ttrav*mod 76 : * on regarde si ip1 est bien en dernier 77 : segdes meleme 78 : nb=ilig(/1) 79 : if (ideux.ge.1) then 80 : if (ilig(nb).eq.ip1) then 81 : * rien a faire 82 : else 83 : * il faut inverser le sens de parcours de la ligne 84 : nc=nb+1 85 : nn=nb/2 86 : do i=1,nn 87 : io=ilig(i) 88 : iu=ilig(nc-i) 89 : ilig(i)=iu 90 : ilig(nc-i)=io 91 : enddo 92 : endif 93 : endif 94 : * on verifie que ip2 est bien le premier point 95 : if (ip2.ne.0) then 96 : if (ilig(1).ne.ip2) then 96 : call erreur (21) 97 : endif 98 : endif 99 : nseg=nb-1 100 : 101 : * sont actifs: melmai (=struct changée en poi1) 102 : 103 : 104 : *********************************************************************** 105 : * debut du travail suivant valeur de icle 106 : if (xcrit.eq.0.d0) then 107 : icle=0 108 : xcrit=1.D10 109 : else 110 : icle=1 111 : * ajout d'un epsilon 1.D-5 de rattrapage 112 : xcrit = 1.00001D0 * xcrit 113 : endif 114 : 115 : if(iimpi.ge.1) then 115 : call gibtem(xkt) 116 : write(ioimp,*)'debut zonage si icle=1=',icle,xcrit,'| xkt=',xkt 117 : endif 118 : 119 : c ***** cas icle=1 *************************************************** 120 : c if (icle.eq.1) then 121 : c 122 : c c write(6,*) '--cas icle=1' 123 : c 124 : c melfis=meleme 125 : c melpoi=melmai 126 : c 127 : c *---- ZONAGE ------------------------------------------------------- 127 : c call zonag2(melfis,melpoi,xcrit,ISEG3,ISEG5,XZONE,NZONE) 128 : c * recup 129 : c ILON=NNMEL(/1) 130 : c NZO=IDEJ(/1) 131 : c Xtomi=XZONE(1,1) 132 : c Ytomi=XZONE(2,1) 133 : c c Ztomi=XZONE(3,1) 134 : c Xpr=XZONE(1,2) 135 : c Ypr=XZONE(2,2) 136 : c c Zpr=XZONE(3,2) 137 : c XZO=XZONE(1,3) 138 : c YZO=XZONE(2,3) 139 : c c ZZO=XZONE(3,3) 140 : c nXzo=NZONE(1) 141 : c nYzo=NZONE(2) 142 : c c nZzo=NZONE(3) 143 : c 144 : c 145 : c *********************************************************************** 146 : c * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT 147 : c *********************************************************************** 148 : c 149 : c * pour chaque segment de la fissure (melfis 2eme maillage lu) 150 : c * on calcule le champ phi et psi des noeuds autour de l'element 151 : c if(iimpi.ge.1) then 151 : c call gibtem(xkt) 152 : c write(ioimp,*)'debut travail | xkt=',xkt 153 : c endif 154 : c * 155 : c meleme=melfis 156 : c c segact,meleme 157 : c ipt3=melpoi 158 : c c segact,ipt3 159 : c npt=ipt3.num(/2) 160 : c * attention on initialise isup avec les valeurs inverses,il faudra 161 : c * le retourner 162 : c segini isup,isupfi,xdis,xdmin,xdmin1,rrlim 163 : c if(ideux.ge.1) segini xdmin2,xdmin3 164 : c if(ip2.ne.0) segini xrr2b,xpsc1b 165 : c igard=0 166 : c igardf=0 167 : c * Boucle sur les elements de la fissure ========================= 168 : c 169 : c do 400 iseg=1,nseg 170 : c c if(idebug.ge.1) 171 : c c $ write(6,*) '====== element fissure ',iseg,' / ',nseg,ifis(iseg) 172 : c 173 : c ia=ilig(iseg) 174 : c ib=ilig(iseg+1) 175 : c x11=xcoor((ia-1)*idim1+1) 176 : c y11=xcoor((ia-1)*idim1+2) 177 : c x22=xcoor((ib-1)*idim1+1) 178 : c y22=xcoor((ib-1)*idim1+2) 179 : c vx=x22-x11 180 : c vy=y22-y11 181 : c * aa= sqrt(vx*vx + vy*vy) + 1.d-20 182 : c * bp: dans ce cas, mieux vaut faire une erreur 183 : c aa= sqrt(vx*vx + vy*vy) 184 : c if (aa.le.XPETIT) then 185 : c write(ioimp,*) 'ERREUR: segment de fissure de longueur nulle!' 185 : c call erreur(26) 186 : c return 187 : c endif 188 : c vx=vx/aa 189 : c vy=vy/aa 190 : c vxyz = max(abs(vx),abs(vy)) 191 : c if (vxyz.le.XPETIT) then 192 : c write(ioimp,*) 'segment de fissure de longueur nulle! ' 193 : c $ ,aa,vx,vy 193 : c call erreur(26) 194 : c return 195 : c endif 196 : c if (x11.gt.x22) then 197 : c xxmi=x22 198 : c xxma=X11 199 : c else 200 : c xxma=x22 201 : c xxmi=x11 202 : c endif 203 : c if (y11.gt.y22) then 204 : c yymi=y22 205 : c yyma=y11 206 : c else 207 : c yyma=y22 208 : c yymi=y11 209 : c endif 210 : c * calcul des zones xmi,ymi et xma,yma 211 : c NIZ1X=INT((Xxmi-XTOMI-XPR)/XZO) +1 212 : c NIZ1Y=INT((yymi-YTOMI-YPR)/YZO)+1 213 : c NIZ2X=INT((xxma-XTOMI+XPR)/XZO) +1 214 : c NIZ2Y=INT((yyma-YTOMI+YPR)/YZO) +1 215 : c * on ajoute une zone de chaque cotés 216 : c NIZ1X=max(1,niz1x-1) 217 : c niz1y=max(niz1y-1,1) 218 : c niz2x=min(nxzo,niz2x+1) 219 : c niz2y=min(nyzo,niz2y+1) 220 : c c write(6,*)'bornes x y z ',niz1x,niz2x,niz1y,niz2y,NIZ1Z,NIZ2Z 221 : c c ncalc = 0 222 : c 223 : c * on explore les zones ------------------------------------------------ 224 : c do 332 m1=niz1y,niz2y 225 : c do 332 m2=niz1x,niz2x 226 : c indzo=m2+(m1-1)*nxzo 227 : c IDEB=NIZO(INDZO) 228 : c IFIN=NIZO(INDZO+1)-1 229 : c if(ideb.gt.ifin) go to 332 230 : c do kk=ideb,ifin 231 : c ipt=NNMEL(KK) 232 : c ia=ipt3.num(1,ipt) 233 : c xx = xcoor((ia-1)*idim1 + 1 ) 234 : c yy = xcoor((ia-1)*idim1 + 2 ) 235 : c * calcul de la distance au premier pt 236 : c rr2=(xx-x11)*(xx-x11)+(yy-y11)*(yy-y11) 237 : c * calcul des produits scals et du produit vectoiel 238 : c psca1= vx*(xx-x11)+vy*(yy-y11) 239 : c psca2= vx*(xx-x22)+vy*(yy-y22) 240 : c pvec=vx*(yy-y11)-vy*(xx-x11) 241 : c xsi=psca1*psca2 242 : c xd= rr2 - psca1*psca1 243 : c if ((iseg.eq.1).and.(ip2.ne.0)) then 244 : c * dans le cas de 2 pointes on sauve la distance a la 1ere pointe 245 : c xrr2b(ipt) = rr2 246 : c xpsc1b(ipt) = psca1 247 : c endif 248 : c if (xsi.le.0.d0) then 249 : c * on est dans le cas ou le point se projette sur le segment 250 : c if(xd.lt.0.D0) xd=0.D0 251 : c if (iseg.eq.1) then 252 : c xdmin(ipt)= xd 253 : c xdmin1(ipt)=pvec 254 : c xdis(ipt) = sign((sqrt(xd)),pvec) 255 : c else 256 : c if (xd.le.xdmin(ipt))then 257 : c xdmin(ipt)= xd 258 : c xdmin1(ipt)=pvec 259 : c xdis(ipt) = sign((sqrt(xd)),pvec) 260 : c endif 261 : c endif 262 : c else 263 : c * le point ne se projette pas sur le segment 264 : c * est-il plus proche du premier ou du dernier pt du segment 265 : c if (psca1.ge.0.D0)then 266 : c rr2=(xx-x22)*(xx-x22)+(yy-y22)*(yy-y22) 267 : c xd= rr2 - psca2*psca2 268 : c endif 269 : c if(xd.lt.0.D0) xd=0.D0 270 : c if (iseg.eq.1) then 271 : c xdmin(ipt)= rr2 272 : c xdmin1(ipt)=pvec 273 : c xdis(ipt) = sign((sqrt(xd)),pvec) 274 : c else 275 : c c if (rr2.le.xdmin(ipt)) then 276 : c if (rr2.lt.0.999999*xdmin(ipt)) then 277 : c xdmin(ipt)= rr2 278 : c xdmin1(ipt)=pvec 279 : c if (iseg.eq.nseg) then 280 : c xdis(ipt) = sign((sqrt(xd)),pvec) 281 : c else 282 : c xdis(ipt) = sign((sqrt(rr2)),pvec) 283 : c endif 284 : c elseif (rr2.lt.1.000001*xdmin(ipt)) then 285 : c ss1=pvec*xdmin(ipt) 286 : c if (ss1.le.0.d0) then 287 : c * bp 20.06.2011 : egalite entre 2 concurrents mais signes opposés ! 288 : c if (abs(pvec).gt.abs(xdmin1(ipt))) then 289 : c xdmin(ipt)= rr2 290 : c xdmin1(ipt)=pvec 291 : c if (iseg.eq.nseg) then 292 : c xdis(ipt) = sign((sqrt(xd)),pvec) 293 : c else 294 : c xdis(ipt) = sign((sqrt(rr2)),pvec) 295 : c endif 296 : c endif 297 : c endif 298 : c endif 299 : c endif 300 : c endif 301 : c if ( isup(ipt).eq.0) then 302 : c * c'est la premiere fois que l'on voit ce point 303 : c isup(ipt)=1 304 : c igard=igard+1 305 : c endif 306 : c * si c'est le dernier segment on stocke les valeurs phi 307 : c if (iseg.eq.nseg) then 308 : c igardf=igardf+1 309 : c isupfi(ipt)=1 310 : c mpova2.vpocha(ipt,1)=xdis(ipt) 311 : c * puis on calcule les valeurs de psi 312 : c * ...relatif a la 1ere pointe de fissure 313 : c if (ideux.ge.1) then 314 : c rr2 = (xx-x22)*(xx-x22)+(yy-y22)*(yy-y22) 315 : c if (psca2.ge.0.) then 316 : c mpova1.vpocha(ipt,1)=psca2 317 : c else 318 : c yd = (rr2 - xdis(ipt)*xdis(ipt)) 319 : c if(yd.lt.0.D0) yd=0.D0 320 : c mpova1.vpocha(ipt,1)= -1.*(yd**0.5) 321 : c endif 322 : c * ...relatif a la 2eme pointe de fissure 323 : c if (ip2.ne.0) then 324 : c if (xpsc1b(ipt).le.0.) then 325 : c psi2= -1.*xpsc1b(ipt) 326 : c else 327 : c psi2= (xrr2b(ipt) - xdis(ipt)*xdis(ipt)) 328 : c if(psi2.lt.0.D0) then 329 : c psi2= 0.D0 330 : c else 331 : c psi2= -1.*(psi2**0.5) 332 : c endif 333 : c endif 334 : c * et on choisit le psi le plus proche 335 : c cbp2013 if(psi2.gt.(mpova1.vpocha(ipt,1))) 336 : c if(rr2.gt.xrr2b(ipt)) 337 : c & mpova1.vpocha(ipt,1) = psi2 338 : c endif 339 : c endif 340 : c endif 341 : c enddo 342 : c 332 continue 343 : c * 344 : c 345 : c c if (idebug.ge.1) then 346 : c c ncalc2 = ncalc2 + ncalc 347 : c c if (mod(iseg,1000).eq.0.or.iseg.eq.nseg) then 347 : c c call gibtem(xkt) 348 : c c write(6,*) '---------------xkt=',ncalc2,xkt 349 : c c ncalc2=0 350 : c c endif 351 : c c endif 352 : c 353 : c 400 continue 354 : c * Fin de Boucle sur les elements de la fissure ========================= 355 : c 356 : c 357 : c * on inverse isup 358 : c iaju=0 359 : c do ia=1,isup(/1) 360 : c if (isup(ia).eq.0) then 361 : c iaju=iaju+1 362 : c isup(ia)=1 363 : c else 364 : c isup(ia)=0 365 : c endif 366 : c enddo 367 : c if (iaju+igard.ne.isup(/1)) then 368 : c write(ioimp,*) 'psiphi : tout va mal ' 369 : c $ ,iaju,'+',igard,'ne',isup(/1) 369 : c call erreur(26) 370 : c return 371 : c endif 372 : c ipt1=ipt3 373 : c c segsup iseg3,iseg4,iseg5 374 : c segsup iseg3,iseg5 375 : c c ipt4=melmai 376 : c c segdes ipt4 377 : c ipt4=melfis 378 : c segdes ipt4 379 : c 380 : c 381 : c ***** cas icle=0 *************************************************** 382 : c else 383 : 384 : c * recherche d'une boite autour de la fissure 385 : c if (ircrit.ne.0) then 386 : c xcrit= xcrit*2.1 387 : c xmin=xgrand 388 : c xmax=-xgrand 389 : c ymin=xgrand 390 : c ymax=-xgrand 391 : c do i=1,nb 392 : c ia=ilig(i) 393 : c xmin= min( xmin,xcoor((ia-1)*idim1+1)) 394 : c ymin= min( ymin,xcoor((ia-1)*idim1+2)) 395 : c xmax= max( xmax,xcoor((ia-1)*idim1+1)) 396 : c ymax= max( ymax,xcoor((ia-1)*idim1+2)) 397 : c enddo 398 : c xmin = xmin-xcrit 399 : c xmax=xmax+xcrit 400 : c ymin=ymin-xcrit 401 : c ymax=ymax+xcrit 402 : c * write(6,*) ' xcrit,xmin,xmax,ymin,ymax' 403 : c * write(6,*) xcrit,xmin,xmax,ymin,ymax 404 : c endif 405 : * 406 : * il faut maintenant remplir les champs 407 : * 408 : * pour chaque pt de ipt1 on calcule le produit scalaire et vectoriel 409 : * par rapport à chaque segment de la ligne 410 : ipt1=melmai 411 : npt=ipt1.num(/2) 412 : segini isup 413 : igard=0 414 : * 415 : * boucle sur les points de ipt1 **************************************** 416 : * 417 : do 100 ipt=1,ipt1.num(/2) 418 : 419 : * recup du ipt ieme point #ia : x 420 : ia= ipt1.num(1,ipt) 421 : xx = xcoor((ia-1)*idim1 + 1 ) 422 : yy = xcoor((ia-1)*idim1 + 2 ) 423 : c write(6,*) 'point',ipt,' #',ia,' x=',xx,' y=',yy 424 : 425 : c qq initialisations 426 : ivu = 0 427 : rr22=XGRAND 428 : psca1b=XGRAND 429 : 430 : * boucle sur les segments ---------------------------------------- 431 : do 101 iseg=1,nseg 432 : 433 : * recup du segment 12 434 : ia=ilig(iseg) 435 : ib=ilig(iseg+1) 436 : x11=xcoor((ia-1)*idim1+1) 437 : y11=xcoor((ia-1)*idim1+2) 438 : x22=xcoor((ib-1)*idim1+1) 439 : y22=xcoor((ib-1)*idim1+2) 440 : c write(6,*)'segment fissure',iseg,' 1:',x11,y11,' 2:',x22,y22 441 : if (icle.eq.1) then 442 : xmin = min(x11,x22) - xcrit 443 : xmax = max(x11,x22) + xcrit 444 : ymin = min(y11,y22) - xcrit 445 : ymax = max(y11,y22) + xcrit 446 : endif 447 : 448 : *bp2013 on ne saute pas le segments du front (=1er et dernier segments) 449 : if(iseg.eq.1.or.iseg.eq.nseg) goto 111 450 : *bp2013 on saute ce segment si point trop loin (critere) 451 : if (icle.eq.1) then 452 : if(xx.lt.xmin) goto 101 453 : if(xx.gt.xmax) goto 101 454 : if(yy.lt.ymin) goto 101 455 : if(yy.gt.ymax) goto 101 456 : endif 457 : 458 : 111 continue 459 : * vecteur tangent v=12/|12| 460 : vx=x22-x11 461 : vy=y22-y11 462 : aa= sqrt(vx*vx + vy*vy) 463 : aa=max(aa,1.d-20) 464 : vx=vx/aa 465 : vy=vy/aa 466 : * calcul de la distance au premier pt 467 : rr2=(xx-x11)*(xx-x11)+(yy-y11)*(yy-y11) 468 : * calcul des produits scals et du produit vectoiel 469 : psca1= vx*(xx-x11)+vy*(yy-y11) 470 : psca2= vx*(xx-x22)+vy*(yy-y22) 471 : pvec = vx*(yy-y11)-vy*(xx-x11) 472 : xsi=psca1*psca2 473 : xd= rr2 - psca1*psca1 474 : c write(6,*) 'xsi,xd,pvec,psca1,rr2=',xsi,xd,pvec,psca1,rr2 475 : 476 : * Cas des pointes (=1er et dernier segments) 477 : if (iseg.eq.1.or.iseg.eq.nseg) then 478 : * on sauve la distance a la 2eme pointe (1er segment) 479 : if (iseg.eq.1) then 480 : rr22 = rr2 481 : psca1b = psca1 482 : endif 483 : * on re-teste ici le critere 484 : if (icle.eq.1) then 485 : if(xx.lt.xmin) goto 101 486 : if(xx.gt.xmax) goto 101 487 : if(yy.lt.ymin) goto 101 488 : if(yy.gt.ymax) goto 101 489 : endif 490 : endif 491 : 492 : * c'est la premiere fois que l'on voit ce point 493 : c write(6,*) 'xminmax=',xmin,xmax,' yminmax=',ymin,ymax 494 : if (isup(ipt).eq.0) then 495 : isup(ipt)=1 496 : igard=igard+1 497 : c write(6,*) 'isup=',isup(ipt),' --> igard=',igard 498 : endif 499 : ivu=ivu+1 500 : 501 : * -on est dans le cas ou le point se projette sur le segment 502 : if (xsi.le.0.d0) then 503 : if(xd.lt.0.D0) xd=0.D0 504 : if (iseg.eq.1) then 505 : dmin= xd 506 : dmin11=pvec 507 : dis = sign((sqrt(xd)),pvec) 508 : else 509 : if (xd.le.dmin)then 510 : dmin= xd 511 : dmin11=pvec 512 : dis = sign((sqrt(xd)),pvec) 513 : endif 514 : endif 515 : * -le point ne se projette pas sur le segment 516 : else 517 : * est-il plus proche du premier ou du dernier pt du segment ? 518 : if (psca1.ge.0.D0) then 519 : rr2=(xx-x22)*(xx-x22)+(yy-y22)*(yy-y22) 520 : xd= rr2 - psca2*psca2 521 : endif 522 : if(xd.lt.0.D0) xd=0.D0 523 : * --si 1ere fois qu on traite ce point 524 : if (ivu.eq.1) then 525 : dmin= rr2 526 : dmin11=pvec 527 : dis = sign((sqrt(xd)),pvec) 528 : * --on a deja vu ce point 529 : else 530 : * ---si, avec ce segment, on est franchement + pres 531 : c if (rr2.le.dmin) then 532 : if (rr2.lt.0.999999*dmin) then 533 : dmin= rr2 534 : dmin11=pvec 535 : if (iseg.eq.nseg) then 536 : dis = sign((sqrt(xd)),pvec) 537 : else 538 : dis = sign((sqrt(rr2)),pvec) 539 : endif 540 : * ---si on est a quasi-egalite avec un autre segment de fissure 541 : elseif (rr2.lt.1.000001*dmin) then 542 : ss1=pvec*dmin 543 : if (ss1.le.0.d0) then 544 : * bp 20.06.2011 : egalite entre 2 concurrents mais signes opposés ! 545 : if (abs(pvec).gt.abs(dmin11)) then 546 : dmin= rr2 547 : dmin11=pvec 548 : if (iseg.eq.nseg) then 549 : dis = sign((sqrt(xd)),pvec) 550 : else 551 : dis = sign((sqrt(rr2)),pvec) 552 : endif 553 : endif 554 : endif 555 : endif 556 : * ---fin de la disctinction selon valeur de rr2 557 : endif 558 : * --fin de la distinction vu ou pas 559 : endif 560 : * -fin de la distinction point se projette ou pas sur le segment 561 : c if(ia.eq.1799) 562 : c write(6,*) ipt,isup(ipt),psca1,xd,rr2,pvec,dis,dmin 563 : 564 : c enddo 565 : 101 continue 566 : * fin de boucle sur segment -------------------------------------- 567 : 568 : * -cas ou on stocke phi (ex-dernier segment) 569 : if (isup(ipt).ne.0) then 570 : c write(6,*) ipt,'stockage de phi =',dis 571 : * stockage de phi 572 : mpova2.vpocha(ipt,1)=dis 573 : * --puis on calcule les valeurs de psi ... 574 : if (ideux.ge.1) then 575 : * ...relatif a la 1ere pointe de fissure (dernier segment) 576 : rr2 = (xx-x22)*(xx-x22)+(yy-y22)*(yy-y22) 577 : if (psca2.ge.0.) then 578 : mpova1.vpocha(ipt,1)=psca2 579 : else 580 : yd = (rr2 - dis*dis) 581 : if(yd.lt.0.D0) yd=0.D0 582 : mpova1.vpocha(ipt,1)= -1.*(yd**0.5) 583 : endif 584 : * ---...relatif a la 2eme pointe de fissure (=le 1er segment) 585 : if (ip2.ne.0) then 586 : if (psca1b.le.0.) then 587 : psi2= -1.*psca1b 588 : else 589 : psi2= (rr22 - dis*dis) 590 : if (psi2.lt.0.D0) then 591 : psi2= 0.D0 592 : else 593 : psi2= -1.*(psi2**0.5) 594 : endif 595 : endif 596 : * et on choisit le psi le plus proche 597 : cbp2013 if(psi2.gt.(mpova1.vpocha(ipt,1))) 598 : cbp2013 & mpova1.vpocha(ipt,1) = psi2 599 : if(rr2.gt.rr22) mpova1.vpocha(ipt,1) = psi2 600 : endif 601 : * ---fin du cas 2 pointes 602 : endif 603 : * --fin du cas ideux.ge.1 (calcul de psi) 604 : endif 605 : * -fin du cas ou on stocke phi (ex-dernier segment) 606 : 607 : 100 continue 608 : * fin de boucle sur point ****************************************** 609 : 610 : c endif 611 : c ***** fin de la distinction icle= 1 / 0 ******************************* 612 : 613 : 614 : *********************************************************************** 615 : * FIN DU TRAVAIL 616 : *********************************************************************** 617 : 618 : 619 : * faut-il ajuster les champs? 620 : iaju=0 621 : do ia=1,isup(/1) 622 : * on inverse isup 623 : if (isup(ia).eq.0) then 624 : iaju=iaju+1 625 : isup(ia)=1 626 : else 627 : isup(ia)=0 628 : endif 629 : enddo 630 : if (iaju+igard.ne.isup(/1)) then 631 : write(ioimp,*) 'psiphi : tout va mal ' 632 : $ ,iaju,'+',igard,'ne',isup(/1) 632 : call erreur(26) 633 : return 634 : endif 635 : 636 : * ajustement des champs 637 : if (iaju.ne.0) then 638 : nbelem=ipt1.num(/2)- iaju 639 : nbnn=1 640 : nbref=0 641 : nbsous=0 642 : segini ipt2 643 : ipt2.itypel=1 644 : nc=1 645 : n=nbelem 646 : segini mpova5 647 : *>>>bp&yt: ajout ligne suivante 648 : if(ideux.ge.1) segini mpova4 649 : if(ideux.ge.2) segini mpova6 650 : iel=0 651 : do i=1,ipt1.num(/2) 652 : if (isup(i).eq.0) then 653 : iel=iel+1 654 : ipt2.num(1,iel)=ipt1.num(1,i) 655 : mpova5.vpocha(iel,1)=mpova2.vpocha(i,1) 656 : *>>>bp&yt: ajout ligne suivante 657 : * on met le meme maillage pour les 2 level set (car toujours utilisées ensemble) 658 : if(ideux.ge.1) mpova4.vpocha(iel,1)=mpova1.vpocha(i,1) 659 : if(ideux.ge.2) mpova6.vpocha(iel,1)=mpova3.vpocha(i,1) 660 : endif 661 : enddo 662 : if (ideux.ge.2) then 663 : c segact msoup3*mod 664 : msoup3.ipoval=mpova6 665 : msoup3.igeoc=ipt2 666 : segdes msoup3,mpova6 667 : segsup,mpova3 668 : endif 669 : if (ideux.ge.1) then 670 : c segact msoup1*mod 671 : msoup1.ipoval=mpova4 672 : msoup1.igeoc=ipt2 673 : segdes msoup1,mpova4 674 : segsup,mpova1 675 : endif 676 : c segact msoup2*mod 677 : msoup2.ipoval=mpova5 678 : msoup2.igeoc=ipt2 679 : segdes msoup2,mpova5,ipt2 680 : segsup,mpova2 681 : segdes ipt1 682 : else 683 : if(ideux.ge.2) segdes,msoup3,mpova3 684 : if(ideux.ge.1) segdes,msoup1,mpova1 685 : segdes msoup2,mpova2,ipt1 686 : endif 687 : segsup ttrav 688 : segsup isup 689 : c if ((idim.eq.2).and.(icle.eq.1)) then 690 : c if (icle.eq.1) then 691 : c segsup xdis,xdmin,xdmin1,rrlim 692 : c if(ideux.ge.1) segsup xdmin2,xdmin3 693 : c if(ip2.ne.0) segsup xrr2b,xpsc1b 694 : c endif 695 : 696 : 697 : return 698 : end 699 : 700 : 701 : 702 : 703 : 704 : 705 : 706 : 707 : 708 : 709 :
© Cast3M 2003 - All rights reserved.
Disclaimer