C PSIP2D SOURCE CB215821 20/11/25 13:37:55 10792 c SUBROUTINE PSIP2D(ideux,melmai,melfis,ip1,ip2,xcrit, & msoup1,msoup2) *********************************************************************** c OPERATEUR : PSIP c c APPELE PAR PSIPHI DANS LE CAS 2D c c CREATION : BP 14/03/2012 c MODIF : BP 18/12/2013 : plus d'appel a zonag2, mais travail sur critere direct c *********************************************************************** IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) -INC PPARAM -INC CCOPTIO -INC SMCOORD -INC SMELEME -INC SMCHPOI -INC CCREEL segment ttrav integer ilig(inoe) endsegment segment isup(npt) c segment xdis(npt) c segment xdmin(npt) c segment xdmin1(npt) c segment rrlim(npt) c segment xdmin2(npt) c segment xdmin3(npt) c segment xrr2b(npt) c segment xpsc1b(npt) c segment isupfi(npt) c SEGMENT ISEG3 c INTEGER NIZO(NZO+1) c ENDSEGMENT c c SEGMENT ISEG4 c c INTEGER NUMZO(NZO) c c ENDSEGMENT c SEGMENT ISEG5 c INTEGER NNMEL(ILON),IDEJ(NZO) c ENDSEGMENT DATA EPSI/1.D-5/ c REAL*8 XZONE(3,3) c INTEGER NZONE(3) idim1=idim+1 *********************************************************************** * Initialisation, Recup, Activation *********************************************************************** idim1=idim+1 idbug = 0 if(msoup2.ne.0) mpova2=msoup2.ipoval if(msoup1.ne.0) mpova1=msoup1.ipoval c activation des maillages / desactivation c melmai : fait dans psiphi / ligne 1325 ou 1329 c melfis : fait ligne 56 / ligne 72 *********************************************************************** * Preparation du maillage de la fissure *********************************************************************** meleme=melfis segact meleme if (lisous(/1).ne.0) then call erreur(25) return endif if (itypel.ne.2.and.itypel.ne.3) then * seg2 et seg3 seulement! call erreur(16) return endif * on ordonne la ligne depuis la pointe de fissure call ligmai(meleme,ttrav,0) if(ierr.ne.0) return segact ttrav*mod * on regarde si ip1 est bien en dernier nb=ilig(/1) if (ideux.ge.1) then if (ilig(nb).eq.ip1) then * rien a faire else * il faut inverser le sens de parcours de la ligne nc=nb+1 nn=nb/2 do i=1,nn io=ilig(i) iu=ilig(nc-i) ilig(i)=iu ilig(nc-i)=io enddo endif endif * on verifie que ip2 est bien le premier point if (ip2.ne.0) then if (ilig(1).ne.ip2) then call erreur (21) endif endif nseg=nb-1 * sont actifs: melmai (=struct changée en poi1) *********************************************************************** * debut du travail suivant valeur de icle if (xcrit.eq.0.d0) then icle=0 xcrit=1.D10 else icle=1 * ajout d'un epsilon 1.D-5 de rattrapage xcrit = 1.00001D0 * xcrit endif if(iimpi.ge.1) then call gibtem(xkt) write(ioimp,*)'debut zonage si icle=1=',icle,xcrit,'| xkt=',xkt endif c ***** cas icle=1 *************************************************** c if (icle.eq.1) then c c c write(6,*) '--cas icle=1' c c melfis=meleme c melpoi=melmai c c *---- ZONAGE ------------------------------------------------------- c call zonag2(melfis,melpoi,xcrit,ISEG3,ISEG5,XZONE,NZONE) c * recup c ILON=NNMEL(/1) c NZO=IDEJ(/1) c Xtomi=XZONE(1,1) c Ytomi=XZONE(2,1) c c Ztomi=XZONE(3,1) c Xpr=XZONE(1,2) c Ypr=XZONE(2,2) c c Zpr=XZONE(3,2) c XZO=XZONE(1,3) c YZO=XZONE(2,3) c c ZZO=XZONE(3,3) c nXzo=NZONE(1) c nYzo=NZONE(2) c c nZzo=NZONE(3) c c c *********************************************************************** c * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT c *********************************************************************** c c * pour chaque segment de la fissure (melfis 2eme maillage lu) c * on calcule le champ phi et psi des noeuds autour de l'element c if(iimpi.ge.1) then c call gibtem(xkt) c write(ioimp,*)'debut travail | xkt=',xkt c endif c * c meleme=melfis c c segact,meleme c ipt3=melpoi c c segact,ipt3 c npt=ipt3.num(/2) c * attention on initialise isup avec les valeurs inverses,il faudra c * le retourner c segini isup,isupfi,xdis,xdmin,xdmin1,rrlim c if(ideux.ge.1) segini xdmin2,xdmin3 c if(ip2.ne.0) segini xrr2b,xpsc1b c igard=0 c igardf=0 c * Boucle sur les elements de la fissure ========================= c c do 400 iseg=1,nseg c c if(idebug.ge.1) c c $ write(6,*) '====== element fissure ',iseg,' / ',nseg,ifis(iseg) c c ia=ilig(iseg) c ib=ilig(iseg+1) c x11=xcoor((ia-1)*idim1+1) c y11=xcoor((ia-1)*idim1+2) c x22=xcoor((ib-1)*idim1+1) c y22=xcoor((ib-1)*idim1+2) c vx=x22-x11 c vy=y22-y11 c * aa= sqrt(vx*vx + vy*vy) + 1.d-20 c * bp: dans ce cas, mieux vaut faire une erreur c aa= sqrt(vx*vx + vy*vy) c if (aa.le.XPETIT) then c write(ioimp,*) 'ERREUR: segment de fissure de longueur nulle!' c call erreur(26) c return c endif c vx=vx/aa c vy=vy/aa c vxyz = max(abs(vx),abs(vy)) c if (vxyz.le.XPETIT) then c write(ioimp,*) 'segment de fissure de longueur nulle! ' c $ ,aa,vx,vy c call erreur(26) c return c endif c if (x11.gt.x22) then c xxmi=x22 c xxma=X11 c else c xxma=x22 c xxmi=x11 c endif c if (y11.gt.y22) then c yymi=y22 c yyma=y11 c else c yyma=y22 c yymi=y11 c endif c * calcul des zones xmi,ymi et xma,yma c NIZ1X=INT((Xxmi-XTOMI-XPR)/XZO) +1 c NIZ1Y=INT((yymi-YTOMI-YPR)/YZO)+1 c NIZ2X=INT((xxma-XTOMI+XPR)/XZO) +1 c NIZ2Y=INT((yyma-YTOMI+YPR)/YZO) +1 c * on ajoute une zone de chaque cotés c NIZ1X=max(1,niz1x-1) c niz1y=max(niz1y-1,1) c niz2x=min(nxzo,niz2x+1) c niz2y=min(nyzo,niz2y+1) c c write(6,*)'bornes x y z ',niz1x,niz2x,niz1y,niz2y,NIZ1Z,NIZ2Z c c ncalc = 0 c c * on explore les zones ------------------------------------------------ c do 332 m1=niz1y,niz2y c do 332 m2=niz1x,niz2x c indzo=m2+(m1-1)*nxzo c IDEB=NIZO(INDZO) c IFIN=NIZO(INDZO+1)-1 c if(ideb.gt.ifin) go to 332 c do kk=ideb,ifin c ipt=NNMEL(KK) c ia=ipt3.num(1,ipt) c xx = xcoor((ia-1)*idim1 + 1 ) c yy = xcoor((ia-1)*idim1 + 2 ) c * calcul de la distance au premier pt c rr2=(xx-x11)*(xx-x11)+(yy-y11)*(yy-y11) c * calcul des produits scals et du produit vectoiel c psca1= vx*(xx-x11)+vy*(yy-y11) c psca2= vx*(xx-x22)+vy*(yy-y22) c pvec=vx*(yy-y11)-vy*(xx-x11) c xsi=psca1*psca2 c xd= rr2 - psca1*psca1 c if ((iseg.eq.1).and.(ip2.ne.0)) then c * dans le cas de 2 pointes on sauve la distance a la 1ere pointe c xrr2b(ipt) = rr2 c xpsc1b(ipt) = psca1 c endif c if (xsi.le.0.d0) then c * on est dans le cas ou le point se projette sur le segment c if(xd.lt.0.D0) xd=0.D0 c if (iseg.eq.1) then c xdmin(ipt)= xd c xdmin1(ipt)=pvec c xdis(ipt) = sign((sqrt(xd)),pvec) c else c if (xd.le.xdmin(ipt))then c xdmin(ipt)= xd c xdmin1(ipt)=pvec c xdis(ipt) = sign((sqrt(xd)),pvec) c endif c endif c else c * le point ne se projette pas sur le segment c * est-il plus proche du premier ou du dernier pt du segment c if (psca1.ge.0.D0)then c rr2=(xx-x22)*(xx-x22)+(yy-y22)*(yy-y22) c xd= rr2 - psca2*psca2 c endif c if(xd.lt.0.D0) xd=0.D0 c if (iseg.eq.1) then c xdmin(ipt)= rr2 c xdmin1(ipt)=pvec c xdis(ipt) = sign((sqrt(xd)),pvec) c else c c if (rr2.le.xdmin(ipt)) then c if (rr2.lt.0.999999*xdmin(ipt)) then c xdmin(ipt)= rr2 c xdmin1(ipt)=pvec c if (iseg.eq.nseg) then c xdis(ipt) = sign((sqrt(xd)),pvec) c else c xdis(ipt) = sign((sqrt(rr2)),pvec) c endif c elseif (rr2.lt.1.000001*xdmin(ipt)) then c ss1=pvec*xdmin(ipt) c if (ss1.le.0.d0) then c * bp 20.06.2011 : egalite entre 2 concurrents mais signes opposés ! c if (abs(pvec).gt.abs(xdmin1(ipt))) then c xdmin(ipt)= rr2 c xdmin1(ipt)=pvec c if (iseg.eq.nseg) then c xdis(ipt) = sign((sqrt(xd)),pvec) c else c xdis(ipt) = sign((sqrt(rr2)),pvec) c endif c endif c endif c endif c endif c endif c if ( isup(ipt).eq.0) then c * c'est la premiere fois que l'on voit ce point c isup(ipt)=1 c igard=igard+1 c endif c * si c'est le dernier segment on stocke les valeurs phi c if (iseg.eq.nseg) then c igardf=igardf+1 c isupfi(ipt)=1 c mpova2.vpocha(ipt,1)=xdis(ipt) c * puis on calcule les valeurs de psi c * ...relatif a la 1ere pointe de fissure c if (ideux.ge.1) then c rr2 = (xx-x22)*(xx-x22)+(yy-y22)*(yy-y22) c if (psca2.ge.0.) then c mpova1.vpocha(ipt,1)=psca2 c else c yd = (rr2 - xdis(ipt)*xdis(ipt)) c if(yd.lt.0.D0) yd=0.D0 c mpova1.vpocha(ipt,1)= -1.*(yd**0.5) c endif c * ...relatif a la 2eme pointe de fissure c if (ip2.ne.0) then c if (xpsc1b(ipt).le.0.) then c psi2= -1.*xpsc1b(ipt) c else c psi2= (xrr2b(ipt) - xdis(ipt)*xdis(ipt)) c if(psi2.lt.0.D0) then c psi2= 0.D0 c else c psi2= -1.*(psi2**0.5) c endif c endif c * et on choisit le psi le plus proche c cbp2013 if(psi2.gt.(mpova1.vpocha(ipt,1))) c if(rr2.gt.xrr2b(ipt)) c & mpova1.vpocha(ipt,1) = psi2 c endif c endif c endif c enddo c 332 continue c * c c c if (idebug.ge.1) then c c ncalc2 = ncalc2 + ncalc c c if (mod(iseg,1000).eq.0.or.iseg.eq.nseg) then c c call gibtem(xkt) c c write(6,*) '---------------xkt=',ncalc2,xkt c c ncalc2=0 c c endif c c endif c c 400 continue c * Fin de Boucle sur les elements de la fissure ========================= c c c * on inverse isup c iaju=0 c do ia=1,isup(/1) c if (isup(ia).eq.0) then c iaju=iaju+1 c isup(ia)=1 c else c isup(ia)=0 c endif c enddo c if (iaju+igard.ne.isup(/1)) then c write(ioimp,*) 'psiphi : tout va mal ' c $ ,iaju,'+',igard,'ne',isup(/1) c call erreur(26) c return c endif c ipt1=ipt3 c c segsup iseg3,iseg4,iseg5 c segsup iseg3,iseg5 c c c ***** cas icle=0 *************************************************** c else c * recherche d'une boite autour de la fissure c if (ircrit.ne.0) then c xcrit= xcrit*2.1 c xmin=xgrand c xmax=-xgrand c ymin=xgrand c ymax=-xgrand c do i=1,nb c ia=ilig(i) c xmin= min( xmin,xcoor((ia-1)*idim1+1)) c ymin= min( ymin,xcoor((ia-1)*idim1+2)) c xmax= max( xmax,xcoor((ia-1)*idim1+1)) c ymax= max( ymax,xcoor((ia-1)*idim1+2)) c enddo c xmin = xmin-xcrit c xmax=xmax+xcrit c ymin=ymin-xcrit c ymax=ymax+xcrit c * write(6,*) ' xcrit,xmin,xmax,ymin,ymax' c * write(6,*) xcrit,xmin,xmax,ymin,ymax c endif * * il faut maintenant remplir les champs * * pour chaque pt de ipt1 on calcule le produit scalaire et vectoriel * par rapport à chaque segment de la ligne ipt1=melmai npt=ipt1.num(/2) segini isup igard=0 * * boucle sur les points de ipt1 **************************************** * do 100 ipt=1,ipt1.num(/2) * recup du ipt ieme point #ia : x ia= ipt1.num(1,ipt) xx = xcoor((ia-1)*idim1 + 1 ) yy = xcoor((ia-1)*idim1 + 2 ) c write(6,*) 'point',ipt,' #',ia,' x=',xx,' y=',yy c qq initialisations ivu = 0 rr22=XGRAND psca1b=XGRAND * boucle sur les segments ---------------------------------------- do 101 iseg=1,nseg * recup du segment 12 ia=ilig(iseg) ib=ilig(iseg+1) x11=xcoor((ia-1)*idim1+1) y11=xcoor((ia-1)*idim1+2) x22=xcoor((ib-1)*idim1+1) y22=xcoor((ib-1)*idim1+2) c write(6,*)'segment fissure',iseg,' 1:',x11,y11,' 2:',x22,y22 if (icle.eq.1) then xmin = min(x11,x22) - xcrit xmax = max(x11,x22) + xcrit ymin = min(y11,y22) - xcrit ymax = max(y11,y22) + xcrit endif *bp2013 on ne saute pas le segments du front (=1er et dernier segments) if(iseg.eq.1.or.iseg.eq.nseg) goto 111 *bp2013 on saute ce segment si point trop loin (critere) if (icle.eq.1) then if(xx.lt.xmin) goto 101 if(xx.gt.xmax) goto 101 if(yy.lt.ymin) goto 101 if(yy.gt.ymax) goto 101 endif 111 continue * vecteur tangent v=12/|12| vx=x22-x11 vy=y22-y11 aa= sqrt(vx*vx + vy*vy) aa=max(aa,1.d-20) vx=vx/aa vy=vy/aa * calcul de la distance au premier pt rr2=(xx-x11)*(xx-x11)+(yy-y11)*(yy-y11) * calcul des produits scals et du produit vectoiel psca1= vx*(xx-x11)+vy*(yy-y11) psca2= vx*(xx-x22)+vy*(yy-y22) pvec = vx*(yy-y11)-vy*(xx-x11) xsi=psca1*psca2 xd= rr2 - psca1*psca1 c write(6,*) 'xsi,xd,pvec,psca1,rr2=',xsi,xd,pvec,psca1,rr2 * Cas des pointes (=1er et dernier segments) if (iseg.eq.1.or.iseg.eq.nseg) then * on sauve la distance a la 2eme pointe (1er segment) if (iseg.eq.1) then rr22 = rr2 psca1b = psca1 endif * on re-teste ici le critere if (icle.eq.1) then if(xx.lt.xmin) goto 101 if(xx.gt.xmax) goto 101 if(yy.lt.ymin) goto 101 if(yy.gt.ymax) goto 101 endif endif * c'est la premiere fois que l'on voit ce point c write(6,*) 'xminmax=',xmin,xmax,' yminmax=',ymin,ymax if (isup(ipt).eq.0) then isup(ipt)=1 igard=igard+1 c write(6,*) 'isup=',isup(ipt),' --> igard=',igard endif ivu=ivu+1 * -on est dans le cas ou le point se projette sur le segment if (xsi.le.0.d0) then if(xd.lt.0.D0) xd=0.D0 if (iseg.eq.1) then dmin= xd dmin11=pvec dis = sign((sqrt(xd)),pvec) else if (xd.le.dmin)then dmin= xd dmin11=pvec dis = sign((sqrt(xd)),pvec) endif endif * -le point ne se projette pas sur le segment else * est-il plus proche du premier ou du dernier pt du segment ? if (psca1.ge.0.D0) then rr2=(xx-x22)*(xx-x22)+(yy-y22)*(yy-y22) xd= rr2 - psca2*psca2 endif if(xd.lt.0.D0) xd=0.D0 * --si 1ere fois qu on traite ce point if (ivu.eq.1) then dmin= rr2 dmin11=pvec dis = sign((sqrt(xd)),pvec) * --on a deja vu ce point else * ---si, avec ce segment, on est franchement + pres c if (rr2.le.dmin) then if (rr2.lt.0.999999*dmin) then dmin= rr2 dmin11=pvec if (iseg.eq.nseg) then dis = sign((sqrt(xd)),pvec) else dis = sign((sqrt(rr2)),pvec) endif * ---si on est a quasi-egalite avec un autre segment de fissure elseif (rr2.lt.1.000001*dmin) then ss1=pvec*dmin if (ss1.le.0.d0) then * bp 20.06.2011 : egalite entre 2 concurrents mais signes opposés ! if (abs(pvec).gt.abs(dmin11)) then dmin= rr2 dmin11=pvec if (iseg.eq.nseg) then dis = sign((sqrt(xd)),pvec) else dis = sign((sqrt(rr2)),pvec) endif endif endif endif * ---fin de la disctinction selon valeur de rr2 endif * --fin de la distinction vu ou pas endif * -fin de la distinction point se projette ou pas sur le segment c if(ia.eq.1799) c write(6,*) ipt,isup(ipt),psca1,xd,rr2,pvec,dis,dmin c enddo 101 continue * fin de boucle sur segment -------------------------------------- * -cas ou on stocke phi (ex-dernier segment) if (isup(ipt).ne.0) then c write(6,*) ipt,'stockage de phi =',dis * stockage de phi mpova2.vpocha(ipt,1)=dis * --puis on calcule les valeurs de psi ... if (ideux.ge.1) then * ...relatif a la 1ere pointe de fissure (dernier segment) rr2 = (xx-x22)*(xx-x22)+(yy-y22)*(yy-y22) if (psca2.ge.0.) then mpova1.vpocha(ipt,1)=psca2 else yd = (rr2 - dis*dis) if(yd.lt.0.D0) yd=0.D0 mpova1.vpocha(ipt,1)= -1.*(yd**0.5) endif * ---...relatif a la 2eme pointe de fissure (=le 1er segment) if (ip2.ne.0) then if (psca1b.le.0.) then psi2= -1.*psca1b else psi2= (rr22 - dis*dis) if (psi2.lt.0.D0) then psi2= 0.D0 else psi2= -1.*(psi2**0.5) endif endif * et on choisit le psi le plus proche cbp2013 if(psi2.gt.(mpova1.vpocha(ipt,1))) cbp2013 & mpova1.vpocha(ipt,1) = psi2 if(rr2.gt.rr22) mpova1.vpocha(ipt,1) = psi2 endif * ---fin du cas 2 pointes endif * --fin du cas ideux.ge.1 (calcul de psi) endif * -fin du cas ou on stocke phi (ex-dernier segment) 100 continue * fin de boucle sur point ****************************************** c endif c ***** fin de la distinction icle= 1 / 0 ******************************* *********************************************************************** * FIN DU TRAVAIL *********************************************************************** * faut-il ajuster les champs? iaju=0 do ia=1,isup(/1) * on inverse isup if (isup(ia).eq.0) then iaju=iaju+1 isup(ia)=1 else isup(ia)=0 endif enddo if (iaju+igard.ne.isup(/1)) then write(ioimp,*) 'psiphi : tout va mal ' $ ,iaju,'+',igard,'ne',isup(/1) call erreur(26) return endif * ajustement des champs if (iaju.ne.0) then nbelem=ipt1.num(/2)- iaju nbnn=1 nbref=0 nbsous=0 segini ipt2 ipt2.itypel=1 nc=1 n=nbelem segini mpova5 *>>>bp&yt: ajout ligne suivante if(ideux.ge.1) segini mpova4 if(ideux.ge.2) segini mpova6 iel=0 do i=1,ipt1.num(/2) if (isup(i).eq.0) then iel=iel+1 ipt2.num(1,iel)=ipt1.num(1,i) mpova5.vpocha(iel,1)=mpova2.vpocha(i,1) *>>>bp&yt: ajout ligne suivante * on met le meme maillage pour les 2 level set (car toujours utilisées ensemble) if(ideux.ge.1) mpova4.vpocha(iel,1)=mpova1.vpocha(i,1) if(ideux.ge.2) mpova6.vpocha(iel,1)=mpova3.vpocha(i,1) endif enddo if (ideux.ge.2) then c segact msoup3*mod msoup3.ipoval=mpova6 msoup3.igeoc=ipt2 segsup,mpova3 endif if (ideux.ge.1) then c segact msoup1*mod msoup1.ipoval=mpova4 msoup1.igeoc=ipt2 segsup,mpova1 endif c segact msoup2*mod msoup2.ipoval=mpova5 msoup2.igeoc=ipt2 segsup,mpova2 endif segsup ttrav segsup isup c if ((idim.eq.2).and.(icle.eq.1)) then c if (icle.eq.1) then c segsup xdis,xdmin,xdmin1,rrlim c if(ideux.ge.1) segsup xdmin2,xdmin3 c if(ip2.ne.0) segsup xrr2b,xpsc1b c endif end