Télécharger psi3d0.eso

Retour à la liste

Numérotation des lignes :

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

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