Télécharger psi3d0.eso

Retour à la liste

Numérotation des lignes :

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

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