Télécharger la plaquette Cast3M


fr en

Télécharger psi3d0.eso

Retour à la liste

Numérotation des lignes : 

psi3d0

  1. C PSI3D0    SOURCE    PV        21/12/18    07:15:11     11240          
  2.  
  3.       SUBROUTINE PSI3D0(ithr,IPOINT)
  4.  
  5. ***********************************************************************
  6. *     Calcul effectif de la (des) distance(s) signee(s)
  7. *     Appele par PSIP3D ou PSIP3Di (//isation)
  8. ***********************************************************************
  9.  
  10. ************************************************************************ 
  11. ***        CALCUL ELEMENTAIRE (PAR RAPPORT AU TRIANGLE abc)          ***
  12. ***        DES LEVELS SET AU NOEUD id                                ***
  13. ***        explication par Benoit Prabel, 2019                       ***
  14. *
  15. *    on cherche :
  16. *    dis = distance entre le noeud et le triangle
  17. *    phi = level set normale au plan
  18. *    psi = level set distance au front
  19. *    
  20. *  
  21. *  **********
  22. *  * vue 3D * 
  23. *  **********
  24. *    x{ii},y{ii},z{ii} = coordonnees des noeuds{i} du triangle 
  25. *    xxp , yyp , zzp   = coordonnees du noeud ou calculer les level set
  26. *    v = vecteur normal au triangle = v12 ∧ -v31 = v12 ∧ v13
  27. *    x = noeud p projete sur le plan du triangle
  28. *  
  29. *                          (vx)
  30. *                       -v (vy)                         (xxp) 
  31. *                      ^   (vz)                     +id (yyp) 
  32. *      (x33)           |          (x11)             |   (zzp)
  33. *      (y33) ic +------|-----+ ia (y11)             |
  34. *      (z33)     \     |    /     (z11)           ~ |xdd   ~  ~  ~  ~ ~
  35. *                  \       /                    /   |                 /
  36. *                    \    /                    /    |  (xx)          /
  37. *                      \ /                    /     +x (yy)         /
  38. *                       +  (x22)             /         (zz)        /
  39. *                       ib (y22)            /                     /
  40. *                          (z22)           /  plan du triangle   / 
  41. *                                          ~  ~  ~  ~  ~  ~  ~ ~
  42. *                     
  43. *  ***************
  44. *  * vue de cote *
  45. *  ***************
  46. *  
  47. *    d{i} = distance noeud {i} du triangle - point
  48. *    rr   = min(d1,d2,d3)
  49. *    
  50. *                      + p 
  51. *                    / |      
  52. *             d1   /   |      dessin identique pour les distances d2 et d3 
  53. *                /     |      aux noeuds 2 et 3
  54. *              /       |      
  55. *            /         |                
  56. *    - - 1 + - - - - - +x - - - Plan du triangle                    
  57. *                |
  58. *                v v 
  59. *                
  60. *  ************************************************* 
  61. *  * vue de dessus (i.e. dans le plan du triangle) *
  62. *  *************************************************
  63. *  
  64. *    v{ij}     = vecteur unitaire du cote du triangle {i}-{j} (avec j=i+1)
  65. *    v{i}x     = vecteur du noeud {i} au noeud x
  66. *    sur       = surface du triangle
  67. *    pv{i}{j}  = v{i}x ∧ v{j}x 
  68. *    ps{ij}v   = pv{i}{j} . v 
  69. *              = produit_mixte(v{i}x, v{j}x, v)
  70. *              < 0 si x appartient au demi-plan defini par la droite (i-j) 
  71. *                  du cote oppose au triangle
  72. *    pour le cote i-j (avec j=i+1), on definit :
  73. *    psc{iij}    = v{i}x . v{ij}  
  74. *    psc{jij}    = v{j}x . v{ij}  
  75. *    
  76. *                                      (v1xx)
  77. *                        v31       v1x=(v1yy)
  78. *             3+----------->+1 _ _     (v1zz)
  79. *            v23\          /       - - _ _      
  80. *                 \       /                - -> + x 
  81. *                   \    /                 
  82. *                     \ / v12  
  83. *                      +
  84. *                      2
  85. *                      
  86. *    DONC : 
  87. *    * ps{ij}v > 0 qqsoit i,j  <=>  x a l'interieur du triangle 123
  88. *    * ps{ij}v = 0  <=>  x appartient a la droite (i-j)
  89. *    * ps{ij}v = 0  ET  (psc{iij} * psc{jij}) < 0  <=> x appartient a [i-j]
  90. *  
  91. *    * si x est a l'exterieur du triangle, 
  92. *      on regarde de quel noeud il est le plus proche (disons {i}) 
  93. *      et on calcule :
  94. *      psc{iij} (tq j=i+1)
  95. *      psc{iki} (tq k=i-1)
  96. *      AINSI :
  97. *      * psc{iij}<0 ET psc{iki}>0 <=> dis = d{i}
  98. *      * ps{ij}v<0  ET psc{iij}>0 ET psc{jij}>0 
  99. *        <=> dis = sqrt(d{i}**2-psc{iij}**2) = sqrt(d{j}**2-psc{jij}**2)
  100. *      * pas d'autres cas en theorie
  101. *      ==> on a calcule dis !
  102. *
  103. *     ... a completer pour psi et phi ...
  104. *  
  105. ************************************************************************ 
  106.  
  107.       IMPLICIT INTEGER(I-N)
  108.       IMPLICIT REAL*8 (A-H,O-Z)
  109.  
  110.  
  111. -INC PPARAM
  112. -INC CCOPTIO
  113. -INC CCREEL
  114. -INC SMCHPOI
  115. -INC SMELEME
  116. -INC SMCOORD
  117. -INC CCASSIS
  118.  
  119.       SEGMENT MTRI3
  120.         INTEGER KNODE(3,NTRI3)
  121.         REAL*8  XNODE(3,NTRI3),YNODE(3,NTRI3),ZNODE(3,NTRI3)
  122. c         REAL*8  XMIMA(NTRI3),YMIMA(NTRI3),ZMIMA(NTRI3)
  123.         REAL*8  VARE12(4,NTRI3),VARE23(3,NTRI3),VARE31(3,NTRI3)
  124.         REAL*8  VNORM(3,NTRI3),SURFAC(NTRI3)
  125.       ENDSEGMENT
  126.  
  127.       segment mxicpr
  128.         real*8 xicpr(nbpts)
  129.       endsegment
  130.       segment isup(npt)
  131.       segment sfis
  132.          integer ifis(nbelfi,2)
  133.       endsegment
  134.  
  135.       SEGMENT SPARAL
  136.         INTEGER NBTHRD
  137.         INTEGER IERROR(NBTHR)
  138.         INTEGER KPT3,KPOVA1,KPOVA2,KPOVA3
  139.         INTEGER KTRI3,KXICPR,KSUP,KFIS
  140.         INTEGER KCLE,KDEUX,KMIL5
  141.         REAL*8  XCRIT1
  142.       ENDSEGMENT
  143.  
  144.       LOGICAL phicor,phi12,phi23,phi31
  145.  
  146.       SPARAL = IPOINT
  147.       NBTHR  = SPARAL.NBTHRD
  148.       IPT3   = SPARAL.KPT3
  149.       mpova1 = SPARAL.KPOVA1
  150.       mpova2 = SPARAL.KPOVA2
  151.       mpova3 = SPARAL.KPOVA3
  152.       MTRI3  = SPARAL.KTRI3
  153.       mxicpr = SPARAL.KXICPR
  154.       isup   = SPARAL.KSUP
  155.       sfis   = SPARAL.KFIS
  156.       icle   = SPARAL.KCLE
  157.       ideux  = SPARAL.KDEUX
  158.       xcrit  = SPARAL.XCRIT1
  159.       imil5  = SPARAL.KMIL5
  160.  
  161.       phicor=.false.
  162.       idim1 = IDIM+1
  163.       npt=ipt3.num(/2)
  164.       nseg=MTRI3.KNODE(/2)
  165. c       idbug=1454
  166. c       idbug=1
  167. c       idbug=-3
  168.  
  169. C     Decoupage pour le travail d'ecriture en parallele
  170.       IRES=MOD(npt,NBTHR)
  171.       IF (IRES.EQ.0) THEN
  172.         ILON  = npt / NBTHR
  173.         IDEB  = (ithr -1)* ILON + 1
  174.       ELSE
  175.         IF (ithr .LE. IRES) THEN
  176.           ILON  = (npt / NBTHR) + 1
  177.           IDEB  = (ithr -1)* ILON + 1
  178.         ELSE
  179.           ILON  = npt / NBTHR
  180.           IDEB  = (IRES * (ILON+1)) + (ithr-IRES-1)* ILON + 1
  181.         ENDIF
  182.       ENDIF
  183.       IFIN  = IDEB + ILON - 1
  184.  
  185. c       print *, "PSI3D0 ith=",ithr,"/",NBTHR," ==>",IDEB,IFIN,"/",npt
  186.  
  187. *=======================================================================
  188. *      Boucle sur les noeuds ======================================
  189. c        do 490 ipt=1,npt
  190.        do 490 ipt=IDEB,IFIN
  191.  
  192.        isupfi=0
  193.        igard=0
  194.        igardf=0
  195.        id=ipt3.num(1,ipt)
  196. c        print *, "NOEUD",ipt,"/",npt," #",id
  197.  
  198. *      recup des coordonnees du noeud en question
  199.        xxp = xcoor((id-1)*idim1 + 1 )
  200.        yyp = xcoor((id-1)*idim1 + 2 )
  201.        zzp = xcoor((id-1)*idim1 + 3 )
  202.  
  203. *=====================================================================
  204. *           Boucle sur les elements de la fissure =========================
  205.             do 400 iseg=1,nseg
  206.  
  207.               dis = xgrand
  208.               iphi2 = 0
  209.  
  210. c             Recup des infos de la fissure (TRI3)
  211.               ia=MTRI3.KNODE(1,iseg)
  212.               ib=MTRI3.KNODE(2,iseg)
  213.               ic=MTRI3.KNODE(3,iseg)
  214. c               if(id.eq.idbug) print *,"TRI",iseg,"/",nseg," #",ia,ib,ic
  215.  
  216. *             recup des coordonnees du triangle
  217.               x11=MTRI3.XNODE(1,iseg)
  218.               x22=MTRI3.XNODE(2,iseg)
  219.               x33=MTRI3.XNODE(3,iseg)
  220.               y11=MTRI3.YNODE(1,iseg)
  221.               y22=MTRI3.YNODE(2,iseg)
  222.               y33=MTRI3.YNODE(3,iseg)
  223.               z11=MTRI3.ZNODE(1,iseg)
  224.               z22=MTRI3.ZNODE(2,iseg)
  225.               z33=MTRI3.ZNODE(3,iseg)
  226. *             recup de la normale au triangle + vecteurs des aretes + surface
  227.               vx=MTRI3.VNORM(1,iseg)
  228.               vy=MTRI3.VNORM(2,iseg)
  229.               vz=MTRI3.VNORM(3,iseg)
  230.               v12x=MTRI3.VARE12(1,iseg)
  231.               v12y=MTRI3.VARE12(2,iseg)
  232.               v12z=MTRI3.VARE12(3,iseg)
  233.               v12lo=MTRI3.VARE12(4,iseg)
  234.               v23x=MTRI3.VARE23(1,iseg)
  235.               v23y=MTRI3.VARE23(2,iseg)
  236.               v23z=MTRI3.VARE23(3,iseg)
  237.               v31x=MTRI3.VARE31(1,iseg)
  238.               v31y=MTRI3.VARE31(2,iseg)
  239.               v31z=MTRI3.VARE31(3,iseg)
  240.               sur =MTRI3.SURFAC(iseg)
  241.               xsur=1.D-8*sur
  242.  
  243. *------------ cas ou le noeud structure id = l un des noeuds des tri3 de la fissure
  244.               if (ia.eq.id.or.ib.eq.id.or.ic.eq.id) then
  245.  
  246.                 dis=0.D0
  247.                 iphi2 = 1
  248.                 phicor=.false.
  249.  
  250. *------------ "autres" cas
  251.               else
  252.  
  253. *bp2013         on saute ce point si trop loin (critere)
  254.                 if (icle.eq.1) then
  255.                   if (x11.ge.x22) then
  256.                     xxmi=min(x22,x33) - xcrit
  257.                     xxma=max(x11,x33) + xcrit
  258.                   else
  259.                     xxmi=min(x11,x33) - xcrit
  260.                     xxma=max(x22,x33) + xcrit
  261.                   endif
  262.                   if (y11.ge.y22) then
  263.                     yymi=min(y22,y33) - xcrit
  264.                     yyma=max(y11,y33) + xcrit
  265.                   else
  266.                     yymi=min(y11,y33) - xcrit
  267.                     yyma=max(y22,y33) + xcrit
  268.                   endif
  269.                   if (z11.ge.z22) then
  270.                     zzmi=min(z22,z33) - xcrit
  271.                     zzma=max(z11,z33) + xcrit
  272.                   else
  273.                     zzmi=min(z11,z33) - xcrit
  274.                     zzma=max(z22,z33) + xcrit
  275.                   endif
  276.                   if(xxp.lt.xxmi) goto 491
  277.                   if(xxp.gt.xxma) goto 491
  278.                   if(yyp.lt.yymi) goto 491
  279.                   if(yyp.gt.yyma) goto 491
  280.                   if(zzp.lt.zzmi) goto 491
  281.                   if(zzp.gt.zzma) goto 491
  282.                 endif
  283.  
  284.                 d1xp=xxp-x11
  285.                 d1yp=yyp-y11
  286.                 d1zp=zzp-z11
  287.                 d2xp=xxp-x22
  288.                 d2yp=yyp-y22
  289.                 d2zp=zzp-z22
  290.                 d3xp=xxp-x33
  291.                 d3yp=yyp-y33
  292.                 d3zp=zzp-z33
  293. * calcul de la distance au triangle
  294.                 d1= sqrt(d1xp**2 + d1yp**2 + d1zp**2)
  295.                 d2= sqrt(d2xp**2 + d2yp**2 + d2zp**2)
  296.                 d3= sqrt(d3xp**2 + d3yp**2 + d3zp**2)
  297.                 if (d1.lt.d2) then
  298.                     rr= min(d1,d3)
  299.                 else
  300.                     rr= min(d2,d3)
  301.                 endif
  302.  
  303. * calcul de la distance au plan du triangle
  304.                 xdd = vx*d1xp + vy*d1yp + vz*d1zp
  305. * determination du pt projetté
  306.                 xx=xxp-xdd*vx
  307.                 yy=yyp-xdd*vy
  308.                 zz=zzp-xdd*vz
  309.  
  310. * ce pt est-il à l'interieur ou sur un coté du triangle?
  311.                 v1xx=xx-x11
  312.                 v1yy=yy-y11
  313.                 v1zz=zz-z11
  314.                 v2xx=xx-x22
  315.                 v2yy=yy-y22
  316.                 v2zz=zz-z22
  317.                 v3xx=xx-x33
  318.                 v3yy=yy-y33
  319.                 v3zz=zz-z33
  320.                 pv12x=v1yy*v2zz-v1zz*v2yy
  321.                 pv12y=v1zz*v2xx-v1xx*v2zz
  322.                 pv12z=v1xx*v2yy-v1yy*v2xx
  323.                 pv23x=v2yy*v3zz-v2zz*v3yy
  324.                 pv23y=v2zz*v3xx-v2xx*v3zz
  325.                 pv23z=v2xx*v3yy-v2yy*v3xx
  326.                 pv31x=v3yy*v1zz-v3zz*v1yy
  327.                 pv31y=v3zz*v1xx-v3xx*v1zz
  328.                 pv31z=v3xx*v1yy-v3yy*v1xx
  329.                 ps12v=pv12x*vx+pv12y*vy+pv12z*vz
  330.                 ps23v=pv23x*vx+pv23y*vy+pv23z*vz
  331.                 ps31v=pv31x*vx+pv31y*vy+pv31z*vz
  332.  
  333. c            if(id.eq.idbug) then
  334. c            write(6,*) '-----------------------------------------------'
  335. c            write(6,*) 'noeud:',id,' (',ipt,') tri3',iseg,':',ia,ib,ic
  336. c            write(6,*) '   (',xx,')    (',x11,')   (',x22,')   (',x33,')'
  337. c            write(6,*) ' x=(',yy,')  1:(',y11,') 2:(',y22,') 3:(',y33,')'
  338. c            write(6,*) '   (',zz,')    (',z11,')   (',z22,')   (',z33,')'
  339. c            write(6,*)'rr=d1,2,3?',(rr.eq.d1),(rr.eq.d2),(rr.eq.d3),rr
  340. c            write(6,*) 'xdd=',xdd
  341. c            write(6,*) 'ps12v,23,31=',ps12v,ps23v,ps31v
  342. c            endif
  343.  
  344. * ce pt est-il à l'interieur ou sur un coté du triangle?
  345.                 if(ps12v.ge.0..and.ps23v.ge.0..and.ps31v.ge.0.) then
  346. * -> on est a l'interieur
  347.                     dis = xdd
  348.                     iphi2 = 1
  349.                 elseif (abs(ps12v).lt.xsur) then
  350. * -> on est situé sur le segment 12 : est-on a l'interieur?
  351.                   psc112=v1xx*v12x+v1yy*v12y+v1zz*v12z
  352.                   psc212=v2xx*v12x+v2yy*v12y+v2zz*v12z
  353.                   if ( psc112*psc212 .lt.0.d0) then
  354.                     dis  =  xdd
  355.                     iphi2 = 1
  356.                   endif
  357.                 elseif(abs(ps23v).lt.xsur)then
  358. * -> on est situé sur le segment 23 : est-on a l'interieur?
  359.                   psc223=v2xx*v23x+v2yy*v23y+v2zz*v23z
  360.                   psc323=v3xx*v23x+v3yy*v23y+v3zz*v23z
  361.                   if ( psc223*psc323 .lt.0.d0) then
  362.                     dis  =  xdd
  363.                     iphi2 = 1
  364.                   endif
  365.                 elseif(abs(ps31v).lt.xsur)then
  366. * -> on est situé sur le segment 31 : est-on a l'interieur?
  367.                   psc331=v3xx*v31x+v3yy*v31y+v3zz*v31z
  368.                   psc131=v1xx*v31x+v1yy*v31y+v1zz*v31z
  369.                   if ( psc331*psc131 .lt.0.d0) then
  370.                     dis  =  xdd
  371.                     iphi2 = 1
  372.                   endif
  373.                 endif
  374.  
  375. * on est a l'exterieur : il faut calculer dis et verifier xdd
  376.                 if (dis.eq.xgrand) then
  377.                   iphi2 = 10
  378.  
  379. * si on est a l'ext du tri3 et en avant du front+contour on change v
  380.  
  381. *              -si on a affaire a un triangle du front (ou du contour) :
  382. *               est on en avant du front (ou du contour)? -> phicor
  383.                 phicor=.false.
  384.                 phi12=.false.
  385.                 phi23=.false.
  386.                 phi31=.false.
  387.                 if (ideux.ge.1) then
  388.                 if (ifis(iseg,1).ge.1) then
  389. *               -Cas ou le segment [1-2] = front ou contour
  390.                  if(ifis(iseg,1).ge.2.or.ifis(iseg,2).eq.1.or.
  391.      $              ifis(iseg,2).eq.3.or.ifis(iseg,2).eq.5) then
  392. c                     phi12=(ps12v.lt.0.d0)
  393.                     phi12=(ps12v.lt.(-xsur))
  394.                     if(phi12) then
  395.                     if(ifis(iseg,1).ge.2) then
  396.                       iphi2=iphi2+5
  397.                     else
  398.                       iphi2=iphi2+1
  399.                     endif
  400.                     endif
  401.                  endif
  402. *               -Cas ou le segment [2-3] =  contour
  403.                  if(ifis(iseg,2).eq.2.or.
  404.      $              ifis(iseg,2).eq.3.or.ifis(iseg,2).eq.6) then
  405. c                     phi23=(ps23v.lt.0.d0)
  406.                     phi23=(ps23v.lt.(-xsur))
  407.                     if(phi23) iphi2=iphi2+1
  408.                  endif
  409. *               -Cas ou le segment [3-1] =  contour
  410.                  if(ifis(iseg,2).eq.4.or.
  411.      &              ifis(iseg,2).eq.5.or.ifis(iseg,2).eq.6) then
  412. c                     phi31=(ps31v.lt.0.d0)
  413.                     phi31=(ps31v.lt.(-xsur))
  414.                     if(phi31) iphi2=iphi2+1
  415.                  endif
  416. *               -Union des test
  417.                  phicor=(ideux.ge.1).and.(phi12.or.phi23.or.phi31)
  418.                 endif
  419.                 endif
  420. *               si phicor, il faudra alors "corriger" phi
  421. *               (et prendre xdd au lieu de dis)
  422. c                 if(id.eq.idbug) print *,"phicor =",phicor
  423.  
  424. c                   psc1=0.d0
  425. c                   psc2=0.d0
  426. c                   psc3=0.d0
  427. c                     dis=sign(rr,xdd)
  428. * bp (05/03/2012) : on remplace ci dessus par formule + precise
  429. * bp,gg (2019-09-25) : correction pour triangle quelconque
  430. *                   produits scalaires determinants si on est proche 
  431. *                   d'un noeud ou d'un segment
  432.                     psc112=v1xx*v12x+v1yy*v12y+v1zz*v12z
  433.                     psc131=v1xx*v31x+v1yy*v31y+v1zz*v31z                   
  434.                     psc223=v2xx*v23x+v2yy*v23y+v2zz*v23z
  435.                     psc212=v2xx*v12x+v2yy*v12y+v2zz*v12z                   
  436.                     psc331=v3xx*v31x+v3yy*v31y+v3zz*v31z
  437.                     psc323=v3xx*v23x+v3yy*v23y+v3zz*v23z
  438. *                 --on teste si on est dans un angle du triangle--
  439.                     if (psc112.le.0.d0.and.psc131.ge.0d0) then
  440. *                       on est definitivement + proche du noeud 1
  441.                         dis = sign(rr,xdd)
  442.                     elseif (psc223.le.0.d0.and.psc212.ge.0d0) then
  443. *                       on est definitivement + proche du noeud 2
  444.                         dis = sign(rr,xdd)
  445.                     elseif (psc331.le.0.d0.and.psc323.ge.0d0) then
  446. *                       on est definitivement + proche du noeud 3
  447.                         dis = sign(rr,xdd)
  448. *                 --on teste quel segment est le plus proche--
  449.                     elseif(ps12v.le.0.d0.and.psc112.ge.0.d0
  450.      &                                  .and.psc212.le.0) then
  451. *                       on est du coté du segment 12
  452.                         dis = sign(sqrt(MAX(d1**2-psc112**2,0.D0)),xdd)                        
  453.                     elseif(ps23v.le.0.d0.and.psc223.ge.0.d0
  454.      &                                  .and.psc323.le.0) then
  455. *                       on est du coté du segment 23
  456.                         dis = sign(sqrt(MAX(d2**2-psc223**2,0.D0)),xdd)                        
  457.                     elseif(ps31v.le.0.d0.and.psc331.ge.0.d0
  458.      &                                  .and.psc131.le.0) then
  459. *                       on est du coté du segment 31
  460.                         dis = sign(sqrt(MAX(d3**2-psc331**2,0.D0)),xdd)                        
  461.                     else
  462. c                      on ne devrait pas arriver ici !
  463.                        print *,"psi3d0 : noeud",id," tri",ia,ib,ic
  464.                        print *,psc112,psc131,psc223,psc212,psc331,psc323
  465. c                        call erreur(5)
  466.                        SPARAL.IERROR(ithr)=5
  467.                        return
  468.                     endif
  469.  
  470.                 endif
  471. * fin du cas ou on est a l'exterieur
  472.  
  473.               endif
  474. *------------ fin des "autres" cas (noeud id different de ceux du tri3)
  475.  
  476. *            -on a calculé (xdd et) dis: s'agit il de phi?
  477. c             on prend xdd plutot que dis pour definir Phi dans le cas ou :
  478. c -on a l option DEUX ou TROI
  479. c -le tri3 utilisé (le + proche) appartient au front ou au contour
  480. c -le point se projette en avant du front ou du contour
  481. *             iphi2 qu on va eventuellement ranger dans isup(ipt) vaut:
  482. *              1 si on se projette a l int du triangle
  483. *             10 si on se projette a l ext du triangle
  484. *             11 si en + on sort du contour
  485. *             15 si en + on est en avant du front
  486. *             16 si en + on est en avant du front et on sort du contour
  487.  
  488. c            if(id.eq.idbug) then
  489. c            write(6,*) 'dis,xdd',dis,xdd
  490. c            write(6,*) 'phicor=',phicor
  491. c            write(6,*) 'phi12,phi23,phi31=',phi12,phi23,phi31
  492. c            endif
  493.  
  494. *            -c'est la premiere fois que l'on voit ce point
  495.               if (isup(ipt).eq.0) then
  496.                   if(phicor) then
  497.                     mpova2.vpocha(ipt,1)=xdd
  498.                   else
  499.                     mpova2.vpocha(ipt,1)=dis
  500.                   endif
  501.                   xdmin =dis
  502.                   xdmin1=xdd
  503.                   isup(ipt)=iphi2
  504.                   iphi2=100
  505.                   igard=igard+1
  506.               else
  507. *             --on est (franchement) + proche que les precedents essais
  508.                 if (abs(dis).lt.0.999999*abs(xdmin)) then
  509.                   if(phicor) then
  510.                     mpova2.vpocha(ipt,1)=xdd
  511.                   else
  512.                     mpova2.vpocha(ipt,1)=dis
  513.                   endif
  514.                    xdmin =dis
  515.                    xdmin1=xdd
  516.                    isup(ipt)=iphi2
  517.                    iphi2=100
  518. *             --on est (quasiment) aussi proche que le meilleur essai
  519.                 elseif (abs(dis).lt.1.000001*abs(xdmin)) then
  520. *                -on a 2 tri3 du front+contour a egalite
  521. *                 et l'un des deux a trouvé le point en avant du front :
  522.                   if ((iseg.le.imil5).and.
  523.      $                (phicor.or.isup(ipt).gt.10)) then
  524. c                   strategie : xdd= max(xdd1 xdd2)
  525.                     if (abs(xdd).gt.abs(xdmin1)) then
  526.                       mpova2.vpocha(ipt,1)=xdd
  527.                       xdmin =dis
  528.                       xdmin1=xdd
  529. c                       isup(ipt)=iphi2
  530.                       if(iphi2.gt.isup(ipt)) isup(ipt)=iphi2
  531.                       iphi2=100
  532. c                   cas ou on avait pas vu qu il fallait prendre xdd
  533.                     elseif(isup(ipt).le.10) then
  534. c                   equivalent à :
  535. c                     elseif(mpova2.vpocha(ipt,1).ne.xdmin1(ipt)) then
  536.                       mpova2.vpocha(ipt,1)=xdmin1
  537.                       isup(ipt)=iphi2
  538.                       iphi2=100
  539.                     endif
  540.  
  541. *                -on avait 1 tri3 du front+contour et ce tri3 est a egalite
  542. *                 et le 1er avait trouvé le point en avant du front :
  543. *                 on laisse le xdd calculé depuis le 1er tri3
  544.                   elseif ((iseg.gt.imil5).and.
  545.      $                    (isup(ipt).gt.10)) then
  546.                     goto 491
  547. *                -quasi-egalite entre 2 concurrents et possibles signes opposés :
  548.                   else
  549. c bp: nouvelle strategie: on moyenne les normales*(y-xp) (soit xdd)
  550.                     xdmin1=xdmin1+xdd
  551.                     mpova2.vpocha(ipt,1)=sign(dis,xdmin1)
  552.                   endif
  553.                 endif
  554. *             --fin de la distinction + proche / aussi proche
  555.               endif
  556. *            -fin de la distinction 1ere fois ou pas
  557.  
  558. c              if(id.eq.idbug) then
  559. c              write(6,*)'iphi2,xdmin,xdmin1,phi',iphi2,
  560. c      $             xdmin,xdmin1,mpova2.vpocha(ipt,1)
  561. c              endif
  562. c              if(iimpi.ge.2) then
  563. c                write(ioimp,*) 'point',id,ipt,'dis,xdd,PHI,isup=',
  564. c      $         dis,xdd,mpova2.vpocha(ipt,1),isup(ipt)
  565. c              endif
  566.  
  567. *            point d'arrivee si hors-boite
  568.   491        continue
  569. *            le tri3 n'est pas un tri3 du front => on saute
  570. *TODO : test a verifier peut etre .eq. et .eq.-1 aussi a ne pas sauter?
  571.              if (ideux.lt.1)        goto 400
  572.              if (ifis(iseg,1).le.1) goto 400
  573.  
  574. c   ------   Calcul de la 2eme Level Set   -------
  575.              xdd1 = 0.D0
  576.  
  577. *            -cas ou le noeud structure id = l un des noeuds des tri3 de la fissure
  578.               if (ia.eq.id.or.ib.eq.id) then
  579.                    dis1=0.D0
  580.                    xd2=0.D0
  581.                    xd3=0.D0
  582.                    if(ideux.ge.2) xtau=xicpr(id)
  583. *            -autres cas
  584.               else
  585. c *              on recupere ce qu on a pas recupéré lors du calcul de phi
  586. c                if (ic.eq.id) then
  587. c                   xxp = xcoor((id-1)*idim1 + 1 )
  588. c                   yyp = xcoor((id-1)*idim1 + 2 )
  589. c                   zzp = xcoor((id-1)*idim1 + 3 )
  590. c *               + coordonnees du triangle+...                
  591. c                endif
  592. cbp, desormais inutile car deplace + haut
  593.  
  594. * calcul du vecteur vn qui est le sens de propagation
  595.                vnx=v12y*vz-v12z*vy
  596.                vny=v12z*vx-v12x*vz
  597.                vnz=v12x*vy-v12y*vx
  598. * il faut orienter vnx pour que le vecteur  soit  dans le sens
  599. *  de la propagation
  600.                qq= vnx*(x33-x11)+vny*(y33-y11)+vnz*(z33-z11)
  601.                if (qq.gt.0.D0) then
  602.                  vnx=-vnx
  603.                  vny=-vny
  604.                  vnz=-vnz
  605.                endif
  606. * normalement vn est normé on le verifie
  607.                qq= sqrt(vnx**2 + vny**2 + vnz**2)
  608.                if  (abs(qq-1.d0).gt.1.d-2) then
  609.                  PRINT *,'PSIP3D0: pb dans le calcul de vn',ithr
  610. c                call erreur(26)
  611.                  SPARAL.IERROR(ithr)=26
  612.                  return
  613.                elseif(abs(qq-1.d0).gt.1.d-6) then
  614.                  vnx=vnx/qq
  615.                  vny=vny/qq
  616.                  vnz=vnz/qq
  617.                endif
  618. * on calcule la distance au plan defini par segment 12 et vecteur vx,vy,vz
  619. * et le point (xx yy zz ) projeté du point sur ce plan
  620.                xdd1=(xxp-x11)*vnx+(yyp-y11)*vny+(zzp-z11)*vnz
  621.                xx=xxp-xdd1*vnx
  622.                yy=yyp-xdd1*vny
  623.                zz=zzp-xdd1*vnz
  624.                psc1=(xx-x11)*v12x+(yy-y11)*v12y+(zz-z11)*v12z
  625.                psc2=(xx-x22)*v12x+(yy-y22)*v12y+(zz-z22)*v12z
  626.                if ((psc1*psc2).le.0.d0) then
  627. C le point se projete sur 12
  628.                    dis1=xdd1
  629.                    pp1x=x11+(psc1*v12x)
  630.                    pp1y=y11+(psc1*v12y)
  631.                    pp1z=z11+(psc1*v12z)
  632.                    xd2 = 0.D0
  633.                    xd3 = sqrt((xxp-pp1x)**2+(yyp-pp1y)**2+(zzp-pp1z)**2)
  634.                    if (ideux.ge.2) then
  635.                      xtau1=xicpr(ia)
  636.                      xtau2=xicpr(ib)
  637.                      stau = psc1 / v12lo
  638.                      xtau=stau*xtau2+(1.d0-stau)*xtau1
  639.                    endif
  640.                else
  641.                  if (psc1.ge.0.D0) then
  642.                    pp2x=xx+((v12lo-psc1)*v12x)
  643.                    pp2y=yy+((v12lo-psc1)*v12y)
  644.                    pp2z=zz+((v12lo-psc1)*v12z)
  645.                    xd2 = abs(v12lo-psc1)
  646.                    xd3 = sqrt((xxp-x22)**2+(yyp-y22)**2+(zzp-z22)**2)
  647. c                    if(ideux.ge.2) xtau=xicpr(ib)
  648.                    if(ideux.ge.2) then
  649.                      if(ib.eq.i2firs.or.ib.eq.i2last) then
  650. c                      write(6,*) 'le noeud',id,'detecte comme apres',ib
  651.                        xtau1=xicpr(ia)
  652.                        xtau2=xicpr(ib)
  653.                        stau = psc1 / v12lo
  654.                        xtau=stau*xtau2+(1.d0-stau)*xtau1
  655.                      else
  656.                        xtau=xicpr(ib)
  657.                      endif
  658.                    endif
  659.                  else
  660.                    pp2x=xx-(psc1*v12x)
  661.                    pp2y=yy-(psc1*v12y)
  662.                    pp2z=zz-(psc1*v12z)
  663.                    xd2 = abs(psc1)
  664.                    xd3 = sqrt((xxp-x11)**2+(yyp-y11)**2+(zzp-z11)**2)
  665. c                    if(ideux.ge.2) xtau=xicpr(ia)
  666.                    if(ideux.ge.2) then
  667.                      if(ia.eq.i2firs.or.ia.eq.i2last) then
  668. c                      write(6,*) 'le noeud',id,'detecte comme avant',ia
  669.                        xtau1=xicpr(ia)
  670.                        xtau2=xicpr(ib)
  671.                        stau = psc1 / v12lo
  672.                        xtau=stau*xtau2+(1.d0-stau)*xtau1
  673.                      else
  674.                        xtau=xicpr(ia)
  675.                      endif
  676.                    endif
  677.                  endif
  678.                  dis1=sqrt((xxp-pp2x)**2+(yyp-pp2y)**2+(zzp-pp2z)**2)
  679.                  dis1=sign(dis1,xdd1)
  680.                endif
  681.  
  682.               endif
  683. *            -fin de la distinction si noeud du tri3 ou pas
  684.  
  685. c                if(id.eq.idbug) print *,"isupfi =",isupfi
  686. *              1ere fois qu on traite ce point?
  687.                if (isupfi.eq.0) then
  688. *                 isupfi(ipt)=1
  689. *                on utilise isupfi pour definir si la valeur inscrite
  690. *                provient du 1er ou dernier seg du front (=2) ou autre(=1)?
  691.                  isupfi=ifis(iseg,1)
  692.                  igardf=igardf+1
  693.                  xdmin2=xd2
  694.                  xdmin3=xd3
  695. cbp2013                 mpova1.vpocha(ipt,1)=dis1
  696.                  if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau
  697.                else
  698. *                 if(abs(dis1).lt.abs(mpova1.vpocha(ipt,1)))then
  699. *                on a trouvé un segment + proche ou = du noeud considéré
  700. *                 if(abs(dis1).le.abs(mpova1.vpocha(ipt,1)))then
  701. **                 if(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))-1.D-8)then
  702. *                 if(xd2.le.xdmin2(ipt)) then
  703.                  if (xd3.lt.(abs(xdmin3)-1.D-8)) then
  704. cbp2013                     mpova1.vpocha(ipt,1)=dis1
  705.                      xdmin2=xd2
  706.                      xdmin3=xd3
  707.                      isupfi=ifis(iseg,1)
  708.                      if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau
  709.                  endif
  710. *                rattrapage dans le cas d'erreur d'arrondi sur dis1
  711. *                 elseif(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))+1.D-8)
  712. *     $           then
  713. *                 if(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))+1.D-8)then
  714.                  if (xd3.lt.(abs(xdmin3)+1.D-8))then
  715. *                  en cas d'égalité,on n'utilise la distance dis1 que si:
  716. *                  - le pt projetté appartient au segment (xd2=0)
  717. *              ou  - le point projeté est + proche du segment
  718.                    if (xd2.lt.xdmin2) then
  719. cbp2013                      mpova1.vpocha(ipt,1)=dis1
  720.                      xdmin2=xd2
  721.                      xdmin3=xd3
  722.                      isupfi=ifis(iseg,1)
  723.                      if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau
  724.                    endif
  725.                  endif
  726.                endif
  727.  
  728.   400       continue
  729. *           Fin de Boucle sur les elements de la fissure ===============
  730. *=======================================================================
  731.  
  732. *------   Petite correction de Psi
  733.           if (ideux.ge.1) then
  734.               if (isup(ipt).eq.0) goto 500
  735. c * pour le cas de fissure debouchante
  736. c * dans le cas ou c'est un segment debouchant qui a servi au calcul de psi,
  737. c * il faut prendre la distance xdd2 (et pas dis) qu'on retrouve avec :
  738. c               if ((isupfi(ipt)).eq.2) then
  739. c                  dis1 = mpova1.vpocha(ipt,1)
  740. c                  xdd1 = sqrt(MAX(dis1**2 - xdmin2(ipt)**2,0.D0))
  741. c                  mpova1.vpocha(ipt,1) = sign(xdd1,dis1)
  742. c               endif
  743. c bp : essai du 16/03/2012 : psi = distance au front**2-phi**2
  744.               xdd1 = (xdmin3**2) - (mpova2.vpocha(ipt,1)**2)
  745.               xdd1 = sqrt(max(xdd1,0.d0))
  746. c            if(ipt3.num(1,ipt).eq.idbug) then
  747. c            write(6,*)'ipt,isup,psi_avant,phi,xdd1,= ',ipt,isup(ipt),
  748. c      $     xdmin3,mpova1.vpocha(ipt,1),mpova2.vpocha(ipt,1),xdd1
  749. c            endif
  750. c               if (isup(ipt).gt.0) then
  751.               if (isup(ipt).ge.15) then
  752. *on est a l exterieur => psi >0
  753.                 mpova1.vpocha(ipt,1) = xdd1
  754.               else
  755.                 mpova1.vpocha(ipt,1) = -1.d0*xdd1
  756.               endif
  757.   500         continue
  758.           endif
  759. *------   Fin de la Petite correction de Psi
  760.  
  761.   490   continue
  762. *       Fin de Boucle sur les noeuds ===================================
  763. *=======================================================================
  764.  
  765.       END
  766.  
  767.  
  768.  
  769.  
  770.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales