C PSI3D0 SOURCE PV 21/12/18 07:15:11 11240 SUBROUTINE PSI3D0(ithr,IPOINT) *********************************************************************** * Calcul effectif de la (des) distance(s) signee(s) * Appele par PSIP3D ou PSIP3Di (//isation) *********************************************************************** ************************************************************************ *** CALCUL ELEMENTAIRE (PAR RAPPORT AU TRIANGLE abc) *** *** DES LEVELS SET AU NOEUD id *** *** explication par Benoit Prabel, 2019 *** * * on cherche : * dis = distance entre le noeud et le triangle * phi = level set normale au plan * psi = level set distance au front * * * ********** * * vue 3D * * ********** * x{ii},y{ii},z{ii} = coordonnees des noeuds{i} du triangle * xxp , yyp , zzp = coordonnees du noeud ou calculer les level set * v = vecteur normal au triangle = v12 ∧ -v31 = v12 ∧ v13 * x = noeud p projete sur le plan du triangle * * (vx) * -v (vy) (xxp) * ^ (vz) +id (yyp) * (x33) | (x11) | (zzp) * (y33) ic +------|-----+ ia (y11) | * (z33) \ | / (z11) ~ |xdd ~ ~ ~ ~ ~ * \ / / | / * \ / / | (xx) / * \ / / +x (yy) / * + (x22) / (zz) / * ib (y22) / / * (z22) / plan du triangle / * ~ ~ ~ ~ ~ ~ ~ ~ * * *************** * * vue de cote * * *************** * * d{i} = distance noeud {i} du triangle - point * rr = min(d1,d2,d3) * * + p * / | * d1 / | dessin identique pour les distances d2 et d3 * / | aux noeuds 2 et 3 * / | * / | * - - 1 + - - - - - +x - - - Plan du triangle * | * v v * * ************************************************* * * vue de dessus (i.e. dans le plan du triangle) * * ************************************************* * * v{ij} = vecteur unitaire du cote du triangle {i}-{j} (avec j=i+1) * v{i}x = vecteur du noeud {i} au noeud x * sur = surface du triangle * pv{i}{j} = v{i}x ∧ v{j}x * ps{ij}v = pv{i}{j} . v * = produit_mixte(v{i}x, v{j}x, v) * < 0 si x appartient au demi-plan defini par la droite (i-j) * du cote oppose au triangle * pour le cote i-j (avec j=i+1), on definit : * psc{iij} = v{i}x . v{ij} * psc{jij} = v{j}x . v{ij} * * (v1xx) * v31 v1x=(v1yy) * 3+----------->+1 _ _ (v1zz) * v23\ / - - _ _ * \ / - -> + x * \ / * \ / v12 * + * 2 * * DONC : * * ps{ij}v > 0 qqsoit i,j <=> x a l'interieur du triangle 123 * * ps{ij}v = 0 <=> x appartient a la droite (i-j) * * ps{ij}v = 0 ET (psc{iij} * psc{jij}) < 0 <=> x appartient a [i-j] * * * si x est a l'exterieur du triangle, * on regarde de quel noeud il est le plus proche (disons {i}) * et on calcule : * psc{iij} (tq j=i+1) * psc{iki} (tq k=i-1) * AINSI : * * psc{iij}<0 ET psc{iki}>0 <=> dis = d{i} * * ps{ij}v<0 ET psc{iij}>0 ET psc{jij}>0 * <=> dis = sqrt(d{i}**2-psc{iij}**2) = sqrt(d{j}**2-psc{jij}**2) * * pas d'autres cas en theorie * ==> on a calcule dis ! * * ... a completer pour psi et phi ... * ************************************************************************ IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SMCHPOI -INC SMELEME -INC SMCOORD -INC CCASSIS SEGMENT MTRI3 INTEGER KNODE(3,NTRI3) REAL*8 XNODE(3,NTRI3),YNODE(3,NTRI3),ZNODE(3,NTRI3) c REAL*8 XMIMA(NTRI3),YMIMA(NTRI3),ZMIMA(NTRI3) REAL*8 VARE12(4,NTRI3),VARE23(3,NTRI3),VARE31(3,NTRI3) REAL*8 VNORM(3,NTRI3),SURFAC(NTRI3) ENDSEGMENT segment mxicpr real*8 xicpr(nbpts) endsegment segment isup(npt) segment sfis integer ifis(nbelfi,2) endsegment SEGMENT SPARAL INTEGER NBTHRD INTEGER IERROR(NBTHR) INTEGER KPT3,KPOVA1,KPOVA2,KPOVA3 INTEGER KTRI3,KXICPR,KSUP,KFIS INTEGER KCLE,KDEUX,KMIL5 REAL*8 XCRIT1 ENDSEGMENT LOGICAL phicor,phi12,phi23,phi31 SPARAL = IPOINT NBTHR = SPARAL.NBTHRD IPT3 = SPARAL.KPT3 mpova1 = SPARAL.KPOVA1 mpova2 = SPARAL.KPOVA2 mpova3 = SPARAL.KPOVA3 MTRI3 = SPARAL.KTRI3 mxicpr = SPARAL.KXICPR isup = SPARAL.KSUP sfis = SPARAL.KFIS icle = SPARAL.KCLE ideux = SPARAL.KDEUX xcrit = SPARAL.XCRIT1 imil5 = SPARAL.KMIL5 phicor=.false. idim1 = IDIM+1 npt=ipt3.num(/2) nseg=MTRI3.KNODE(/2) c idbug=1454 c idbug=1 c idbug=-3 C Decoupage pour le travail d'ecriture en parallele IRES=MOD(npt,NBTHR) IF (IRES.EQ.0) THEN ILON = npt / NBTHR IDEB = (ithr -1)* ILON + 1 ELSE IF (ithr .LE. IRES) THEN ILON = (npt / NBTHR) + 1 IDEB = (ithr -1)* ILON + 1 ELSE ILON = npt / NBTHR IDEB = (IRES * (ILON+1)) + (ithr-IRES-1)* ILON + 1 ENDIF ENDIF IFIN = IDEB + ILON - 1 c print *, "PSI3D0 ith=",ithr,"/",NBTHR," ==>",IDEB,IFIN,"/",npt *======================================================================= * Boucle sur les noeuds ====================================== c do 490 ipt=1,npt do 490 ipt=IDEB,IFIN isupfi=0 igard=0 igardf=0 id=ipt3.num(1,ipt) c print *, "NOEUD",ipt,"/",npt," #",id * recup des coordonnees du noeud en question xxp = xcoor((id-1)*idim1 + 1 ) yyp = xcoor((id-1)*idim1 + 2 ) zzp = xcoor((id-1)*idim1 + 3 ) *===================================================================== * Boucle sur les elements de la fissure ========================= do 400 iseg=1,nseg dis = xgrand iphi2 = 0 c Recup des infos de la fissure (TRI3) ia=MTRI3.KNODE(1,iseg) ib=MTRI3.KNODE(2,iseg) ic=MTRI3.KNODE(3,iseg) c if(id.eq.idbug) print *,"TRI",iseg,"/",nseg," #",ia,ib,ic * recup des coordonnees du triangle x11=MTRI3.XNODE(1,iseg) x22=MTRI3.XNODE(2,iseg) x33=MTRI3.XNODE(3,iseg) y11=MTRI3.YNODE(1,iseg) y22=MTRI3.YNODE(2,iseg) y33=MTRI3.YNODE(3,iseg) z11=MTRI3.ZNODE(1,iseg) z22=MTRI3.ZNODE(2,iseg) z33=MTRI3.ZNODE(3,iseg) * recup de la normale au triangle + vecteurs des aretes + surface vx=MTRI3.VNORM(1,iseg) vy=MTRI3.VNORM(2,iseg) vz=MTRI3.VNORM(3,iseg) v12x=MTRI3.VARE12(1,iseg) v12y=MTRI3.VARE12(2,iseg) v12z=MTRI3.VARE12(3,iseg) v12lo=MTRI3.VARE12(4,iseg) v23x=MTRI3.VARE23(1,iseg) v23y=MTRI3.VARE23(2,iseg) v23z=MTRI3.VARE23(3,iseg) v31x=MTRI3.VARE31(1,iseg) v31y=MTRI3.VARE31(2,iseg) v31z=MTRI3.VARE31(3,iseg) sur =MTRI3.SURFAC(iseg) xsur=1.D-8*sur *------------ 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 phicor=.false. *------------ "autres" cas else *bp2013 on saute ce point si trop loin (critere) if (icle.eq.1) then 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 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**2 + d1yp**2 + d1zp**2) d2= sqrt(d2xp**2 + d2yp**2 + d2zp**2) d3= sqrt(d3xp**2 + d3yp**2 + d3zp**2) if (d1.lt.d2) then rr= min(d1,d3) else rr= min(d2,d3) 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? v1xx=xx-x11 v1yy=yy-y11 v1zz=zz-z11 v2xx=xx-x22 v2yy=yy-y22 v2zz=zz-z22 v3xx=xx-x33 v3yy=yy-y33 v3zz=zz-z33 pv12x=v1yy*v2zz-v1zz*v2yy pv12y=v1zz*v2xx-v1xx*v2zz pv12z=v1xx*v2yy-v1yy*v2xx pv23x=v2yy*v3zz-v2zz*v3yy pv23y=v2zz*v3xx-v2xx*v3zz pv23z=v2xx*v3yy-v2yy*v3xx pv31x=v3yy*v1zz-v3zz*v1yy pv31y=v3zz*v1xx-v3xx*v1zz pv31z=v3xx*v1yy-v3yy*v1xx ps12v=pv12x*vx+pv12y*vy+pv12z*vz ps23v=pv23x*vx+pv23y*vy+pv23z*vz ps31v=pv31x*vx+pv31y*vy+pv31z*vz c if(id.eq.idbug) then c write(6,*) '-----------------------------------------------' c write(6,*) 'noeud:',id,' (',ipt,') tri3',iseg,':',ia,ib,ic c write(6,*) ' (',xx,') (',x11,') (',x22,') (',x33,')' c write(6,*) ' x=(',yy,') 1:(',y11,') 2:(',y22,') 3:(',y33,')' c write(6,*) ' (',zz,') (',z11,') (',z22,') (',z33,')' c write(6,*)'rr=d1,2,3?',(rr.eq.d1),(rr.eq.d2),(rr.eq.d3),rr c write(6,*) 'xdd=',xdd c write(6,*) 'ps12v,23,31=',ps12v,ps23v,ps31v c endif * 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? psc112=v1xx*v12x+v1yy*v12y+v1zz*v12z psc212=v2xx*v12x+v2yy*v12y+v2zz*v12z if ( psc112*psc212 .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? psc223=v2xx*v23x+v2yy*v23y+v2zz*v23z psc323=v3xx*v23x+v3yy*v23y+v3zz*v23z if ( psc223*psc323 .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? psc331=v3xx*v31x+v3yy*v31y+v3zz*v31z psc131=v1xx*v31x+v1yy*v31y+v1zz*v31z if ( psc331*psc131 .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) c if(id.eq.idbug) print *,"phicor =",phicor c psc1=0.d0 c psc2=0.d0 c psc3=0.d0 c dis=sign(rr,xdd) * bp (05/03/2012) : on remplace ci dessus par formule + precise * bp,gg (2019-09-25) : correction pour triangle quelconque * produits scalaires determinants si on est proche * d'un noeud ou d'un segment psc112=v1xx*v12x+v1yy*v12y+v1zz*v12z psc131=v1xx*v31x+v1yy*v31y+v1zz*v31z psc223=v2xx*v23x+v2yy*v23y+v2zz*v23z psc212=v2xx*v12x+v2yy*v12y+v2zz*v12z psc331=v3xx*v31x+v3yy*v31y+v3zz*v31z psc323=v3xx*v23x+v3yy*v23y+v3zz*v23z * --on teste si on est dans un angle du triangle-- if (psc112.le.0.d0.and.psc131.ge.0d0) then * on est definitivement + proche du noeud 1 dis = sign(rr,xdd) elseif (psc223.le.0.d0.and.psc212.ge.0d0) then * on est definitivement + proche du noeud 2 dis = sign(rr,xdd) elseif (psc331.le.0.d0.and.psc323.ge.0d0) then * on est definitivement + proche du noeud 3 dis = sign(rr,xdd) * --on teste quel segment est le plus proche-- elseif(ps12v.le.0.d0.and.psc112.ge.0.d0 & .and.psc212.le.0) then * on est du coté du segment 12 dis = sign(sqrt(MAX(d1**2-psc112**2,0.D0)),xdd) elseif(ps23v.le.0.d0.and.psc223.ge.0.d0 & .and.psc323.le.0) then * on est du coté du segment 23 dis = sign(sqrt(MAX(d2**2-psc223**2,0.D0)),xdd) elseif(ps31v.le.0.d0.and.psc331.ge.0.d0 & .and.psc131.le.0) then * on est du coté du segment 31 dis = sign(sqrt(MAX(d3**2-psc331**2,0.D0)),xdd) else c on ne devrait pas arriver ici ! print *,"psi3d0 : noeud",id," tri",ia,ib,ic print *,psc112,psc131,psc223,psc212,psc331,psc323 c call erreur(5) SPARAL.IERROR(ithr)=5 return 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 c if(id.eq.idbug) then c write(6,*) 'dis,xdd',dis,xdd c write(6,*) 'phicor=',phicor c write(6,*) 'phi12,phi23,phi31=',phi12,phi23,phi31 c 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 =dis xdmin1=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)) then if(phicor) then mpova2.vpocha(ipt,1)=xdd else mpova2.vpocha(ipt,1)=dis endif xdmin =dis xdmin1=xdd isup(ipt)=iphi2 iphi2=100 * --on est (quasiment) aussi proche que le meilleur essai elseif (abs(dis).lt.1.000001*abs(xdmin)) 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)) then mpova2.vpocha(ipt,1)=xdd xdmin =dis xdmin1=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 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 et possibles signes opposés : else c bp: nouvelle strategie: on moyenne les normales*(y-xp) (soit xdd) xdmin1=xdmin1+xdd mpova2.vpocha(ipt,1)=sign(dis,xdmin1) endif endif * --fin de la distinction + proche / aussi proche endif * -fin de la distinction 1ere fois ou pas c if(id.eq.idbug) then c write(6,*)'iphi2,xdmin,xdmin1,phi',iphi2, c $ xdmin,xdmin1,mpova2.vpocha(ipt,1) c endif c if(iimpi.ge.2) then c write(ioimp,*) 'point',id,ipt,'dis,xdd,PHI,isup=', c $ dis,xdd,mpova2.vpocha(ipt,1),isup(ipt) c endif * point d'arrivee si hors-boite 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 400 if (ifis(iseg,1).le.1) goto 400 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 c * on recupere ce qu on a pas recupéré lors du calcul de phi c if (ic.eq.id) then c xxp = xcoor((id-1)*idim1 + 1 ) c yyp = xcoor((id-1)*idim1 + 2 ) c zzp = xcoor((id-1)*idim1 + 3 ) c * + coordonnees du triangle+... c endif cbp, desormais inutile car deplace + haut * 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.gt.0.D0) then vnx=-vnx vny=-vny vnz=-vnz endif * normalement vn est normé on le verifie qq= sqrt(vnx**2 + vny**2 + vnz**2) if (abs(qq-1.d0).gt.1.d-2) then PRINT *,'PSIP3D0: pb dans le calcul de vn',ithr c call erreur(26) SPARAL.IERROR(ithr)=26 return elseif(abs(qq-1.d0).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 c if(id.eq.idbug) print *,"isupfi =",isupfi * 1ere fois qu on traite ce point? if (isupfi.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=ifis(iseg,1) igardf=igardf+1 xdmin2=xd2 xdmin3=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)-1.D-8)) then cbp2013 mpova1.vpocha(ipt,1)=dis1 xdmin2=xd2 xdmin3=xd3 isupfi=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)+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) then cbp2013 mpova1.vpocha(ipt,1)=dis1 xdmin2=xd2 xdmin3=xd3 isupfi=ifis(iseg,1) if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau endif endif endif 400 continue * Fin de Boucle sur les elements de la fissure =============== *======================================================================= *------ Petite correction de Psi if (ideux.ge.1) then 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(MAX(dis1**2 - xdmin2(ipt)**2,0.D0)) 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**2) - (mpova2.vpocha(ipt,1)**2) xdd1 = sqrt(max(xdd1,0.d0)) c if(ipt3.num(1,ipt).eq.idbug) then c write(6,*)'ipt,isup,psi_avant,phi,xdd1,= ',ipt,isup(ipt), c $ xdmin3,mpova1.vpocha(ipt,1),mpova2.vpocha(ipt,1),xdd1 c 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 490 continue * Fin de Boucle sur les noeuds =================================== *======================================================================= END