C PSIP3D SOURCE BP208322 14/01/29 21:15:14 7921 c SUBROUTINE PSIP3D(ideux,melmai,melfis,melfro,xcrit, & msoup1,msoup2,msoup3) *********************************************************************** c OPERATEUR : PSIP c APPELE PAR PSIPHI DANS LE CAS 3D c c CREATION : BP, 14/03/2012 c MODIF : BP, 19/12/2013 plus d'appel au zonage c *********************************************************************** IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) -INC CCOPTIO -INC SMCOORD -INC SMELEME -INC SMCHPOI -INC CCREEL segment icpr(xcoor(/1)/(idim+1)) segment xicpr(xcoor(/1)/(idim+1)) segment isup(npt) segment xdis(npt) segment xdmin(npt) segment xdmin1(npt) segment rrlim(npt) segment xdmin2(npt) segment xdmin3(npt) segment isupfi(npt) segment sfis integer ifis(nbelfi,2) endsegment DATA EPSI/1.D-5/ logical phicor,phi12,phi23,phi31 integer ooperm(7,2) DATA ooperm/ 2, 4, 6, 1, 3, 5, 7, $ 4, 1, 5, 2, 6, 3, 7/ *********************************************************************** * Initialisation, Recup, Activation *********************************************************************** if(iimpi.ge.1) then call gibtem(xkt) write(ioimp,*)'===== ENTREE DANS PSIP3D ===== | xkt=',xkt endif idim1=idim+1 idbug = 0 c des chpoints (les msoup sont deja actifs) if(msoup2.gt.0) mpova2=msoup2.ipoval if(msoup1.gt.0) mpova1=msoup1.ipoval if(msoup3.gt.0) mpova3=msoup3.ipoval c activation des maillages / desactivation c melmai : fait dans psiphi / ligne 1325 ou 1329 c melfis : fait ligne 64 / ligne 1274 c melfro : fait (si ideux>=1) ligne 78 / ligne 356 ipt3=melmai c segact,ipt3 npt=ipt3.num(/2) c variables parfois testées imil5=0 *********************************************************************** * Preparation du maillage de la fissure *********************************************************************** meleme=melfis segact meleme * transfo du front en TRI3 if(itypel.ne.4) call change(melfis,4) meleme=melfis segact meleme c rem: inutile car change laisse melfis actif nseg=num(/2) *------ Cas ou l'on calcule phi ET psi -------------------------------- if (ideux.ge.1) then * on rearrange le melfis de facon à prevenir dans ifis(iel,1) quels * element touchent le front de fissure. Ceux ci sont ordonnées de facon * que le front de fissure soit du noeud 1 vers le noeud 2 de l'element ipt2=melfro segact ipt2 segini icpr nbelem2 = ipt2.num(/2) nbnode2 = ipt2.num(/1) * petite numerotation locale * icpr(noeud)>0 si le noeud appartient au front de fissure * =numero du 1er element de melfro ou on a vu ce noeud * + on en profite pour calculer l'abscisse curviligne du front si ideux.eq.2 * icpr(noeud)<0 si le noeud appartient au reste du contour de la fissure * =-numero du 1er element de melcon ou on a vu ce noeud if (ideux.ge.2) then segini,xicpr xtau=0.d0 ic=ipt2.num(1,1) x11= xcoor((ic-1)*idim1 +1) y11= xcoor((ic-1)*idim1 +2) z11= xcoor((ic-1)*idim1 +3) do ia=1,nbelem2 do ib=1,nbnode2 ic=ipt2.num(ib,ia) if (icpr(ic).eq.0) then icpr(ic)=ia x22= xcoor((ic-1)*idim1 +1) y22= xcoor((ic-1)*idim1 +2) z22= xcoor((ic-1)*idim1 +3) dxtau= sqrt( (x22-x11)*(x22-x11)+(y22-y11)*(y22-y11) $ +(z22-z11)*(z22-z11) ) xtau=xtau+dxtau xicpr(ic)=xtau x11=x22 y11=y22 z11=z22 c if(idebug.ge.1) c $ write(6,*) 'ia,ib,ic,xicpr(ic)=',ia,ib,ic,xicpr(ic) endif enddo enddo else do ia=1,nbelem2 do ib=1,nbnode2 ic=ipt2.num(ib,ia) if (icpr(ic).eq.0) then icpr(ic)=ia endif enddo enddo endif * bp (23/02/2012) : calcul du contour et fin de remplissage de icpr CALL ECROBJ('MAILLAGE',melfis) CALL PRCONT CALL LIROBJ('MAILLAGE',melcon,1,IRETOU) IF (IERR.NE.0) RETURN meleme = melfis IPT6 = melcon SEGACT,meleme,IPT6 nbelCON = IPT6.num(/2) nbnoCON = IPT6.num(/1) do ia=1,nbelCON do ib=1,nbnoCON ic=IPT6.num(ib,ia) if(icpr(ic).eq.0) icpr(ic)=-ia enddo enddo * nombre de segments appartenant au contour sans le front nbelem6 = nbelCON - nbelem2 * remplissage effectif de ifis contenu dans le segment sfis i2firs = ipt2.num(1,1) i2last = ipt2.num(nbnode2,nbelem2) ivui2f = 0 ivui2l = 0 c write(6,*) 'i2firs,i2last=',i2firs,i2last nbelfi=nseg nbelem=nseg nbnn=3 nbsous=0 nbref=0 c segini,ipt5=meleme * bp 20/12/2011 : on propose de remplir ipt5 avec : c + des elements du front (et du contour) d abord, c + des elements avec 1 point appartenant au front (et au contour) ensuite, c + les autres elements enfin. c bp 06/03/2012 : on propose de remplir ifis avec : c ifis(iseg,1) =2 si le segment [1-2] du iseg^ieme triangle appartient c au front de fissure, =3 si en plus c'est l extrémité du front, c =1 si l un des segments appartient au contour c =-1 si l'un des noeuds appartient au front c ifis(iseg,2) =+1+2+4 si les segments [1-2] [2-3] [3-1] du iseg^ieme triangle appartiennent au contour segini,ipt5 ideb5 = 0 imil5 = nbelem2 ifin5 = nseg+1 segini sfis *--------boucle sur les elements de la fissure ---------------------- do iel=1,nseg * io est incrémenté de +1,+2,+4 si node 1,2,3 appartient au front * ioo est incrémenté de +1,+2,+4 si node 1,2,3 appartient au contour io=0 ioo=0 * ---boucle sur les noeuds des elements de la fissure ---------- do ib=1,num(/1) ia=ib if(ib.eq.3) ia = 4 ic = num(ib,iel) if(icpr(ic).gt.0) io=io+ia if(icpr(ic).lt.0) ioo=ioo+ia * cas particulier d un segment avec 1 point appartenant au contour et un autre if(ic.eq.i2firs.or.ic.eq.i2last) ioo=ioo+ia enddo * ---fin de boucle sur les noeuds des elements de la fissure --- * remplissage de ipt5(i5,...) et ifis(i5,1) i5=0 iperm=0 * ---cas ou 1 segment appartient au front if (io.eq.3.or.io.gt.4) then ideb5=ideb5+1 i5=ideb5 * (eventuelle) permutation des noeuds pour avoir segment 1-2 = front if (io.eq.3) then ipt5.num(1,i5)=num(1,iel) ipt5.num(2,i5)=num(2,iel) ipt5.num(3,i5)=num(3,iel) elseif(io.eq.5) then ipt5.num(1,i5)=num(3,iel) ipt5.num(2,i5)=num(1,iel) ipt5.num(3,i5)=num(2,iel) iperm=1 elseif(io.eq.6) then ipt5.num(1,i5)=num(2,iel) ipt5.num(2,i5)=num(3,iel) ipt5.num(3,i5)=num(1,iel) iperm=2 elseif(io.eq.7) then write(ioimp,*) 'ERREUR: ' $ ,'Plusieurs cotes d un meme triangle du maillage de ' $ ,'la fissure converti en TRI3 appartiennent au front!' return endif * ifis( ,1)=3 si c'est une extremité if ((ipt5.num(1,i5)).eq.i2firs) then ifis(i5,1)=3 ivui2f = ivui2f+1 elseif((ipt5.num(2,i5)).eq.i2last) then ifis(i5,1)=3 ivui2l = ivui2l+1 else * ifis( ,1)=2 si c'est un tri3 du front avec des voisins ifis(i5,1)=2 endif * ---cas ou 1 (ou 2) segment(s) appartient au contour elseif(ioo.eq.3.or.ioo.gt.4) then imil5=imil5+1 i5=imil5 * pas de permutation des noeuds ipt5.num(1,i5)=num(1,iel) ipt5.num(2,i5)=num(2,iel) ipt5.num(3,i5)=num(3,iel) * ifis( ,1)=1 si un segment appartient au contour ifis(i5,1)=1 * ---cas ou seulement 1 point appartient au front elseif(io.eq.1.or.io.eq.2.or.io.eq.4) then ifin5=ifin5-1 i5=ifin5 * pas de permutation des noeuds ipt5.num(1,i5)=num(1,iel) ipt5.num(2,i5)=num(2,iel) ipt5.num(3,i5)=num(3,iel) * ifis( ,1)=-1 si 1 seul point du tri3 appartient au front ifis(i5,1)=-1 * ---autres cas else ifin5=ifin5-1 i5=ifin5 * pas de permutation des noeuds ipt5.num(1,i5)=num(1,iel) ipt5.num(2,i5)=num(2,iel) ipt5.num(3,i5)=num(3,iel) endif * remplissage de ifis(i5,2) * ---uniquement si 1 (ou 2) segment(s) appartient au contour if(ioo.eq.3.or.ioo.gt.4) then * attention a l eventuelle permutation effectuee if(iperm.ne.0) ioo=ooperm(ioo,iperm) * -ioo=3 => noeuds 1-2 = segment 1 appartient au contour if (ioo.eq.3) then ifis(i5,2)=1 * -ioo=5 => noeuds 3-1 = segment 3 appartient au contour elseif(ioo.eq.5) then ifis(i5,2)=4 * -ioo=6 => noeuds 2-3 = segment 2 appartient au contour elseif(ioo.eq.6) then ifis(i5,2)=2 * -cas ou 3 points appartiennent au contour => travail en + elseif(ioo.eq.7) then c write(6,*)'tri3 avec 2 segments dans le contour!',iel,i5 * rem: si les 3 points appartiennent au contour, on n'a pas de segments * appartenant au front ii1 = ipt5.num(1,i5) ii2 = ipt5.num(2,i5) ii3 = ipt5.num(3,i5) c il faut trouver quels sont les 2 segments appartenant au contour ! ia1 = -icpr(ii1) ia2 = -icpr(ii2) ia3 = -icpr(ii3) ia1min = min(ia1,ia2) ia1min = min(ia1min,ia3) ia1seg = 0 do 11 ia1=ia1min,nbelCON ifis1 = 0 do ib=1,nbnoCON if(ii1.eq.IPT6.num(ib,ia1)) ifis1=ifis1+1 if(ii2.eq.IPT6.num(ib,ia1)) ifis1=ifis1+2 if(ii3.eq.IPT6.num(ib,ia1)) ifis1=ifis1+4 enddo if(ifis1.le.2) goto 11 if(ifis1.eq.4) goto 11 if(ifis1.eq.3) ifis(i5,2)=ifis(i5,2)+1 if(ifis1.eq.6) ifis(i5,2)=ifis(i5,2)+2 if(ifis1.eq.5) ifis(i5,2)=ifis(i5,2)+4 if(ifis1.ge.7) then write(ioimp,*) 'tri3 avec 3 segments dans le contour !' call erreur(26) endif if(ifis1.eq.3.or.ifis1.ge.5) ia1seg=ia1seg+1 if(ia1seg.eq.2) goto 12 11 continue c on n a pas trouvé les 2 segments appartenant au contour write(ioimp,*)'on n a pas trouvé les 2 segments', $ ' du tri3 devant appartenant au contour !' call erreur(26) c on a trouvé les 2 segments appartenant au contour 12 continue if(ifis(i5,2).lt.3.or.ifis(i5,2).ge.7) then write(ioimp,*) 'impossible de detecter les 2 segments ', $ 'appartenant au contour avec les noeuds ',ii1,ii2,ii3 write(ioimp,*) 'ifis(ideb5,2)=',ifis(ideb5,2) write(ioimp,*) 'IPT6 pour les segments :' write(ioimp,*) ' ia1',ia1,(IPT6.num(iou,ia1),iou=1,2) write(ioimp,*) ' ia2',ia2,(IPT6.num(iou,ia2),iou=1,2) write(ioimp,*) ' ia3',ia3,(IPT6.num(iou,ia3),iou=1,2) call erreur(26) return endif endif endif enddo * fin de boucle sur les elements de la fissure ----------------- c if ((ifin5-ideb5).ne.1) then if (ideb5.ne.nbelem2.or.(ifin5-imil5).ne.1) then write(ioimp,*) 'pb avec le rearrangement de elements de fissure' write(ioimp,*) ideb5,nbelem2,imil5,ifin5 call erreur(26) return endif * si le front n est pas ferme, il faut avoir detecte ses 2 extremites * sinon il est probable que la ligne du front melfro ne soit pas * orientée comme le bord du maillage melfis if ((ivui2f.ne.1).and.(ivui2l.ne.1)) then write(ioimp,*) 'ERREUR: Les extremites du front (noeuds', $ i2firs,i2last,') ont été detectees ',ivui2f, 'et ',ivui2l $ ,' fois au lieu de 1 fois attendue' write(ioimp,*) 'Verifiez la coherence avec l orientation' $ ,'du maillage de la fissure fourni' $ ,' et INVErser le sens du front le cas échéant.' c write(6,*) 'front de fissure=' c CALL ECMAIL(ipt2,JENT2) c write(6,*) 'surface de fissure=' c CALL ECMAIL(ipt5,JENT5) call erreur(26) return endif c segdes,meleme,IPT6 segdes,meleme,IPT6,ipt2 meleme=ipt5 melfis=meleme segsup icpr endif *------ fin du Cas ou l'on calcule phi ET psi -------------------------- * sont actifs: melmai(=struct changée en poi1), ipt5=meleme=melfis, * sont inactifs: ipt2=melfro *********************************************************************** * 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 melfis=meleme melpoi=melmai *********************************************************************** * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT *********************************************************************** * pour chaque segment de la fissure (melfis 2eme maillage lu) * on calcule le champ phi et psi des noeuds autour de l'element if(iimpi.ge.1) then call gibtem(xkt) write(ioimp,*)'debut travail | xkt=',xkt endif * meleme=melfis c segact,meleme ipt3=melpoi c segact,ipt3 npt=ipt3.num(/2) * attention on initialise isup avec les valeurs inverses,il faudra * le retourner segini isup,isupfi,xdis,xdmin,xdmin1,rrlim if(ideux.ge.1) segini xdmin2,xdmin3 igard=0 igardf=0 *===================================================================== * Boucle sur les elements de la fissure ========================= do 400 iseg=1,nseg c coordonnes du tri3 ia=num(1,iseg) x11= xcoor((ia-1)*idim1 +1) y11= xcoor((ia-1)*idim1 +2) z11= xcoor((ia-1)*idim1 +3) ib=num(2,iseg) x22= xcoor((ib-1)*idim1 +1) y22= xcoor((ib-1)*idim1 +2) z22= xcoor((ib-1)*idim1 +3) ic=num(3,iseg) x33= xcoor((ic-1)*idim1 +1) y33= xcoor((ic-1)*idim1 +2) z33= xcoor((ic-1)*idim1 +3) if (x11.ge.x22) then xxmi=min(x22,x33) - xcrit xxma=max(x11,x33) + xcrit else xxmi=min(x11,x33) - xcrit xxma=max(x22,x33) + xcrit endif if (y11.ge.y22) then yymi=min(y22,y33) - xcrit yyma=max(y11,y33) + xcrit else yymi=min(y11,y33) - xcrit yyma=max(y22,y33) + xcrit endif if (z11.ge.z22) then zzmi=min(z22,z33) - xcrit zzma=max(z11,z33) + xcrit else zzmi=min(z11,z33) - xcrit zzma=max(z22,z33) + xcrit endif c vecteurs des aretes v12x=x22-x11 v12y=y22-y11 v12z=z22-z11 v23x=x33-x22 v23y=y33-y22 v23z=z33-z22 v31x=x11-x33 v31y=y11-y33 v31z=z11-z33 c vecteur normal vx=v12y*(z33-z11)-v12z*(y33-y11) vy=v12z*(x33-x11)-v12x*(z33-z11) vz=v12x*(y33-y11)-v12y*(x33-x11) bb=vx*vx + vy*vy+ vz*vz aa= sqrt(bb) if (aa.le.XPETIT) then write(ioimp,*) 'element fissure de longueur nulle ',aa write(ioimp,*) 'tri3.',iseg,'.',ia,'=',x11,y11,z11 write(ioimp,*) 'tri3.',iseg,'.',ib,'=',x22,y22,z22 write(ioimp,*) 'tri3.',iseg,'.',ic,'=',x33,y33,z33 call erreur(26) return endif sur = bb/2.d0 vx=vx/aa vy=vy/aa vz=vz/aa vxyz = max(abs(vx),abs(vy)) vxyz = max(vxyz,abs(vz)) if (vxyz.le.XPETIT) then write(ioimp,*)'element fissure de longueur nulle',aa,vx,vy,vz call erreur(26) return endif xsur = 1.D-8*sur c vecteurs des aretes (suite) v12lo= sqrt(v12x*v12x+v12y*v12y+v12z*v12z) v23lo= sqrt(v23x*v23x+v23y*v23y+v23z*v23z) v31lo= sqrt(v31x*v31x+v31y*v31y+v31z*v31z) v12x= v12x/v12lo v12y= v12y/v12lo v12z= v12z/v12lo v23x= v23x/v23lo v23y= v23y/v23lo v23z= v23z/v23lo v31x= v31x/v31lo v31y= v31y/v31lo v31z= v31z/v31lo if(iimpi.ge.1) $ write(ioimp,*)'-> triangle ',iseg,'/',nseg,' noeuds ',ia,ib,ic *======================================================================= * Boucle sur les noeuds ====================================== do 490 ipt=1,npt id=ipt3.num(1,ipt) dis = xgrand iphi2 = 0 *------------ cas ou le noeud structure id = l un des noeuds des tri3 de la fissure if (ia.eq.id.or.ib.eq.id.or.ic.eq.id) then dis=0.D0 iphi2 = 1 *------------ "autres" cas else * recup des coordonnes du noeud en question xxp = xcoor((id-1)*idim1 + 1 ) yyp = xcoor((id-1)*idim1 + 2 ) zzp = xcoor((id-1)*idim1 + 3 ) *bp2013 on saute ce point si trop loin (critere) if (icle.eq.1) then if(xxp.lt.xxmi) goto 491 if(xxp.gt.xxma) goto 491 if(yyp.lt.yymi) goto 491 if(yyp.gt.yyma) goto 491 if(zzp.lt.zzmi) goto 491 if(zzp.gt.zzma) goto 491 endif d1xp=xxp-x11 d1yp=yyp-y11 d1zp=zzp-z11 d2xp=xxp-x22 d2yp=yyp-y22 d2zp=zzp-z22 d3xp=xxp-x33 d3yp=yyp-y33 d3zp=zzp-z33 * calcul de la distance au triangle d1= sqrt(d1xp*d1xp + d1yp*d1yp + d1zp*d1zp) d2= sqrt(d2xp*d2xp + d2yp*d2yp + d2zp*d2zp) d3= sqrt(d3xp*d3xp + d3yp*d3yp + d3zp*d3zp) if (d1.lt.d2) then rr= min(d1,d3) else rr= min(d2,d3) endif * test pour eviter plein de calculs => on saute au noeud suivant if (isup(ipt).eq.0) then rrlim(ipt) = rr else c if(rr.gt.(1.1d0*rrlim(ipt))) goto 491 if(rr.gt.(2.d0*rrlim(ipt))) goto 491 if(rr.lt.rrlim(ipt)) rrlim(ipt)=rr endif * calcul de la distance au plan du triangle xdd = vx*d1xp + vy*d1yp + vz*d1zp * determination du pt projetté xx=xxp-xdd*vx yy=yyp-xdd*vy zz=zzp-xdd*vz * ce pt est-il à l'interieur ou sur un coté du triangle? d1x=xx-x11 d1y=yy-y11 d1z=zz-z11 d2x=xx-x22 d2y=yy-y22 d2z=zz-z22 d3x=xx-x33 d3y=yy-y33 d3z=zz-z33 pv12x=d1y*d2z-d1z*d2y pv12y=d1z*d2x-d1x*d2z pv12z=d1x*d2y-d1y*d2x pv23x=d2y*d3z-d2z*d3y pv23y=d2z*d3x-d2x*d3z pv23z=d2x*d3y-d2y*d3x pv31x=d3y*d1z-d3z*d1y pv31y=d3z*d1x-d3x*d1z pv31z=d3x*d1y-d3y*d1x ps12v=pv12X*vx+pv12y*vy+pv12z*vz ps23v=pv23x*vx+pv23y*vy+pv23z*vz ps31v=pv31x*vx+pv31y*vy+pv31z*vz * ce pt est-il à l'interieur ou sur un coté du triangle? if(ps12v.ge.0..and.ps23v.ge.0..and.ps31v.ge.0.) then * -> on est a l'interieur dis = xdd iphi2 = 1 elseif (abs(ps12v).lt.xsur) then * -> on est situé sur le segment 12 : est-on a l'interieur? psc1=v12x*d1x+v12y*d1y+v12z*d1z psc2=v12x*d2x+v12y*d2y+v12z*d2z if ( psc1*psc2 .lt.0.d0) then dis = xdd iphi2 = 1 endif elseif(abs(ps23v).lt.xsur)then * -> on est situé sur le segment 23 : est-on a l'interieur? psc2=v23x*d2x+v23y*d2y+v23z*d2z psc3=v23x*d3x+v23y*d3y+v23z*d3z if ( psc2*psc3 .lt.0.d0) then dis = xdd iphi2 = 1 endif elseif(abs(ps31v).lt.xsur)then * -> on est situé sur le segment 31 : est-on a l'interieur? psc3=v31x*d3x+v31y*d3y+v31z*d3z psc1=v31x*d1x+v31y*d1y+v31z*d1z if ( psc3*psc1 .lt.0.d0) then dis = xdd iphi2 = 1 endif endif * on est a l'exterieur : il faut calculer dis et verifier xdd if (dis.eq.xgrand) then iphi2 = 10 * si on est a l'ext du tri3 et en avant du front+contour on change v * -si on a affaire a un triangle du front (ou du contour) : * est on en avant du front (ou du contour)? -> phicor phicor=.false. phi12=.false. phi23=.false. phi31=.false. if (ideux.ge.1) then if (ifis(iseg,1).ge.1) then * -Cas ou le segment [1-2] = front ou contour if(ifis(iseg,1).ge.2.or.ifis(iseg,2).eq.1.or. $ ifis(iseg,2).eq.3.or.ifis(iseg,2).eq.5) then c phi12=(ps12v.lt.0.d0) phi12=(ps12v.lt.(-xsur)) if(phi12) then if(ifis(iseg,1).ge.2) then iphi2=iphi2+5 else iphi2=iphi2+1 endif endif endif * -Cas ou le segment [2-3] = contour if(ifis(iseg,2).eq.2.or. $ ifis(iseg,2).eq.3.or.ifis(iseg,2).eq.6) then c phi23=(ps23v.lt.0.d0) phi23=(ps23v.lt.(-xsur)) if(phi23) iphi2=iphi2+1 endif * -Cas ou le segment [3-1] = contour if(ifis(iseg,2).eq.4.or. $ ifis(iseg,2).eq.5.or.ifis(iseg,2).eq.6) then c phi31=(ps31v.lt.0.d0) phi31=(ps31v.lt.(-xsur)) if(phi31) iphi2=iphi2+1 endif * -Union des test phicor=(ideux.ge.1).and.(phi12.or.phi23.or.phi31) endif endif * si phicor, il faudra alors "corriger" phi * (et prendre xdd au lieu de dis) psc1=0.d0 psc2=0.d0 psc3=0.d0 c dis=sign(rr,xdd) * bp (05/03/2012) : on remplace ci dessus par formule + precise * --on teste autour du noeud 1-- if (rr.eq.d1) then if (ps12v.lt.0.d0.and.ps31v.lt.0d0) then * on est definitivement + proche du noeud 1 dis = sign(rr,xdd) else if (ps12v.lt.0.d0) then * on est du coté du segment 12 psc1 = v12x*d1x + v12y*d1y + v12z*d1z elseif (ps31v.lt.0.d0) then * on est du coté du segment 13 psc1 = -1.d0*(v31x*d1x + v31y*d1y + v31z*d1z) else * on ne devrait pas passer par la psc1=0.d0 endif if (psc1.le.0d0) then dis = sign(rr,xdd) else dis = sign(sqrt(rr*rr-psc1*psc1),xdd) endif endif * --on teste autour du noeud 2-- elseif (rr.eq.d2) then if (ps12v.lt.0.d0.and.ps23v.lt.0d0) then * on est definitivement + proche du noeud 2 dis = sign(rr,xdd) else if (ps12v.lt.0.d0) then * on est du coté du segment 21 psc2 = -1.d0*(v12x*d2x + v12y*d2y + v12z*d2z) elseif (ps31v.lt.0.d0) then * on est du coté du segment 23 psc2 = v23x*d2x + v23y*d2y + v23z*d2z else * on ne devrait pas passer par la psc2=0.d0 endif if (psc2.le.0d0) then dis = sign(rr,xdd) else dis = sign(sqrt(rr*rr-psc2*psc2),xdd) endif endif * --on teste autour du noeud 3-- elseif (rr.eq.d3) then if (ps31v.lt.0.d0.and.ps23v.lt.0d0) then * on est definitivement + proche du noeud 2 dis = sign(rr,xdd) else if (ps31v.lt.0.d0) then * on est du coté du segment 31 psc3 = v31x*d3x + v31y*d3y + v31z*d3z elseif (ps31v.lt.0.d0) then * on est du coté du segment 32 psc3 = -1.d0*(v23x*d3x + v23y*d3y + v23z*d3z) else * on ne devrait pas passer par la psc3=0.d0 endif if (psc3.le.0d0) then dis = sign(rr,xdd) else dis = sign(sqrt(rr*rr-psc3*psc3),xdd) endif endif endif endif * fin du cas ou on est a l'exterieur endif *------------ fin des "autres" cas (noeud id different de ceux du tri3) * -on a calculé (xdd et) dis: s'agit il de phi? c on prend xdd plutot que dis pour definir Phi dans le cas ou : c -on a l option DEUX ou TROI c -le tri3 utilisé (le + proche) appartient au front ou au contour c -le point se projette en avant du front ou du contour * iphi2 qu on va eventuellement ranger dans isup(ipt) vaut: * 1 si on se projette a l int du triangle * 10 si on se projette a l ext du triangle * 11 si en + on sort du contour * 15 si en + on est en avant du front * 16 si en + on est en avant du front et on sort du contour if(id.eq.idbug) then c if(phicor) then write(6,*) 'noeud:',id,' (',ipt,') tri3',iseg,':',ia,ib,ic write(6,*) 'iphi2,isup(ipt),dis,xdd',iphi2,isup(ipt),dis,xdd write(6,*) 'phi12,phicor=',phi12,phicor endif * -c'est la premiere fois que l'on voit ce point if (isup(ipt).eq.0) then if(phicor) then mpova2.vpocha(ipt,1)=xdd else mpova2.vpocha(ipt,1)=dis endif xdmin(ipt) =dis xdmin1(ipt)=xdd isup(ipt)=iphi2 iphi2=100 igard=igard+1 else * --on est (franchement) + proche que les precedents essais if (abs(dis).lt.0.999999*abs(xdmin(ipt))) then if(phicor) then mpova2.vpocha(ipt,1)=xdd else mpova2.vpocha(ipt,1)=dis endif xdmin(ipt) =dis xdmin1(ipt)=xdd isup(ipt)=iphi2 iphi2=100 * --on est (quasiment) aussi proche que le meilleur essai elseif (abs(dis).lt.1.000001*abs(xdmin(ipt))) then * -on a 2 tri3 du front+contour a egalite * et l'un des deux a trouvé le point en avant du front : if ((iseg.le.imil5).and. $ (phicor.or.isup(ipt).gt.10)) then c strategie : xdd= max(xdd1 xdd2) if (abs(xdd).gt.abs(xdmin1(ipt))) then mpova2.vpocha(ipt,1)=xdd xdmin(ipt) =dis xdmin1(ipt)=xdd c isup(ipt)=iphi2 if(iphi2.gt.isup(ipt)) isup(ipt)=iphi2 iphi2=100 c cas ou on avait pas vu qu il fallait prendre xdd elseif(isup(ipt).le.10) then c equivalent à : c elseif(mpova2.vpocha(ipt,1).ne.xdmin1(ipt)) then mpova2.vpocha(ipt,1)=xdmin1(ipt) isup(ipt)=iphi2 iphi2=100 endif * -on avait 1 tri3 du front+contour et ce tri3 est a egalite * et le 1er avait trouvé le point en avant du front : * on laisse le xdd calculé depuis le 1er tri3 elseif ((iseg.gt.imil5).and. $ (isup(ipt).gt.10)) then goto 491 * -quasi-egalite entre 2 concurrents mais signes opposés : else ss1=dis*xdmin(ipt) if (ss1.le.0.d0) then * on choisit le dis associé au xdd le + marqué if (abs(xdd).gt.abs(xdmin1(ipt))) then mpova2.vpocha(ipt,1)=dis xdmin(ipt) =dis xdmin1(ipt)=xdd c isup(ipt)=iphi2 if(iphi2.gt.isup(ipt)) isup(ipt)=iphi2 iphi2=100 endif endif endif endif * --fin de la distinction + proche / aussi proche endif * -fin de la distinction 1ere fois ou pas if(id.eq.idbug) then write(6,*)'iphi2,xdmin,xdmin1,phi',iphi2, $ xdmin(ipt),xdmin1(ipt),mpova2.vpocha(ipt,1) endif if(iimpi.ge.2) then write(ioimp,*) 'point',id,ipt,'dis,xdd,PHI,isup=', $ dis,xdd,mpova2.vpocha(ipt,1),isup(ipt) endif 491 continue * le tri3 n'est pas un tri3 du front => on saute *TODO : test a verifier peut etre .eq. et .eq.-1 aussi a ne pas sauter? if (ideux.lt.1) goto 490 if (ifis(iseg,1).le.1) goto 490 c ------ Calcul de la 2eme Level Set ------- xdd1 = 0.D0 * -cas ou le noeud structure id = l un des noeuds des tri3 de la fissure if (ia.eq.id.or.ib.eq.id) then dis1=0.D0 xd2=0.D0 xd3=0.D0 if(ideux.ge.2) xtau=xicpr(id) * -autres cas else * on recupere ce qu on a pas recupéré lors du calcul de phi if (ic.eq.id) then xxp = xcoor((id-1)*idim1 + 1 ) yyp = xcoor((id-1)*idim1 + 2 ) zzp = xcoor((id-1)*idim1 + 3 ) endif * calcul du vecteur vn qui est le sens de propagation vnx=v12y*vz-v12z*vy vny=v12z*vx-v12x*vz vnz=v12x*vy-v12y*vx * il faut orienter vnx pour que le vecteur soit dans le sens * de la propagation qq= vnx*(x33-x11)+vny*(y33-y11)+vnz*(z33-z11) if (qq.ge.0.D0) then vnx=-vnx vny=-vny vnz=-vnz endif * normalement vn est normé on le verifie qq= sqrt(vnx*vnx + vny*vny + vnz*vnz) if (abs(qq-1.d0).gt.1.d-2) then write(ioimp,*) ' pb dans le calcul de vn' call erreur(26) return elseif(abs(qq-1.).gt.1.d-6) then vnx=vnx/qq vny=vny/qq vnz=vnz/qq endif * on calcule la distance au plan defini par segment 12 et vecteur vx,vy,vz * et le point (xx yy zz ) projeté du point sur ce plan xdd1=(xxp-x11)*vnx+(yyp-y11)*vny+(zzp-z11)*vnz xx=xxp-xdd1*vnx yy=yyp-xdd1*vny zz=zzp-xdd1*vnz psc1=(xx-x11)*v12x+(yy-y11)*v12y+(zz-z11)*v12z psc2=(xx-x22)*v12x+(yy-y22)*v12y+(zz-z22)*v12z if ((psc1*psc2).le.0.d0) then C le point se projete sur 12 dis1=xdd1 pp1x=x11+(psc1*v12x) pp1y=y11+(psc1*v12y) pp1z=z11+(psc1*v12z) xd2 = 0.D0 xd3 = sqrt((xxp-pp1x)**2+(yyp-pp1y)**2+(zzp-pp1z)**2) if (ideux.ge.2) then xtau1=xicpr(ia) xtau2=xicpr(ib) stau = psc1 / v12lo xtau=stau*xtau2+(1.d0-stau)*xtau1 endif else if (psc1.ge.0.D0) then pp2x=xx+((v12lo-psc1)*v12x) pp2y=yy+((v12lo-psc1)*v12y) pp2z=zz+((v12lo-psc1)*v12z) xd2 = abs(v12lo-psc1) xd3 = sqrt((xxp-x22)**2+(yyp-y22)**2+(zzp-z22)**2) c if(ideux.ge.2) xtau=xicpr(ib) if(ideux.ge.2) then if(ib.eq.i2firs.or.ib.eq.i2last) then c write(6,*) 'le noeud',id,'detecte comme apres',ib xtau1=xicpr(ia) xtau2=xicpr(ib) stau = psc1 / v12lo xtau=stau*xtau2+(1.d0-stau)*xtau1 else xtau=xicpr(ib) endif endif else pp2x=xx-(psc1*v12x) pp2y=yy-(psc1*v12y) pp2z=zz-(psc1*v12z) xd2 = abs(psc1) xd3 = sqrt((xxp-x11)**2+(yyp-y11)**2+(zzp-z11)**2) c if(ideux.ge.2) xtau=xicpr(ia) if(ideux.ge.2) then if(ia.eq.i2firs.or.ia.eq.i2last) then c write(6,*) 'le noeud',id,'detecte comme avant',ia xtau1=xicpr(ia) xtau2=xicpr(ib) stau = psc1 / v12lo xtau=stau*xtau2+(1.d0-stau)*xtau1 else xtau=xicpr(ia) endif endif endif dis1=sqrt((xxp-pp2x)**2+(yyp-pp2y)**2+(zzp-pp2z)**2) dis1=sign(dis1,xdd1) endif endif * -fin de la distinction si noeud du tri3 ou pas * 1ere fois qu on traite ce point? if (isupfi(ipt).eq.0) then * isupfi(ipt)=1 * on utilise isupfi pour definir si la valeur inscrite * provient du 1er ou dernier seg du front (=2) ou autre(=1)? isupfi(ipt)=ifis(iseg,1) igardf=igardf+1 xdmin2(ipt)= xd2 xdmin3(ipt)= xd3 cbp2013 mpova1.vpocha(ipt,1)=dis1 if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau else * if(abs(dis1).lt.abs(mpova1.vpocha(ipt,1)))then * on a trouvé un segment + proche ou = du noeud considéré * if(abs(dis1).le.abs(mpova1.vpocha(ipt,1)))then ** if(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))-1.D-8)then * if(xd2.le.xdmin2(ipt)) then if (xd3.lt.(abs(xdmin3(ipt))-1.D-8)) then cbp2013 mpova1.vpocha(ipt,1)=dis1 xdmin2(ipt)=xd2 xdmin3(ipt)=xd3 isupfi(ipt)=ifis(iseg,1) if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau endif * rattrapage dans le cas d'erreur d'arrondi sur dis1 * elseif(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))+1.D-8) * $ then * if(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))+1.D-8)then if (xd3.lt.(abs(xdmin3(ipt))+1.D-8))then * en cas d'égalité,on n'utilise la distance dis1 que si: * - le pt projetté appartient au segment (xd2=0) * ou - le point projeté est + proche du segment if (xd2.lt.xdmin2(ipt)) then cbp2013 mpova1.vpocha(ipt,1)=dis1 xdmin2(ipt)=xd2 xdmin3(ipt)=xd3 isupfi(ipt)=ifis(iseg,1) if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau endif endif endif 490 continue * Fin de Boucle sur les noeuds =============================== *======================================================================= 400 continue * Fin de Boucle sur les elements de la fissure ========================= *======================================================================= *------Petite correction de Psi if (ideux.ge.1) then * Boucle sur les noeuds =============================== do 500 ipt=1,npt if (isup(ipt).eq.0) goto 500 c * pour le cas de fissure debouchante c * dans le cas ou c'est un segment debouchant qui a servi au calcul de psi, c * il faut prendre la distance xdd2 (et pas dis) qu'on retrouve avec : c if ((isupfi(ipt)).eq.2) then c dis1 = mpova1.vpocha(ipt,1) c xdd1 = sqrt( (dis1**2) - (xdmin2(ipt)**2) ) c mpova1.vpocha(ipt,1) = sign(xdd1,dis1) c endif c bp : essai du 16/03/2012 : psi = distance au front**2-phi**2 xdd1 = (xdmin3(ipt)**2) - (mpova2.vpocha(ipt,1)**2) xdd1 = sqrt(max(xdd1,0.d0)) if(ipt3.num(1,ipt).eq.idbug) then write(6,*)'ipt,isup,psi_avant,phi,xdd1,= ',ipt,isup(ipt), $ xdmin3(ipt),mpova1.vpocha(ipt,1),mpova2.vpocha(ipt,1),xdd1 endif c if (isup(ipt).gt.0) then if (isup(ipt).ge.15) then *on est a l exterieur => psi >0 mpova1.vpocha(ipt,1) = xdd1 else mpova1.vpocha(ipt,1) = -1.d0*xdd1 endif 500 continue endif *------Fin de la Petite correction de Psi *********************************************************************** * FIN DU TRAVAIL *********************************************************************** if(iimpi.ge.1) then call gibtem(xkt) write(ioimp,*)'fin du travail | xkt=',xkt endif * on inverse isup * AVANT : isup(i) =0 veut dire que le point na pas été traité * APRES : isup(i) =1 veut dire que le point peut etre enleve (=0 sinon) iaju=0 do ia=1,isup(/1) 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 ipt1=ipt3 ipt4=melfis segdes ipt4 * iaju est le nombre de points a éliminer * isup(i) =1 veut dire que le point peut etre enleve (=0 sinon) * faut-il ajuster les champs? if (iaju.ne.0) then * si oui on recree les bons chpoints 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 segdes msoup3,mpova6 segsup,mpova3 endif if (ideux.ge.1) then c segact msoup1*mod msoup1.ipoval=mpova4 msoup1.igeoc=ipt2 segdes msoup1,mpova4 segsup,mpova1 endif c segact msoup2*mod msoup2.ipoval=mpova5 msoup2.igeoc=ipt2 segdes msoup2,mpova5,ipt2 segsup,mpova2 segdes ipt1 else if(ideux.ge.2) segdes,msoup3,mpova3 if(ideux.ge.1) segdes,msoup1,mpova1 segdes msoup2,mpova2,ipt1 endif if (ideux.ge.1) then segsup,icpr,sfis if(ideux.ge.2) segsup,xicpr endif if (icle.eq.1) then segsup isup,isupfi,xdis,xdmin,xdmin1,rrlim if(ideux.ge.1) segsup xdmin2,xdmin3 endif RETURN END