Télécharger psip3d.eso

Retour à la liste

Numérotation des lignes :

  1. C PSIP3D SOURCE BP208322 14/01/29 21:15:14 7921
  2. c
  3. SUBROUTINE PSIP3D(ideux,melmai,melfis,melfro,xcrit,
  4. & msoup1,msoup2,msoup3)
  5.  
  6. ***********************************************************************
  7. c OPERATEUR : PSIP
  8.  
  9. c APPELE PAR PSIPHI DANS LE CAS 3D
  10. c
  11. c CREATION : BP, 14/03/2012
  12. c MODIF : BP, 19/12/2013 plus d'appel au zonage
  13. c
  14. ***********************************************************************
  15. IMPLICIT REAL*8(A-H,O-Z)
  16. IMPLICIT INTEGER(I-N)
  17. -INC CCOPTIO
  18. -INC SMCOORD
  19. -INC SMELEME
  20. -INC SMCHPOI
  21. -INC CCREEL
  22. segment icpr(xcoor(/1)/(idim+1))
  23. segment xicpr(xcoor(/1)/(idim+1))
  24. segment isup(npt)
  25. segment xdis(npt)
  26. segment xdmin(npt)
  27. segment xdmin1(npt)
  28. segment rrlim(npt)
  29. segment xdmin2(npt)
  30. segment xdmin3(npt)
  31. segment isupfi(npt)
  32. segment sfis
  33. integer ifis(nbelfi,2)
  34. endsegment
  35. DATA EPSI/1.D-5/
  36. logical phicor,phi12,phi23,phi31
  37. integer ooperm(7,2)
  38. DATA ooperm/ 2, 4, 6, 1, 3, 5, 7,
  39. $ 4, 1, 5, 2, 6, 3, 7/
  40.  
  41.  
  42. ***********************************************************************
  43. * Initialisation, Recup, Activation
  44. ***********************************************************************
  45. if(iimpi.ge.1) then
  46. call gibtem(xkt)
  47. write(ioimp,*)'===== ENTREE DANS PSIP3D ===== | xkt=',xkt
  48. endif
  49.  
  50. idim1=idim+1
  51. idbug = 0
  52.  
  53. c des chpoints (les msoup sont deja actifs)
  54. if(msoup2.gt.0) mpova2=msoup2.ipoval
  55. if(msoup1.gt.0) mpova1=msoup1.ipoval
  56. if(msoup3.gt.0) mpova3=msoup3.ipoval
  57. c activation des maillages / desactivation
  58. c melmai : fait dans psiphi / ligne 1325 ou 1329
  59. c melfis : fait ligne 64 / ligne 1274
  60. c melfro : fait (si ideux>=1) ligne 78 / ligne 356
  61. ipt3=melmai
  62. c segact,ipt3
  63. npt=ipt3.num(/2)
  64.  
  65. c variables parfois testées
  66. imil5=0
  67.  
  68.  
  69. ***********************************************************************
  70. * Preparation du maillage de la fissure
  71. ***********************************************************************
  72. meleme=melfis
  73. segact meleme
  74. * transfo du front en TRI3
  75. if(itypel.ne.4) call change(melfis,4)
  76. meleme=melfis
  77. segact meleme
  78. c rem: inutile car change laisse melfis actif
  79. nseg=num(/2)
  80.  
  81. *------ Cas ou l'on calcule phi ET psi --------------------------------
  82. if (ideux.ge.1) then
  83. * on rearrange le melfis de facon à prevenir dans ifis(iel,1) quels
  84. * element touchent le front de fissure. Ceux ci sont ordonnées de facon
  85. * que le front de fissure soit du noeud 1 vers le noeud 2 de l'element
  86. ipt2=melfro
  87. segact ipt2
  88. segini icpr
  89. nbelem2 = ipt2.num(/2)
  90. nbnode2 = ipt2.num(/1)
  91. * petite numerotation locale
  92. * icpr(noeud)>0 si le noeud appartient au front de fissure
  93. * =numero du 1er element de melfro ou on a vu ce noeud
  94. * + on en profite pour calculer l'abscisse curviligne du front si ideux.eq.2
  95. * icpr(noeud)<0 si le noeud appartient au reste du contour de la fissure
  96. * =-numero du 1er element de melcon ou on a vu ce noeud
  97. if (ideux.ge.2) then
  98. segini,xicpr
  99. xtau=0.d0
  100. ic=ipt2.num(1,1)
  101. x11= xcoor((ic-1)*idim1 +1)
  102. y11= xcoor((ic-1)*idim1 +2)
  103. z11= xcoor((ic-1)*idim1 +3)
  104. do ia=1,nbelem2
  105. do ib=1,nbnode2
  106. ic=ipt2.num(ib,ia)
  107. if (icpr(ic).eq.0) then
  108. icpr(ic)=ia
  109. x22= xcoor((ic-1)*idim1 +1)
  110. y22= xcoor((ic-1)*idim1 +2)
  111. z22= xcoor((ic-1)*idim1 +3)
  112. dxtau= sqrt( (x22-x11)*(x22-x11)+(y22-y11)*(y22-y11)
  113. $ +(z22-z11)*(z22-z11) )
  114. xtau=xtau+dxtau
  115. xicpr(ic)=xtau
  116. x11=x22
  117. y11=y22
  118. z11=z22
  119. c if(idebug.ge.1)
  120. c $ write(6,*) 'ia,ib,ic,xicpr(ic)=',ia,ib,ic,xicpr(ic)
  121. endif
  122. enddo
  123. enddo
  124. else
  125. do ia=1,nbelem2
  126. do ib=1,nbnode2
  127. ic=ipt2.num(ib,ia)
  128. if (icpr(ic).eq.0) then
  129. icpr(ic)=ia
  130. endif
  131. enddo
  132. enddo
  133. endif
  134. * bp (23/02/2012) : calcul du contour et fin de remplissage de icpr
  135. CALL ECROBJ('MAILLAGE',melfis)
  136. CALL PRCONT
  137. CALL LIROBJ('MAILLAGE',melcon,1,IRETOU)
  138. IF (IERR.NE.0) RETURN
  139. meleme = melfis
  140. IPT6 = melcon
  141. SEGACT,meleme,IPT6
  142. nbelCON = IPT6.num(/2)
  143. nbnoCON = IPT6.num(/1)
  144. do ia=1,nbelCON
  145. do ib=1,nbnoCON
  146. ic=IPT6.num(ib,ia)
  147. if(icpr(ic).eq.0) icpr(ic)=-ia
  148. enddo
  149. enddo
  150. * nombre de segments appartenant au contour sans le front
  151. nbelem6 = nbelCON - nbelem2
  152. * remplissage effectif de ifis contenu dans le segment sfis
  153. i2firs = ipt2.num(1,1)
  154. i2last = ipt2.num(nbnode2,nbelem2)
  155. ivui2f = 0
  156. ivui2l = 0
  157. c write(6,*) 'i2firs,i2last=',i2firs,i2last
  158. nbelfi=nseg
  159. nbelem=nseg
  160. nbnn=3
  161. nbsous=0
  162. nbref=0
  163. c segini,ipt5=meleme
  164. * bp 20/12/2011 : on propose de remplir ipt5 avec :
  165. c + des elements du front (et du contour) d abord,
  166. c + des elements avec 1 point appartenant au front (et au contour) ensuite,
  167. c + les autres elements enfin.
  168. c bp 06/03/2012 : on propose de remplir ifis avec :
  169. c ifis(iseg,1) =2 si le segment [1-2] du iseg^ieme triangle appartient
  170. c au front de fissure, =3 si en plus c'est l extrémité du front,
  171. c =1 si l un des segments appartient au contour
  172. c =-1 si l'un des noeuds appartient au front
  173. c ifis(iseg,2) =+1+2+4 si les segments [1-2] [2-3] [3-1] du iseg^ieme triangle appartiennent au contour
  174. segini,ipt5
  175. ideb5 = 0
  176. imil5 = nbelem2
  177. ifin5 = nseg+1
  178. segini sfis
  179.  
  180. *--------boucle sur les elements de la fissure ----------------------
  181. do iel=1,nseg
  182.  
  183. * io est incrémenté de +1,+2,+4 si node 1,2,3 appartient au front
  184. * ioo est incrémenté de +1,+2,+4 si node 1,2,3 appartient au contour
  185. io=0
  186. ioo=0
  187. * ---boucle sur les noeuds des elements de la fissure ----------
  188. do ib=1,num(/1)
  189. ia=ib
  190. if(ib.eq.3) ia = 4
  191. ic = num(ib,iel)
  192. if(icpr(ic).gt.0) io=io+ia
  193. if(icpr(ic).lt.0) ioo=ioo+ia
  194. * cas particulier d un segment avec 1 point appartenant au contour et un autre
  195. if(ic.eq.i2firs.or.ic.eq.i2last) ioo=ioo+ia
  196. enddo
  197. * ---fin de boucle sur les noeuds des elements de la fissure ---
  198.  
  199. * remplissage de ipt5(i5,...) et ifis(i5,1)
  200. i5=0
  201. iperm=0
  202. * ---cas ou 1 segment appartient au front
  203. if (io.eq.3.or.io.gt.4) then
  204. ideb5=ideb5+1
  205. i5=ideb5
  206. * (eventuelle) permutation des noeuds pour avoir segment 1-2 = front
  207. if (io.eq.3) then
  208. ipt5.num(1,i5)=num(1,iel)
  209. ipt5.num(2,i5)=num(2,iel)
  210. ipt5.num(3,i5)=num(3,iel)
  211. elseif(io.eq.5) then
  212. ipt5.num(1,i5)=num(3,iel)
  213. ipt5.num(2,i5)=num(1,iel)
  214. ipt5.num(3,i5)=num(2,iel)
  215. iperm=1
  216. elseif(io.eq.6) then
  217. ipt5.num(1,i5)=num(2,iel)
  218. ipt5.num(2,i5)=num(3,iel)
  219. ipt5.num(3,i5)=num(1,iel)
  220. iperm=2
  221. elseif(io.eq.7) then
  222. write(ioimp,*) 'ERREUR: '
  223. $ ,'Plusieurs cotes d un meme triangle du maillage de '
  224. $ ,'la fissure converti en TRI3 appartiennent au front!'
  225. return
  226. endif
  227. * ifis( ,1)=3 si c'est une extremité
  228. if ((ipt5.num(1,i5)).eq.i2firs) then
  229. ifis(i5,1)=3
  230. ivui2f = ivui2f+1
  231. elseif((ipt5.num(2,i5)).eq.i2last) then
  232. ifis(i5,1)=3
  233. ivui2l = ivui2l+1
  234. else
  235. * ifis( ,1)=2 si c'est un tri3 du front avec des voisins
  236. ifis(i5,1)=2
  237. endif
  238. * ---cas ou 1 (ou 2) segment(s) appartient au contour
  239. elseif(ioo.eq.3.or.ioo.gt.4) then
  240. imil5=imil5+1
  241. i5=imil5
  242. * pas de permutation des noeuds
  243. ipt5.num(1,i5)=num(1,iel)
  244. ipt5.num(2,i5)=num(2,iel)
  245. ipt5.num(3,i5)=num(3,iel)
  246. * ifis( ,1)=1 si un segment appartient au contour
  247. ifis(i5,1)=1
  248. * ---cas ou seulement 1 point appartient au front
  249. elseif(io.eq.1.or.io.eq.2.or.io.eq.4) then
  250. ifin5=ifin5-1
  251. i5=ifin5
  252. * pas de permutation des noeuds
  253. ipt5.num(1,i5)=num(1,iel)
  254. ipt5.num(2,i5)=num(2,iel)
  255. ipt5.num(3,i5)=num(3,iel)
  256. * ifis( ,1)=-1 si 1 seul point du tri3 appartient au front
  257. ifis(i5,1)=-1
  258. * ---autres cas
  259. else
  260. ifin5=ifin5-1
  261. i5=ifin5
  262. * pas de permutation des noeuds
  263. ipt5.num(1,i5)=num(1,iel)
  264. ipt5.num(2,i5)=num(2,iel)
  265. ipt5.num(3,i5)=num(3,iel)
  266. endif
  267.  
  268. * remplissage de ifis(i5,2)
  269. * ---uniquement si 1 (ou 2) segment(s) appartient au contour
  270. if(ioo.eq.3.or.ioo.gt.4) then
  271. * attention a l eventuelle permutation effectuee
  272. if(iperm.ne.0) ioo=ooperm(ioo,iperm)
  273. * -ioo=3 => noeuds 1-2 = segment 1 appartient au contour
  274. if (ioo.eq.3) then
  275. ifis(i5,2)=1
  276. * -ioo=5 => noeuds 3-1 = segment 3 appartient au contour
  277. elseif(ioo.eq.5) then
  278. ifis(i5,2)=4
  279. * -ioo=6 => noeuds 2-3 = segment 2 appartient au contour
  280. elseif(ioo.eq.6) then
  281. ifis(i5,2)=2
  282. * -cas ou 3 points appartiennent au contour => travail en +
  283. elseif(ioo.eq.7) then
  284. c write(6,*)'tri3 avec 2 segments dans le contour!',iel,i5
  285. * rem: si les 3 points appartiennent au contour, on n'a pas de segments
  286. * appartenant au front
  287. ii1 = ipt5.num(1,i5)
  288. ii2 = ipt5.num(2,i5)
  289. ii3 = ipt5.num(3,i5)
  290. c il faut trouver quels sont les 2 segments appartenant au contour !
  291. ia1 = -icpr(ii1)
  292. ia2 = -icpr(ii2)
  293. ia3 = -icpr(ii3)
  294. ia1min = min(ia1,ia2)
  295. ia1min = min(ia1min,ia3)
  296. ia1seg = 0
  297. do 11 ia1=ia1min,nbelCON
  298. ifis1 = 0
  299. do ib=1,nbnoCON
  300. if(ii1.eq.IPT6.num(ib,ia1)) ifis1=ifis1+1
  301. if(ii2.eq.IPT6.num(ib,ia1)) ifis1=ifis1+2
  302. if(ii3.eq.IPT6.num(ib,ia1)) ifis1=ifis1+4
  303. enddo
  304. if(ifis1.le.2) goto 11
  305. if(ifis1.eq.4) goto 11
  306. if(ifis1.eq.3) ifis(i5,2)=ifis(i5,2)+1
  307. if(ifis1.eq.6) ifis(i5,2)=ifis(i5,2)+2
  308. if(ifis1.eq.5) ifis(i5,2)=ifis(i5,2)+4
  309. if(ifis1.ge.7) then
  310. write(ioimp,*) 'tri3 avec 3 segments dans le contour !'
  311. call erreur(26)
  312. endif
  313. if(ifis1.eq.3.or.ifis1.ge.5) ia1seg=ia1seg+1
  314. if(ia1seg.eq.2) goto 12
  315. 11 continue
  316. c on n a pas trouvé les 2 segments appartenant au contour
  317. write(ioimp,*)'on n a pas trouvé les 2 segments',
  318. $ ' du tri3 devant appartenant au contour !'
  319. call erreur(26)
  320. c on a trouvé les 2 segments appartenant au contour
  321. 12 continue
  322. if(ifis(i5,2).lt.3.or.ifis(i5,2).ge.7) then
  323. write(ioimp,*) 'impossible de detecter les 2 segments ',
  324. $ 'appartenant au contour avec les noeuds ',ii1,ii2,ii3
  325. write(ioimp,*) 'ifis(ideb5,2)=',ifis(ideb5,2)
  326. write(ioimp,*) 'IPT6 pour les segments :'
  327. write(ioimp,*) ' ia1',ia1,(IPT6.num(iou,ia1),iou=1,2)
  328. write(ioimp,*) ' ia2',ia2,(IPT6.num(iou,ia2),iou=1,2)
  329. write(ioimp,*) ' ia3',ia3,(IPT6.num(iou,ia3),iou=1,2)
  330. call erreur(26)
  331. return
  332. endif
  333. endif
  334. endif
  335.  
  336.  
  337.  
  338. enddo
  339. * fin de boucle sur les elements de la fissure -----------------
  340.  
  341.  
  342. c if ((ifin5-ideb5).ne.1) then
  343. if (ideb5.ne.nbelem2.or.(ifin5-imil5).ne.1) then
  344. write(ioimp,*) 'pb avec le rearrangement de elements de fissure'
  345. write(ioimp,*) ideb5,nbelem2,imil5,ifin5
  346. call erreur(26)
  347. return
  348. endif
  349. * si le front n est pas ferme, il faut avoir detecte ses 2 extremites
  350. * sinon il est probable que la ligne du front melfro ne soit pas
  351. * orientée comme le bord du maillage melfis
  352. if ((ivui2f.ne.1).and.(ivui2l.ne.1)) then
  353. write(ioimp,*) 'ERREUR: Les extremites du front (noeuds',
  354. $ i2firs,i2last,') ont été detectees ',ivui2f, 'et ',ivui2l
  355. $ ,' fois au lieu de 1 fois attendue'
  356. write(ioimp,*) 'Verifiez la coherence avec l orientation'
  357. $ ,'du maillage de la fissure fourni'
  358. $ ,' et INVErser le sens du front le cas échéant.'
  359. c write(6,*) 'front de fissure='
  360. c CALL ECMAIL(ipt2,JENT2)
  361. c write(6,*) 'surface de fissure='
  362. c CALL ECMAIL(ipt5,JENT5)
  363. call erreur(26)
  364. return
  365. endif
  366. c segdes,meleme,IPT6
  367. segdes,meleme,IPT6,ipt2
  368. meleme=ipt5
  369. melfis=meleme
  370. segsup icpr
  371.  
  372. endif
  373. *------ fin du Cas ou l'on calcule phi ET psi --------------------------
  374.  
  375. * sont actifs: melmai(=struct changée en poi1), ipt5=meleme=melfis,
  376. * sont inactifs: ipt2=melfro
  377.  
  378.  
  379.  
  380. ***********************************************************************
  381. * debut du travail suivant valeur de icle
  382.  
  383. if (xcrit.eq.0.d0) then
  384. icle=0
  385. xcrit=1.D10
  386. else
  387. icle=1
  388. * ajout d'un epsilon 1.D-5 de rattrapage
  389. xcrit = 1.00001D0 * xcrit
  390. endif
  391. if(iimpi.ge.1) then
  392. call gibtem(xkt)
  393. write(ioimp,*)'debut zonage si icle=1=',icle,xcrit,'| xkt=',xkt
  394. endif
  395.  
  396. melfis=meleme
  397. melpoi=melmai
  398.  
  399. ***********************************************************************
  400. * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT
  401. ***********************************************************************
  402.  
  403. * pour chaque segment de la fissure (melfis 2eme maillage lu)
  404. * on calcule le champ phi et psi des noeuds autour de l'element
  405. if(iimpi.ge.1) then
  406. call gibtem(xkt)
  407. write(ioimp,*)'debut travail | xkt=',xkt
  408. endif
  409. *
  410. meleme=melfis
  411. c segact,meleme
  412. ipt3=melpoi
  413. c segact,ipt3
  414. npt=ipt3.num(/2)
  415. * attention on initialise isup avec les valeurs inverses,il faudra
  416. * le retourner
  417. segini isup,isupfi,xdis,xdmin,xdmin1,rrlim
  418. if(ideux.ge.1) segini xdmin2,xdmin3
  419. igard=0
  420. igardf=0
  421.  
  422.  
  423. *=====================================================================
  424. * Boucle sur les elements de la fissure =========================
  425.  
  426. do 400 iseg=1,nseg
  427.  
  428. c coordonnes du tri3
  429. ia=num(1,iseg)
  430. x11= xcoor((ia-1)*idim1 +1)
  431. y11= xcoor((ia-1)*idim1 +2)
  432. z11= xcoor((ia-1)*idim1 +3)
  433. ib=num(2,iseg)
  434. x22= xcoor((ib-1)*idim1 +1)
  435. y22= xcoor((ib-1)*idim1 +2)
  436. z22= xcoor((ib-1)*idim1 +3)
  437. ic=num(3,iseg)
  438. x33= xcoor((ic-1)*idim1 +1)
  439. y33= xcoor((ic-1)*idim1 +2)
  440. z33= xcoor((ic-1)*idim1 +3)
  441. if (x11.ge.x22) then
  442. xxmi=min(x22,x33) - xcrit
  443. xxma=max(x11,x33) + xcrit
  444. else
  445. xxmi=min(x11,x33) - xcrit
  446. xxma=max(x22,x33) + xcrit
  447. endif
  448. if (y11.ge.y22) then
  449. yymi=min(y22,y33) - xcrit
  450. yyma=max(y11,y33) + xcrit
  451. else
  452. yymi=min(y11,y33) - xcrit
  453. yyma=max(y22,y33) + xcrit
  454. endif
  455. if (z11.ge.z22) then
  456. zzmi=min(z22,z33) - xcrit
  457. zzma=max(z11,z33) + xcrit
  458. else
  459. zzmi=min(z11,z33) - xcrit
  460. zzma=max(z22,z33) + xcrit
  461. endif
  462.  
  463. c vecteurs des aretes
  464. v12x=x22-x11
  465. v12y=y22-y11
  466. v12z=z22-z11
  467. v23x=x33-x22
  468. v23y=y33-y22
  469. v23z=z33-z22
  470. v31x=x11-x33
  471. v31y=y11-y33
  472. v31z=z11-z33
  473.  
  474. c vecteur normal
  475. vx=v12y*(z33-z11)-v12z*(y33-y11)
  476. vy=v12z*(x33-x11)-v12x*(z33-z11)
  477. vz=v12x*(y33-y11)-v12y*(x33-x11)
  478. bb=vx*vx + vy*vy+ vz*vz
  479. aa= sqrt(bb)
  480. if (aa.le.XPETIT) then
  481. write(ioimp,*) 'element fissure de longueur nulle ',aa
  482. write(ioimp,*) 'tri3.',iseg,'.',ia,'=',x11,y11,z11
  483. write(ioimp,*) 'tri3.',iseg,'.',ib,'=',x22,y22,z22
  484. write(ioimp,*) 'tri3.',iseg,'.',ic,'=',x33,y33,z33
  485. call erreur(26)
  486. return
  487. endif
  488. sur = bb/2.d0
  489. vx=vx/aa
  490. vy=vy/aa
  491. vz=vz/aa
  492. vxyz = max(abs(vx),abs(vy))
  493. vxyz = max(vxyz,abs(vz))
  494. if (vxyz.le.XPETIT) then
  495. write(ioimp,*)'element fissure de longueur nulle',aa,vx,vy,vz
  496. call erreur(26)
  497. return
  498. endif
  499. xsur = 1.D-8*sur
  500.  
  501. c vecteurs des aretes (suite)
  502. v12lo= sqrt(v12x*v12x+v12y*v12y+v12z*v12z)
  503. v23lo= sqrt(v23x*v23x+v23y*v23y+v23z*v23z)
  504. v31lo= sqrt(v31x*v31x+v31y*v31y+v31z*v31z)
  505. v12x= v12x/v12lo
  506. v12y= v12y/v12lo
  507. v12z= v12z/v12lo
  508. v23x= v23x/v23lo
  509. v23y= v23y/v23lo
  510. v23z= v23z/v23lo
  511. v31x= v31x/v31lo
  512. v31y= v31y/v31lo
  513. v31z= v31z/v31lo
  514.  
  515. if(iimpi.ge.1)
  516. $ write(ioimp,*)'-> triangle ',iseg,'/',nseg,' noeuds ',ia,ib,ic
  517.  
  518.  
  519. *=======================================================================
  520. * Boucle sur les noeuds ======================================
  521. do 490 ipt=1,npt
  522.  
  523. id=ipt3.num(1,ipt)
  524. dis = xgrand
  525. iphi2 = 0
  526.  
  527. *------------ cas ou le noeud structure id = l un des noeuds des tri3 de la fissure
  528. if (ia.eq.id.or.ib.eq.id.or.ic.eq.id) then
  529.  
  530. dis=0.D0
  531. iphi2 = 1
  532.  
  533. *------------ "autres" cas
  534. else
  535.  
  536. * recup des coordonnes du noeud en question
  537. xxp = xcoor((id-1)*idim1 + 1 )
  538. yyp = xcoor((id-1)*idim1 + 2 )
  539. zzp = xcoor((id-1)*idim1 + 3 )
  540.  
  541. *bp2013 on saute ce point si trop loin (critere)
  542. if (icle.eq.1) then
  543. if(xxp.lt.xxmi) goto 491
  544. if(xxp.gt.xxma) goto 491
  545. if(yyp.lt.yymi) goto 491
  546. if(yyp.gt.yyma) goto 491
  547. if(zzp.lt.zzmi) goto 491
  548. if(zzp.gt.zzma) goto 491
  549. endif
  550.  
  551. d1xp=xxp-x11
  552. d1yp=yyp-y11
  553. d1zp=zzp-z11
  554. d2xp=xxp-x22
  555. d2yp=yyp-y22
  556. d2zp=zzp-z22
  557. d3xp=xxp-x33
  558. d3yp=yyp-y33
  559. d3zp=zzp-z33
  560. * calcul de la distance au triangle
  561. d1= sqrt(d1xp*d1xp + d1yp*d1yp + d1zp*d1zp)
  562. d2= sqrt(d2xp*d2xp + d2yp*d2yp + d2zp*d2zp)
  563. d3= sqrt(d3xp*d3xp + d3yp*d3yp + d3zp*d3zp)
  564. if (d1.lt.d2) then
  565. rr= min(d1,d3)
  566. else
  567. rr= min(d2,d3)
  568. endif
  569.  
  570. * test pour eviter plein de calculs => on saute au noeud suivant
  571. if (isup(ipt).eq.0) then
  572. rrlim(ipt) = rr
  573. else
  574. c if(rr.gt.(1.1d0*rrlim(ipt))) goto 491
  575. if(rr.gt.(2.d0*rrlim(ipt))) goto 491
  576. if(rr.lt.rrlim(ipt)) rrlim(ipt)=rr
  577. endif
  578.  
  579. * calcul de la distance au plan du triangle
  580. xdd = vx*d1xp + vy*d1yp + vz*d1zp
  581. * determination du pt projetté
  582. xx=xxp-xdd*vx
  583. yy=yyp-xdd*vy
  584. zz=zzp-xdd*vz
  585. * ce pt est-il à l'interieur ou sur un coté du triangle?
  586. d1x=xx-x11
  587. d1y=yy-y11
  588. d1z=zz-z11
  589. d2x=xx-x22
  590. d2y=yy-y22
  591. d2z=zz-z22
  592. d3x=xx-x33
  593. d3y=yy-y33
  594. d3z=zz-z33
  595. pv12x=d1y*d2z-d1z*d2y
  596. pv12y=d1z*d2x-d1x*d2z
  597. pv12z=d1x*d2y-d1y*d2x
  598. pv23x=d2y*d3z-d2z*d3y
  599. pv23y=d2z*d3x-d2x*d3z
  600. pv23z=d2x*d3y-d2y*d3x
  601. pv31x=d3y*d1z-d3z*d1y
  602. pv31y=d3z*d1x-d3x*d1z
  603. pv31z=d3x*d1y-d3y*d1x
  604. ps12v=pv12X*vx+pv12y*vy+pv12z*vz
  605. ps23v=pv23x*vx+pv23y*vy+pv23z*vz
  606. ps31v=pv31x*vx+pv31y*vy+pv31z*vz
  607.  
  608. * ce pt est-il à l'interieur ou sur un coté du triangle?
  609. if(ps12v.ge.0..and.ps23v.ge.0..and.ps31v.ge.0.) then
  610. * -> on est a l'interieur
  611. dis = xdd
  612. iphi2 = 1
  613. elseif (abs(ps12v).lt.xsur) then
  614. * -> on est situé sur le segment 12 : est-on a l'interieur?
  615. psc1=v12x*d1x+v12y*d1y+v12z*d1z
  616. psc2=v12x*d2x+v12y*d2y+v12z*d2z
  617. if ( psc1*psc2 .lt.0.d0) then
  618. dis = xdd
  619. iphi2 = 1
  620. endif
  621. elseif(abs(ps23v).lt.xsur)then
  622. * -> on est situé sur le segment 23 : est-on a l'interieur?
  623. psc2=v23x*d2x+v23y*d2y+v23z*d2z
  624. psc3=v23x*d3x+v23y*d3y+v23z*d3z
  625. if ( psc2*psc3 .lt.0.d0) then
  626. dis = xdd
  627. iphi2 = 1
  628. endif
  629. elseif(abs(ps31v).lt.xsur)then
  630. * -> on est situé sur le segment 31 : est-on a l'interieur?
  631. psc3=v31x*d3x+v31y*d3y+v31z*d3z
  632. psc1=v31x*d1x+v31y*d1y+v31z*d1z
  633. if ( psc3*psc1 .lt.0.d0) then
  634. dis = xdd
  635. iphi2 = 1
  636. endif
  637. endif
  638.  
  639. * on est a l'exterieur : il faut calculer dis et verifier xdd
  640. if (dis.eq.xgrand) then
  641. iphi2 = 10
  642.  
  643. * si on est a l'ext du tri3 et en avant du front+contour on change v
  644.  
  645. * -si on a affaire a un triangle du front (ou du contour) :
  646. * est on en avant du front (ou du contour)? -> phicor
  647. phicor=.false.
  648. phi12=.false.
  649. phi23=.false.
  650. phi31=.false.
  651. if (ideux.ge.1) then
  652. if (ifis(iseg,1).ge.1) then
  653. * -Cas ou le segment [1-2] = front ou contour
  654. if(ifis(iseg,1).ge.2.or.ifis(iseg,2).eq.1.or.
  655. $ ifis(iseg,2).eq.3.or.ifis(iseg,2).eq.5) then
  656. c phi12=(ps12v.lt.0.d0)
  657. phi12=(ps12v.lt.(-xsur))
  658. if(phi12) then
  659. if(ifis(iseg,1).ge.2) then
  660. iphi2=iphi2+5
  661. else
  662. iphi2=iphi2+1
  663. endif
  664. endif
  665. endif
  666. * -Cas ou le segment [2-3] = contour
  667. if(ifis(iseg,2).eq.2.or.
  668. $ ifis(iseg,2).eq.3.or.ifis(iseg,2).eq.6) then
  669. c phi23=(ps23v.lt.0.d0)
  670. phi23=(ps23v.lt.(-xsur))
  671. if(phi23) iphi2=iphi2+1
  672. endif
  673. * -Cas ou le segment [3-1] = contour
  674. if(ifis(iseg,2).eq.4.or.
  675. $ ifis(iseg,2).eq.5.or.ifis(iseg,2).eq.6) then
  676. c phi31=(ps31v.lt.0.d0)
  677. phi31=(ps31v.lt.(-xsur))
  678. if(phi31) iphi2=iphi2+1
  679. endif
  680. * -Union des test
  681. phicor=(ideux.ge.1).and.(phi12.or.phi23.or.phi31)
  682. endif
  683. endif
  684. * si phicor, il faudra alors "corriger" phi
  685. * (et prendre xdd au lieu de dis)
  686.  
  687. psc1=0.d0
  688. psc2=0.d0
  689. psc3=0.d0
  690. c dis=sign(rr,xdd)
  691. * bp (05/03/2012) : on remplace ci dessus par formule + precise
  692. * --on teste autour du noeud 1--
  693. if (rr.eq.d1) then
  694. if (ps12v.lt.0.d0.and.ps31v.lt.0d0) then
  695. * on est definitivement + proche du noeud 1
  696. dis = sign(rr,xdd)
  697. else
  698. if (ps12v.lt.0.d0) then
  699. * on est du coté du segment 12
  700. psc1 = v12x*d1x + v12y*d1y + v12z*d1z
  701. elseif (ps31v.lt.0.d0) then
  702. * on est du coté du segment 13
  703. psc1 = -1.d0*(v31x*d1x + v31y*d1y + v31z*d1z)
  704. else
  705. * on ne devrait pas passer par la
  706. psc1=0.d0
  707. endif
  708. if (psc1.le.0d0) then
  709. dis = sign(rr,xdd)
  710. else
  711. dis = sign(sqrt(rr*rr-psc1*psc1),xdd)
  712. endif
  713. endif
  714. * --on teste autour du noeud 2--
  715. elseif (rr.eq.d2) then
  716. if (ps12v.lt.0.d0.and.ps23v.lt.0d0) then
  717. * on est definitivement + proche du noeud 2
  718. dis = sign(rr,xdd)
  719. else
  720. if (ps12v.lt.0.d0) then
  721. * on est du coté du segment 21
  722. psc2 = -1.d0*(v12x*d2x + v12y*d2y + v12z*d2z)
  723. elseif (ps31v.lt.0.d0) then
  724. * on est du coté du segment 23
  725. psc2 = v23x*d2x + v23y*d2y + v23z*d2z
  726. else
  727. * on ne devrait pas passer par la
  728. psc2=0.d0
  729. endif
  730. if (psc2.le.0d0) then
  731. dis = sign(rr,xdd)
  732. else
  733. dis = sign(sqrt(rr*rr-psc2*psc2),xdd)
  734. endif
  735. endif
  736. * --on teste autour du noeud 3--
  737. elseif (rr.eq.d3) then
  738. if (ps31v.lt.0.d0.and.ps23v.lt.0d0) then
  739. * on est definitivement + proche du noeud 2
  740. dis = sign(rr,xdd)
  741. else
  742. if (ps31v.lt.0.d0) then
  743. * on est du coté du segment 31
  744. psc3 = v31x*d3x + v31y*d3y + v31z*d3z
  745. elseif (ps31v.lt.0.d0) then
  746. * on est du coté du segment 32
  747. psc3 = -1.d0*(v23x*d3x + v23y*d3y + v23z*d3z)
  748. else
  749. * on ne devrait pas passer par la
  750. psc3=0.d0
  751. endif
  752. if (psc3.le.0d0) then
  753. dis = sign(rr,xdd)
  754. else
  755. dis = sign(sqrt(rr*rr-psc3*psc3),xdd)
  756. endif
  757. endif
  758. endif
  759.  
  760. endif
  761. * fin du cas ou on est a l'exterieur
  762.  
  763. endif
  764. *------------ fin des "autres" cas (noeud id different de ceux du tri3)
  765.  
  766. * -on a calculé (xdd et) dis: s'agit il de phi?
  767. c on prend xdd plutot que dis pour definir Phi dans le cas ou :
  768. c -on a l option DEUX ou TROI
  769. c -le tri3 utilisé (le + proche) appartient au front ou au contour
  770. c -le point se projette en avant du front ou du contour
  771. * iphi2 qu on va eventuellement ranger dans isup(ipt) vaut:
  772. * 1 si on se projette a l int du triangle
  773. * 10 si on se projette a l ext du triangle
  774. * 11 si en + on sort du contour
  775. * 15 si en + on est en avant du front
  776. * 16 si en + on est en avant du front et on sort du contour
  777.  
  778. if(id.eq.idbug) then
  779. c if(phicor) then
  780. write(6,*) 'noeud:',id,' (',ipt,') tri3',iseg,':',ia,ib,ic
  781. write(6,*) 'iphi2,isup(ipt),dis,xdd',iphi2,isup(ipt),dis,xdd
  782. write(6,*) 'phi12,phicor=',phi12,phicor
  783. endif
  784.  
  785. * -c'est la premiere fois que l'on voit ce point
  786. if (isup(ipt).eq.0) then
  787. if(phicor) then
  788. mpova2.vpocha(ipt,1)=xdd
  789. else
  790. mpova2.vpocha(ipt,1)=dis
  791. endif
  792. xdmin(ipt) =dis
  793. xdmin1(ipt)=xdd
  794. isup(ipt)=iphi2
  795. iphi2=100
  796. igard=igard+1
  797. else
  798. * --on est (franchement) + proche que les precedents essais
  799. if (abs(dis).lt.0.999999*abs(xdmin(ipt))) then
  800. if(phicor) then
  801. mpova2.vpocha(ipt,1)=xdd
  802. else
  803. mpova2.vpocha(ipt,1)=dis
  804. endif
  805. xdmin(ipt) =dis
  806. xdmin1(ipt)=xdd
  807. isup(ipt)=iphi2
  808. iphi2=100
  809. * --on est (quasiment) aussi proche que le meilleur essai
  810. elseif (abs(dis).lt.1.000001*abs(xdmin(ipt))) then
  811. * -on a 2 tri3 du front+contour a egalite
  812. * et l'un des deux a trouvé le point en avant du front :
  813. if ((iseg.le.imil5).and.
  814. $ (phicor.or.isup(ipt).gt.10)) then
  815. c strategie : xdd= max(xdd1 xdd2)
  816. if (abs(xdd).gt.abs(xdmin1(ipt))) then
  817. mpova2.vpocha(ipt,1)=xdd
  818. xdmin(ipt) =dis
  819. xdmin1(ipt)=xdd
  820. c isup(ipt)=iphi2
  821. if(iphi2.gt.isup(ipt)) isup(ipt)=iphi2
  822. iphi2=100
  823. c cas ou on avait pas vu qu il fallait prendre xdd
  824. elseif(isup(ipt).le.10) then
  825. c equivalent à :
  826. c elseif(mpova2.vpocha(ipt,1).ne.xdmin1(ipt)) then
  827. mpova2.vpocha(ipt,1)=xdmin1(ipt)
  828. isup(ipt)=iphi2
  829. iphi2=100
  830. endif
  831.  
  832. * -on avait 1 tri3 du front+contour et ce tri3 est a egalite
  833. * et le 1er avait trouvé le point en avant du front :
  834. * on laisse le xdd calculé depuis le 1er tri3
  835. elseif ((iseg.gt.imil5).and.
  836. $ (isup(ipt).gt.10)) then
  837. goto 491
  838. * -quasi-egalite entre 2 concurrents mais signes opposés :
  839. else
  840. ss1=dis*xdmin(ipt)
  841. if (ss1.le.0.d0) then
  842. * on choisit le dis associé au xdd le + marqué
  843. if (abs(xdd).gt.abs(xdmin1(ipt))) then
  844. mpova2.vpocha(ipt,1)=dis
  845. xdmin(ipt) =dis
  846. xdmin1(ipt)=xdd
  847. c isup(ipt)=iphi2
  848. if(iphi2.gt.isup(ipt)) isup(ipt)=iphi2
  849. iphi2=100
  850. endif
  851. endif
  852. endif
  853. endif
  854. * --fin de la distinction + proche / aussi proche
  855. endif
  856. * -fin de la distinction 1ere fois ou pas
  857.  
  858. if(id.eq.idbug) then
  859. write(6,*)'iphi2,xdmin,xdmin1,phi',iphi2,
  860. $ xdmin(ipt),xdmin1(ipt),mpova2.vpocha(ipt,1)
  861. endif
  862. if(iimpi.ge.2) then
  863. write(ioimp,*) 'point',id,ipt,'dis,xdd,PHI,isup=',
  864. $ dis,xdd,mpova2.vpocha(ipt,1),isup(ipt)
  865. endif
  866.  
  867.  
  868. 491 continue
  869. * le tri3 n'est pas un tri3 du front => on saute
  870. *TODO : test a verifier peut etre .eq. et .eq.-1 aussi a ne pas sauter?
  871. if (ideux.lt.1) goto 490
  872. if (ifis(iseg,1).le.1) goto 490
  873.  
  874. c ------ Calcul de la 2eme Level Set -------
  875. xdd1 = 0.D0
  876.  
  877. * -cas ou le noeud structure id = l un des noeuds des tri3 de la fissure
  878. if (ia.eq.id.or.ib.eq.id) then
  879. dis1=0.D0
  880. xd2=0.D0
  881. xd3=0.D0
  882. if(ideux.ge.2) xtau=xicpr(id)
  883. * -autres cas
  884. else
  885. * on recupere ce qu on a pas recupéré lors du calcul de phi
  886. if (ic.eq.id) then
  887. xxp = xcoor((id-1)*idim1 + 1 )
  888. yyp = xcoor((id-1)*idim1 + 2 )
  889. zzp = xcoor((id-1)*idim1 + 3 )
  890. endif
  891. * calcul du vecteur vn qui est le sens de propagation
  892. vnx=v12y*vz-v12z*vy
  893. vny=v12z*vx-v12x*vz
  894. vnz=v12x*vy-v12y*vx
  895. * il faut orienter vnx pour que le vecteur soit dans le sens
  896. * de la propagation
  897. qq= vnx*(x33-x11)+vny*(y33-y11)+vnz*(z33-z11)
  898. if (qq.ge.0.D0) then
  899. vnx=-vnx
  900. vny=-vny
  901. vnz=-vnz
  902. endif
  903. * normalement vn est normé on le verifie
  904. qq= sqrt(vnx*vnx + vny*vny + vnz*vnz)
  905. if (abs(qq-1.d0).gt.1.d-2) then
  906. write(ioimp,*) ' pb dans le calcul de vn'
  907. call erreur(26)
  908. return
  909. elseif(abs(qq-1.).gt.1.d-6) then
  910. vnx=vnx/qq
  911. vny=vny/qq
  912. vnz=vnz/qq
  913. endif
  914. * on calcule la distance au plan defini par segment 12 et vecteur vx,vy,vz
  915. * et le point (xx yy zz ) projeté du point sur ce plan
  916. xdd1=(xxp-x11)*vnx+(yyp-y11)*vny+(zzp-z11)*vnz
  917. xx=xxp-xdd1*vnx
  918. yy=yyp-xdd1*vny
  919. zz=zzp-xdd1*vnz
  920. psc1=(xx-x11)*v12x+(yy-y11)*v12y+(zz-z11)*v12z
  921. psc2=(xx-x22)*v12x+(yy-y22)*v12y+(zz-z22)*v12z
  922. if ((psc1*psc2).le.0.d0) then
  923. C le point se projete sur 12
  924. dis1=xdd1
  925. pp1x=x11+(psc1*v12x)
  926. pp1y=y11+(psc1*v12y)
  927. pp1z=z11+(psc1*v12z)
  928. xd2 = 0.D0
  929. xd3 = sqrt((xxp-pp1x)**2+(yyp-pp1y)**2+(zzp-pp1z)**2)
  930. if (ideux.ge.2) then
  931. xtau1=xicpr(ia)
  932. xtau2=xicpr(ib)
  933. stau = psc1 / v12lo
  934. xtau=stau*xtau2+(1.d0-stau)*xtau1
  935. endif
  936. else
  937. if (psc1.ge.0.D0) then
  938. pp2x=xx+((v12lo-psc1)*v12x)
  939. pp2y=yy+((v12lo-psc1)*v12y)
  940. pp2z=zz+((v12lo-psc1)*v12z)
  941. xd2 = abs(v12lo-psc1)
  942. xd3 = sqrt((xxp-x22)**2+(yyp-y22)**2+(zzp-z22)**2)
  943. c if(ideux.ge.2) xtau=xicpr(ib)
  944. if(ideux.ge.2) then
  945. if(ib.eq.i2firs.or.ib.eq.i2last) then
  946. c write(6,*) 'le noeud',id,'detecte comme apres',ib
  947. xtau1=xicpr(ia)
  948. xtau2=xicpr(ib)
  949. stau = psc1 / v12lo
  950. xtau=stau*xtau2+(1.d0-stau)*xtau1
  951. else
  952. xtau=xicpr(ib)
  953. endif
  954. endif
  955. else
  956. pp2x=xx-(psc1*v12x)
  957. pp2y=yy-(psc1*v12y)
  958. pp2z=zz-(psc1*v12z)
  959. xd2 = abs(psc1)
  960. xd3 = sqrt((xxp-x11)**2+(yyp-y11)**2+(zzp-z11)**2)
  961. c if(ideux.ge.2) xtau=xicpr(ia)
  962. if(ideux.ge.2) then
  963. if(ia.eq.i2firs.or.ia.eq.i2last) then
  964. c write(6,*) 'le noeud',id,'detecte comme avant',ia
  965. xtau1=xicpr(ia)
  966. xtau2=xicpr(ib)
  967. stau = psc1 / v12lo
  968. xtau=stau*xtau2+(1.d0-stau)*xtau1
  969. else
  970. xtau=xicpr(ia)
  971. endif
  972. endif
  973. endif
  974. dis1=sqrt((xxp-pp2x)**2+(yyp-pp2y)**2+(zzp-pp2z)**2)
  975. dis1=sign(dis1,xdd1)
  976. endif
  977.  
  978. endif
  979. * -fin de la distinction si noeud du tri3 ou pas
  980.  
  981. * 1ere fois qu on traite ce point?
  982. if (isupfi(ipt).eq.0) then
  983. * isupfi(ipt)=1
  984. * on utilise isupfi pour definir si la valeur inscrite
  985. * provient du 1er ou dernier seg du front (=2) ou autre(=1)?
  986. isupfi(ipt)=ifis(iseg,1)
  987. igardf=igardf+1
  988. xdmin2(ipt)= xd2
  989. xdmin3(ipt)= xd3
  990. cbp2013 mpova1.vpocha(ipt,1)=dis1
  991. if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau
  992. else
  993. * if(abs(dis1).lt.abs(mpova1.vpocha(ipt,1)))then
  994. * on a trouvé un segment + proche ou = du noeud considéré
  995. * if(abs(dis1).le.abs(mpova1.vpocha(ipt,1)))then
  996. ** if(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))-1.D-8)then
  997. * if(xd2.le.xdmin2(ipt)) then
  998. if (xd3.lt.(abs(xdmin3(ipt))-1.D-8)) then
  999. cbp2013 mpova1.vpocha(ipt,1)=dis1
  1000. xdmin2(ipt)=xd2
  1001. xdmin3(ipt)=xd3
  1002. isupfi(ipt)=ifis(iseg,1)
  1003. if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau
  1004. endif
  1005. * rattrapage dans le cas d'erreur d'arrondi sur dis1
  1006. * elseif(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))+1.D-8)
  1007. * $ then
  1008. * if(abs(dis1).lt.(abs(mpova1.vpocha(ipt,1)))+1.D-8)then
  1009. if (xd3.lt.(abs(xdmin3(ipt))+1.D-8))then
  1010. * en cas d'égalité,on n'utilise la distance dis1 que si:
  1011. * - le pt projetté appartient au segment (xd2=0)
  1012. * ou - le point projeté est + proche du segment
  1013. if (xd2.lt.xdmin2(ipt)) then
  1014. cbp2013 mpova1.vpocha(ipt,1)=dis1
  1015. xdmin2(ipt)=xd2
  1016. xdmin3(ipt)=xd3
  1017. isupfi(ipt)=ifis(iseg,1)
  1018. if(ideux.ge.2) mpova3.vpocha(ipt,1)=xtau
  1019. endif
  1020. endif
  1021. endif
  1022.  
  1023. 490 continue
  1024. * Fin de Boucle sur les noeuds ===============================
  1025. *=======================================================================
  1026.  
  1027. 400 continue
  1028. * Fin de Boucle sur les elements de la fissure =========================
  1029. *=======================================================================
  1030.  
  1031.  
  1032. *------Petite correction de Psi
  1033. if (ideux.ge.1) then
  1034. * Boucle sur les noeuds ===============================
  1035. do 500 ipt=1,npt
  1036. if (isup(ipt).eq.0) goto 500
  1037. c * pour le cas de fissure debouchante
  1038. c * dans le cas ou c'est un segment debouchant qui a servi au calcul de psi,
  1039. c * il faut prendre la distance xdd2 (et pas dis) qu'on retrouve avec :
  1040. c if ((isupfi(ipt)).eq.2) then
  1041. c dis1 = mpova1.vpocha(ipt,1)
  1042. c xdd1 = sqrt( (dis1**2) - (xdmin2(ipt)**2) )
  1043. c mpova1.vpocha(ipt,1) = sign(xdd1,dis1)
  1044. c endif
  1045. c bp : essai du 16/03/2012 : psi = distance au front**2-phi**2
  1046. xdd1 = (xdmin3(ipt)**2) - (mpova2.vpocha(ipt,1)**2)
  1047. xdd1 = sqrt(max(xdd1,0.d0))
  1048. if(ipt3.num(1,ipt).eq.idbug) then
  1049. write(6,*)'ipt,isup,psi_avant,phi,xdd1,= ',ipt,isup(ipt),
  1050. $ xdmin3(ipt),mpova1.vpocha(ipt,1),mpova2.vpocha(ipt,1),xdd1
  1051. endif
  1052. c if (isup(ipt).gt.0) then
  1053. if (isup(ipt).ge.15) then
  1054. *on est a l exterieur => psi >0
  1055. mpova1.vpocha(ipt,1) = xdd1
  1056. else
  1057. mpova1.vpocha(ipt,1) = -1.d0*xdd1
  1058. endif
  1059. 500 continue
  1060. endif
  1061. *------Fin de la Petite correction de Psi
  1062.  
  1063.  
  1064. ***********************************************************************
  1065. * FIN DU TRAVAIL
  1066. ***********************************************************************
  1067. if(iimpi.ge.1) then
  1068. call gibtem(xkt)
  1069. write(ioimp,*)'fin du travail | xkt=',xkt
  1070. endif
  1071.  
  1072. * on inverse isup
  1073. * AVANT : isup(i) =0 veut dire que le point na pas été traité
  1074. * APRES : isup(i) =1 veut dire que le point peut etre enleve (=0 sinon)
  1075. iaju=0
  1076. do ia=1,isup(/1)
  1077. if (isup(ia).eq.0) then
  1078. iaju=iaju+1
  1079. isup(ia)=1
  1080. else
  1081. isup(ia)=0
  1082. endif
  1083. enddo
  1084. if (iaju+igard.ne.isup(/1)) then
  1085. write(ioimp,*) 'psiphi : tout va mal '
  1086. $ ,iaju,'+',igard,'ne',isup(/1)
  1087. call erreur(26)
  1088. return
  1089. endif
  1090. ipt1=ipt3
  1091. ipt4=melfis
  1092. segdes ipt4
  1093. * iaju est le nombre de points a éliminer
  1094. * isup(i) =1 veut dire que le point peut etre enleve (=0 sinon)
  1095.  
  1096. * faut-il ajuster les champs?
  1097. if (iaju.ne.0) then
  1098. * si oui on recree les bons chpoints
  1099. nbelem=ipt1.num(/2)- iaju
  1100. nbnn=1
  1101. nbref=0
  1102. nbsous=0
  1103. segini ipt2
  1104. ipt2.itypel=1
  1105. nc=1
  1106. n=nbelem
  1107. segini mpova5
  1108. *>>>bp&yt: ajout ligne suivante
  1109. if(ideux.ge.1) segini mpova4
  1110. if(ideux.ge.2) segini mpova6
  1111. iel=0
  1112. do i=1,ipt1.num(/2)
  1113. if (isup(i).eq.0) then
  1114. iel=iel+1
  1115. ipt2.num(1,iel)=ipt1.num(1,i)
  1116. mpova5.vpocha(iel,1)=mpova2.vpocha(i,1)
  1117. *>>>bp&yt: ajout ligne suivante
  1118. * on met le meme maillage pour les 2 level set (car toujours utilisées ensemble)
  1119. if(ideux.ge.1) mpova4.vpocha(iel,1)=mpova1.vpocha(i,1)
  1120. if(ideux.ge.2) mpova6.vpocha(iel,1)=mpova3.vpocha(i,1)
  1121. endif
  1122. enddo
  1123. if (ideux.ge.2) then
  1124. c segact msoup3*mod
  1125. msoup3.ipoval=mpova6
  1126. msoup3.igeoc=ipt2
  1127. segdes msoup3,mpova6
  1128. segsup,mpova3
  1129. endif
  1130. if (ideux.ge.1) then
  1131. c segact msoup1*mod
  1132. msoup1.ipoval=mpova4
  1133. msoup1.igeoc=ipt2
  1134. segdes msoup1,mpova4
  1135. segsup,mpova1
  1136. endif
  1137. c segact msoup2*mod
  1138. msoup2.ipoval=mpova5
  1139. msoup2.igeoc=ipt2
  1140. segdes msoup2,mpova5,ipt2
  1141. segsup,mpova2
  1142. segdes ipt1
  1143. else
  1144. if(ideux.ge.2) segdes,msoup3,mpova3
  1145. if(ideux.ge.1) segdes,msoup1,mpova1
  1146. segdes msoup2,mpova2,ipt1
  1147. endif
  1148. if (ideux.ge.1) then
  1149. segsup,icpr,sfis
  1150. if(ideux.ge.2) segsup,xicpr
  1151. endif
  1152. if (icle.eq.1) then
  1153. segsup isup,isupfi,xdis,xdmin,xdmin1,rrlim
  1154. if(ideux.ge.1) segsup xdmin2,xdmin3
  1155. endif
  1156.  
  1157. RETURN
  1158. END
  1159.  
  1160.  
  1161.  
  1162.  
  1163.  
  1164.  

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