Télécharger psip2d.eso

Retour à la liste

Numérotation des lignes :

  1. C PSIP2D SOURCE CB215821 19/08/20 21:21:06 10287
  2. c
  3. SUBROUTINE PSIP2D(ideux,melmai,melfis,ip1,ip2,xcrit,
  4. & msoup1,msoup2)
  5. ***********************************************************************
  6. c OPERATEUR : PSIP
  7. c
  8. c APPELE PAR PSIPHI DANS LE CAS 2D
  9. c
  10. c CREATION : BP 14/03/2012
  11. c MODIF : BP 18/12/2013 : plus d'appel a zonag2, mais travail sur critere direct
  12. c
  13. ***********************************************************************
  14. IMPLICIT REAL*8(A-H,O-Z)
  15. IMPLICIT INTEGER(I-N)
  16. -INC CCOPTIO
  17. -INC SMCOORD
  18. -INC SMELEME
  19. -INC SMCHPOI
  20. -INC CCREEL
  21. segment ttrav
  22. integer ilig(inoe)
  23. endsegment
  24. segment isup(npt)
  25. c segment xdis(npt)
  26. c segment xdmin(npt)
  27. c segment xdmin1(npt)
  28. c segment rrlim(npt)
  29. c segment xdmin2(npt)
  30. c segment xdmin3(npt)
  31. c segment xrr2b(npt)
  32. c segment xpsc1b(npt)
  33. c segment isupfi(npt)
  34. c SEGMENT ISEG3
  35. c INTEGER NIZO(NZO+1)
  36. c ENDSEGMENT
  37. c c SEGMENT ISEG4
  38. c c INTEGER NUMZO(NZO)
  39. c c ENDSEGMENT
  40. c SEGMENT ISEG5
  41. c INTEGER NNMEL(ILON),IDEJ(NZO)
  42. c ENDSEGMENT
  43. DATA EPSI/1.D-5/
  44. c REAL*8 XZONE(3,3)
  45. c INTEGER NZONE(3)
  46. idim1=idim+1
  47.  
  48.  
  49. ***********************************************************************
  50. * Initialisation, Recup, Activation
  51. ***********************************************************************
  52. idim1=idim+1
  53. idbug = 0
  54. if(msoup2.ne.0) mpova2=msoup2.ipoval
  55. if(msoup1.ne.0) mpova1=msoup1.ipoval
  56. c activation des maillages / desactivation
  57. c melmai : fait dans psiphi / ligne 1325 ou 1329
  58. c melfis : fait ligne 56 / ligne 72
  59.  
  60. ***********************************************************************
  61. * Preparation du maillage de la fissure
  62. ***********************************************************************
  63. meleme=melfis
  64. segact meleme
  65.  
  66. if (lisous(/1).ne.0) then
  67. call erreur(25)
  68. return
  69. endif
  70. if (itypel.ne.2.and.itypel.ne.3) then
  71. * seg2 et seg3 seulement!
  72. call erreur(16)
  73. return
  74. endif
  75. * on ordonne la ligne depuis la pointe de fissure
  76. call ligmai(meleme,ttrav,0)
  77. if(ierr.ne.0) return
  78. segact ttrav*mod
  79. * on regarde si ip1 est bien en dernier
  80. nb=ilig(/1)
  81. if (ideux.ge.1) then
  82. if (ilig(nb).eq.ip1) then
  83. * rien a faire
  84. else
  85. * il faut inverser le sens de parcours de la ligne
  86. nc=nb+1
  87. nn=nb/2
  88. do i=1,nn
  89. io=ilig(i)
  90. iu=ilig(nc-i)
  91. ilig(i)=iu
  92. ilig(nc-i)=io
  93. enddo
  94. endif
  95. endif
  96. * on verifie que ip2 est bien le premier point
  97. if (ip2.ne.0) then
  98. if (ilig(1).ne.ip2) then
  99. call erreur (21)
  100. endif
  101. endif
  102. nseg=nb-1
  103.  
  104. * sont actifs: melmai (=struct changée en poi1)
  105.  
  106.  
  107. ***********************************************************************
  108. * debut du travail suivant valeur de icle
  109. if (xcrit.eq.0.d0) then
  110. icle=0
  111. xcrit=1.D10
  112. else
  113. icle=1
  114. * ajout d'un epsilon 1.D-5 de rattrapage
  115. xcrit = 1.00001D0 * xcrit
  116. endif
  117.  
  118. if(iimpi.ge.1) then
  119. call gibtem(xkt)
  120. write(ioimp,*)'debut zonage si icle=1=',icle,xcrit,'| xkt=',xkt
  121. endif
  122.  
  123. c ***** cas icle=1 ***************************************************
  124. c if (icle.eq.1) then
  125. c
  126. c c write(6,*) '--cas icle=1'
  127. c
  128. c melfis=meleme
  129. c melpoi=melmai
  130. c
  131. c *---- ZONAGE -------------------------------------------------------
  132. c call zonag2(melfis,melpoi,xcrit,ISEG3,ISEG5,XZONE,NZONE)
  133. c * recup
  134. c ILON=NNMEL(/1)
  135. c NZO=IDEJ(/1)
  136. c Xtomi=XZONE(1,1)
  137. c Ytomi=XZONE(2,1)
  138. c c Ztomi=XZONE(3,1)
  139. c Xpr=XZONE(1,2)
  140. c Ypr=XZONE(2,2)
  141. c c Zpr=XZONE(3,2)
  142. c XZO=XZONE(1,3)
  143. c YZO=XZONE(2,3)
  144. c c ZZO=XZONE(3,3)
  145. c nXzo=NZONE(1)
  146. c nYzo=NZONE(2)
  147. c c nZzo=NZONE(3)
  148. c
  149. c
  150. c ***********************************************************************
  151. c * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT
  152. c ***********************************************************************
  153. c
  154. c * pour chaque segment de la fissure (melfis 2eme maillage lu)
  155. c * on calcule le champ phi et psi des noeuds autour de l'element
  156. c if(iimpi.ge.1) then
  157. c call gibtem(xkt)
  158. c write(ioimp,*)'debut travail | xkt=',xkt
  159. c endif
  160. c *
  161. c meleme=melfis
  162. c c segact,meleme
  163. c ipt3=melpoi
  164. c c segact,ipt3
  165. c npt=ipt3.num(/2)
  166. c * attention on initialise isup avec les valeurs inverses,il faudra
  167. c * le retourner
  168. c segini isup,isupfi,xdis,xdmin,xdmin1,rrlim
  169. c if(ideux.ge.1) segini xdmin2,xdmin3
  170. c if(ip2.ne.0) segini xrr2b,xpsc1b
  171. c igard=0
  172. c igardf=0
  173. c * Boucle sur les elements de la fissure =========================
  174. c
  175. c do 400 iseg=1,nseg
  176. c c if(idebug.ge.1)
  177. c c $ write(6,*) '====== element fissure ',iseg,' / ',nseg,ifis(iseg)
  178. c
  179. c ia=ilig(iseg)
  180. c ib=ilig(iseg+1)
  181. c x11=xcoor((ia-1)*idim1+1)
  182. c y11=xcoor((ia-1)*idim1+2)
  183. c x22=xcoor((ib-1)*idim1+1)
  184. c y22=xcoor((ib-1)*idim1+2)
  185. c vx=x22-x11
  186. c vy=y22-y11
  187. c * aa= sqrt(vx*vx + vy*vy) + 1.d-20
  188. c * bp: dans ce cas, mieux vaut faire une erreur
  189. c aa= sqrt(vx*vx + vy*vy)
  190. c if (aa.le.XPETIT) then
  191. c write(ioimp,*) 'ERREUR: segment de fissure de longueur nulle!'
  192. c call erreur(26)
  193. c return
  194. c endif
  195. c vx=vx/aa
  196. c vy=vy/aa
  197. c vxyz = max(abs(vx),abs(vy))
  198. c if (vxyz.le.XPETIT) then
  199. c write(ioimp,*) 'segment de fissure de longueur nulle! '
  200. c $ ,aa,vx,vy
  201. c call erreur(26)
  202. c return
  203. c endif
  204. c if (x11.gt.x22) then
  205. c xxmi=x22
  206. c xxma=X11
  207. c else
  208. c xxma=x22
  209. c xxmi=x11
  210. c endif
  211. c if (y11.gt.y22) then
  212. c yymi=y22
  213. c yyma=y11
  214. c else
  215. c yyma=y22
  216. c yymi=y11
  217. c endif
  218. c * calcul des zones xmi,ymi et xma,yma
  219. c NIZ1X=INT((Xxmi-XTOMI-XPR)/XZO) +1
  220. c NIZ1Y=INT((yymi-YTOMI-YPR)/YZO)+1
  221. c NIZ2X=INT((xxma-XTOMI+XPR)/XZO) +1
  222. c NIZ2Y=INT((yyma-YTOMI+YPR)/YZO) +1
  223. c * on ajoute une zone de chaque cotés
  224. c NIZ1X=max(1,niz1x-1)
  225. c niz1y=max(niz1y-1,1)
  226. c niz2x=min(nxzo,niz2x+1)
  227. c niz2y=min(nyzo,niz2y+1)
  228. c c write(6,*)'bornes x y z ',niz1x,niz2x,niz1y,niz2y,NIZ1Z,NIZ2Z
  229. c c ncalc = 0
  230. c
  231. c * on explore les zones ------------------------------------------------
  232. c do 332 m1=niz1y,niz2y
  233. c do 332 m2=niz1x,niz2x
  234. c indzo=m2+(m1-1)*nxzo
  235. c IDEB=NIZO(INDZO)
  236. c IFIN=NIZO(INDZO+1)-1
  237. c if(ideb.gt.ifin) go to 332
  238. c do kk=ideb,ifin
  239. c ipt=NNMEL(KK)
  240. c ia=ipt3.num(1,ipt)
  241. c xx = xcoor((ia-1)*idim1 + 1 )
  242. c yy = xcoor((ia-1)*idim1 + 2 )
  243. c * calcul de la distance au premier pt
  244. c rr2=(xx-x11)*(xx-x11)+(yy-y11)*(yy-y11)
  245. c * calcul des produits scals et du produit vectoiel
  246. c psca1= vx*(xx-x11)+vy*(yy-y11)
  247. c psca2= vx*(xx-x22)+vy*(yy-y22)
  248. c pvec=vx*(yy-y11)-vy*(xx-x11)
  249. c xsi=psca1*psca2
  250. c xd= rr2 - psca1*psca1
  251. c if ((iseg.eq.1).and.(ip2.ne.0)) then
  252. c * dans le cas de 2 pointes on sauve la distance a la 1ere pointe
  253. c xrr2b(ipt) = rr2
  254. c xpsc1b(ipt) = psca1
  255. c endif
  256. c if (xsi.le.0.d0) then
  257. c * on est dans le cas ou le point se projette sur le segment
  258. c if(xd.lt.0.D0) xd=0.D0
  259. c if (iseg.eq.1) then
  260. c xdmin(ipt)= xd
  261. c xdmin1(ipt)=pvec
  262. c xdis(ipt) = sign((sqrt(xd)),pvec)
  263. c else
  264. c if (xd.le.xdmin(ipt))then
  265. c xdmin(ipt)= xd
  266. c xdmin1(ipt)=pvec
  267. c xdis(ipt) = sign((sqrt(xd)),pvec)
  268. c endif
  269. c endif
  270. c else
  271. c * le point ne se projette pas sur le segment
  272. c * est-il plus proche du premier ou du dernier pt du segment
  273. c if (psca1.ge.0.D0)then
  274. c rr2=(xx-x22)*(xx-x22)+(yy-y22)*(yy-y22)
  275. c xd= rr2 - psca2*psca2
  276. c endif
  277. c if(xd.lt.0.D0) xd=0.D0
  278. c if (iseg.eq.1) then
  279. c xdmin(ipt)= rr2
  280. c xdmin1(ipt)=pvec
  281. c xdis(ipt) = sign((sqrt(xd)),pvec)
  282. c else
  283. c c if (rr2.le.xdmin(ipt)) then
  284. c if (rr2.lt.0.999999*xdmin(ipt)) then
  285. c xdmin(ipt)= rr2
  286. c xdmin1(ipt)=pvec
  287. c if (iseg.eq.nseg) then
  288. c xdis(ipt) = sign((sqrt(xd)),pvec)
  289. c else
  290. c xdis(ipt) = sign((sqrt(rr2)),pvec)
  291. c endif
  292. c elseif (rr2.lt.1.000001*xdmin(ipt)) then
  293. c ss1=pvec*xdmin(ipt)
  294. c if (ss1.le.0.d0) then
  295. c * bp 20.06.2011 : egalite entre 2 concurrents mais signes opposés !
  296. c if (abs(pvec).gt.abs(xdmin1(ipt))) then
  297. c xdmin(ipt)= rr2
  298. c xdmin1(ipt)=pvec
  299. c if (iseg.eq.nseg) then
  300. c xdis(ipt) = sign((sqrt(xd)),pvec)
  301. c else
  302. c xdis(ipt) = sign((sqrt(rr2)),pvec)
  303. c endif
  304. c endif
  305. c endif
  306. c endif
  307. c endif
  308. c endif
  309. c if ( isup(ipt).eq.0) then
  310. c * c'est la premiere fois que l'on voit ce point
  311. c isup(ipt)=1
  312. c igard=igard+1
  313. c endif
  314. c * si c'est le dernier segment on stocke les valeurs phi
  315. c if (iseg.eq.nseg) then
  316. c igardf=igardf+1
  317. c isupfi(ipt)=1
  318. c mpova2.vpocha(ipt,1)=xdis(ipt)
  319. c * puis on calcule les valeurs de psi
  320. c * ...relatif a la 1ere pointe de fissure
  321. c if (ideux.ge.1) then
  322. c rr2 = (xx-x22)*(xx-x22)+(yy-y22)*(yy-y22)
  323. c if (psca2.ge.0.) then
  324. c mpova1.vpocha(ipt,1)=psca2
  325. c else
  326. c yd = (rr2 - xdis(ipt)*xdis(ipt))
  327. c if(yd.lt.0.D0) yd=0.D0
  328. c mpova1.vpocha(ipt,1)= -1.*(yd**0.5)
  329. c endif
  330. c * ...relatif a la 2eme pointe de fissure
  331. c if (ip2.ne.0) then
  332. c if (xpsc1b(ipt).le.0.) then
  333. c psi2= -1.*xpsc1b(ipt)
  334. c else
  335. c psi2= (xrr2b(ipt) - xdis(ipt)*xdis(ipt))
  336. c if(psi2.lt.0.D0) then
  337. c psi2= 0.D0
  338. c else
  339. c psi2= -1.*(psi2**0.5)
  340. c endif
  341. c endif
  342. c * et on choisit le psi le plus proche
  343. c cbp2013 if(psi2.gt.(mpova1.vpocha(ipt,1)))
  344. c if(rr2.gt.xrr2b(ipt))
  345. c & mpova1.vpocha(ipt,1) = psi2
  346. c endif
  347. c endif
  348. c endif
  349. c enddo
  350. c 332 continue
  351. c *
  352. c
  353. c c if (idebug.ge.1) then
  354. c c ncalc2 = ncalc2 + ncalc
  355. c c if (mod(iseg,1000).eq.0.or.iseg.eq.nseg) then
  356. c c call gibtem(xkt)
  357. c c write(6,*) '---------------xkt=',ncalc2,xkt
  358. c c ncalc2=0
  359. c c endif
  360. c c endif
  361. c
  362. c 400 continue
  363. c * Fin de Boucle sur les elements de la fissure =========================
  364. c
  365. c
  366. c * on inverse isup
  367. c iaju=0
  368. c do ia=1,isup(/1)
  369. c if (isup(ia).eq.0) then
  370. c iaju=iaju+1
  371. c isup(ia)=1
  372. c else
  373. c isup(ia)=0
  374. c endif
  375. c enddo
  376. c if (iaju+igard.ne.isup(/1)) then
  377. c write(ioimp,*) 'psiphi : tout va mal '
  378. c $ ,iaju,'+',igard,'ne',isup(/1)
  379. c call erreur(26)
  380. c return
  381. c endif
  382. c ipt1=ipt3
  383. c c segsup iseg3,iseg4,iseg5
  384. c segsup iseg3,iseg5
  385. c
  386. c
  387. c ***** cas icle=0 ***************************************************
  388. c else
  389.  
  390. c * recherche d'une boite autour de la fissure
  391. c if (ircrit.ne.0) then
  392. c xcrit= xcrit*2.1
  393. c xmin=xgrand
  394. c xmax=-xgrand
  395. c ymin=xgrand
  396. c ymax=-xgrand
  397. c do i=1,nb
  398. c ia=ilig(i)
  399. c xmin= min( xmin,xcoor((ia-1)*idim1+1))
  400. c ymin= min( ymin,xcoor((ia-1)*idim1+2))
  401. c xmax= max( xmax,xcoor((ia-1)*idim1+1))
  402. c ymax= max( ymax,xcoor((ia-1)*idim1+2))
  403. c enddo
  404. c xmin = xmin-xcrit
  405. c xmax=xmax+xcrit
  406. c ymin=ymin-xcrit
  407. c ymax=ymax+xcrit
  408. c * write(6,*) ' xcrit,xmin,xmax,ymin,ymax'
  409. c * write(6,*) xcrit,xmin,xmax,ymin,ymax
  410. c endif
  411. *
  412. * il faut maintenant remplir les champs
  413. *
  414. * pour chaque pt de ipt1 on calcule le produit scalaire et vectoriel
  415. * par rapport à chaque segment de la ligne
  416. ipt1=melmai
  417. npt=ipt1.num(/2)
  418. segini isup
  419. igard=0
  420. *
  421. * boucle sur les points de ipt1 ****************************************
  422. *
  423. do 100 ipt=1,ipt1.num(/2)
  424.  
  425. * recup du ipt ieme point #ia : x
  426. ia= ipt1.num(1,ipt)
  427. xx = xcoor((ia-1)*idim1 + 1 )
  428. yy = xcoor((ia-1)*idim1 + 2 )
  429. c write(6,*) 'point',ipt,' #',ia,' x=',xx,' y=',yy
  430.  
  431. c qq initialisations
  432. ivu = 0
  433. rr22=XGRAND
  434. psca1b=XGRAND
  435.  
  436. * boucle sur les segments ----------------------------------------
  437. do 101 iseg=1,nseg
  438.  
  439. * recup du segment 12
  440. ia=ilig(iseg)
  441. ib=ilig(iseg+1)
  442. x11=xcoor((ia-1)*idim1+1)
  443. y11=xcoor((ia-1)*idim1+2)
  444. x22=xcoor((ib-1)*idim1+1)
  445. y22=xcoor((ib-1)*idim1+2)
  446. c write(6,*)'segment fissure',iseg,' 1:',x11,y11,' 2:',x22,y22
  447. if (icle.eq.1) then
  448. xmin = min(x11,x22) - xcrit
  449. xmax = max(x11,x22) + xcrit
  450. ymin = min(y11,y22) - xcrit
  451. ymax = max(y11,y22) + xcrit
  452. endif
  453.  
  454. *bp2013 on ne saute pas le segments du front (=1er et dernier segments)
  455. if(iseg.eq.1.or.iseg.eq.nseg) goto 111
  456. *bp2013 on saute ce segment si point trop loin (critere)
  457. if (icle.eq.1) then
  458. if(xx.lt.xmin) goto 101
  459. if(xx.gt.xmax) goto 101
  460. if(yy.lt.ymin) goto 101
  461. if(yy.gt.ymax) goto 101
  462. endif
  463.  
  464. 111 continue
  465. * vecteur tangent v=12/|12|
  466. vx=x22-x11
  467. vy=y22-y11
  468. aa= sqrt(vx*vx + vy*vy)
  469. aa=max(aa,1.d-20)
  470. vx=vx/aa
  471. vy=vy/aa
  472. * calcul de la distance au premier pt
  473. rr2=(xx-x11)*(xx-x11)+(yy-y11)*(yy-y11)
  474. * calcul des produits scals et du produit vectoiel
  475. psca1= vx*(xx-x11)+vy*(yy-y11)
  476. psca2= vx*(xx-x22)+vy*(yy-y22)
  477. pvec = vx*(yy-y11)-vy*(xx-x11)
  478. xsi=psca1*psca2
  479. xd= rr2 - psca1*psca1
  480. c write(6,*) 'xsi,xd,pvec,psca1,rr2=',xsi,xd,pvec,psca1,rr2
  481.  
  482. * Cas des pointes (=1er et dernier segments)
  483. if (iseg.eq.1.or.iseg.eq.nseg) then
  484. * on sauve la distance a la 2eme pointe (1er segment)
  485. if (iseg.eq.1) then
  486. rr22 = rr2
  487. psca1b = psca1
  488. endif
  489. * on re-teste ici le critere
  490. if (icle.eq.1) then
  491. if(xx.lt.xmin) goto 101
  492. if(xx.gt.xmax) goto 101
  493. if(yy.lt.ymin) goto 101
  494. if(yy.gt.ymax) goto 101
  495. endif
  496. endif
  497.  
  498. * c'est la premiere fois que l'on voit ce point
  499. c write(6,*) 'xminmax=',xmin,xmax,' yminmax=',ymin,ymax
  500. if (isup(ipt).eq.0) then
  501. isup(ipt)=1
  502. igard=igard+1
  503. c write(6,*) 'isup=',isup(ipt),' --> igard=',igard
  504. endif
  505. ivu=ivu+1
  506.  
  507. * -on est dans le cas ou le point se projette sur le segment
  508. if (xsi.le.0.d0) then
  509. if(xd.lt.0.D0) xd=0.D0
  510. if (iseg.eq.1) then
  511. dmin= xd
  512. dmin11=pvec
  513. dis = sign((sqrt(xd)),pvec)
  514. else
  515. if (xd.le.dmin)then
  516. dmin= xd
  517. dmin11=pvec
  518. dis = sign((sqrt(xd)),pvec)
  519. endif
  520. endif
  521. * -le point ne se projette pas sur le segment
  522. else
  523. * est-il plus proche du premier ou du dernier pt du segment ?
  524. if (psca1.ge.0.D0) then
  525. rr2=(xx-x22)*(xx-x22)+(yy-y22)*(yy-y22)
  526. xd= rr2 - psca2*psca2
  527. endif
  528. if(xd.lt.0.D0) xd=0.D0
  529. * --si 1ere fois qu on traite ce point
  530. if (ivu.eq.1) then
  531. dmin= rr2
  532. dmin11=pvec
  533. dis = sign((sqrt(xd)),pvec)
  534. * --on a deja vu ce point
  535. else
  536. * ---si, avec ce segment, on est franchement + pres
  537. c if (rr2.le.dmin) then
  538. if (rr2.lt.0.999999*dmin) then
  539. dmin= rr2
  540. dmin11=pvec
  541. if (iseg.eq.nseg) then
  542. dis = sign((sqrt(xd)),pvec)
  543. else
  544. dis = sign((sqrt(rr2)),pvec)
  545. endif
  546. * ---si on est a quasi-egalite avec un autre segment de fissure
  547. elseif (rr2.lt.1.000001*dmin) then
  548. ss1=pvec*dmin
  549. if (ss1.le.0.d0) then
  550. * bp 20.06.2011 : egalite entre 2 concurrents mais signes opposés !
  551. if (abs(pvec).gt.abs(dmin11)) then
  552. dmin= rr2
  553. dmin11=pvec
  554. if (iseg.eq.nseg) then
  555. dis = sign((sqrt(xd)),pvec)
  556. else
  557. dis = sign((sqrt(rr2)),pvec)
  558. endif
  559. endif
  560. endif
  561. endif
  562. * ---fin de la disctinction selon valeur de rr2
  563. endif
  564. * --fin de la distinction vu ou pas
  565. endif
  566. * -fin de la distinction point se projette ou pas sur le segment
  567. c if(ia.eq.1799)
  568. c write(6,*) ipt,isup(ipt),psca1,xd,rr2,pvec,dis,dmin
  569.  
  570. c enddo
  571. 101 continue
  572. * fin de boucle sur segment --------------------------------------
  573.  
  574. * -cas ou on stocke phi (ex-dernier segment)
  575. if (isup(ipt).ne.0) then
  576. c write(6,*) ipt,'stockage de phi =',dis
  577. * stockage de phi
  578. mpova2.vpocha(ipt,1)=dis
  579. * --puis on calcule les valeurs de psi ...
  580. if (ideux.ge.1) then
  581. * ...relatif a la 1ere pointe de fissure (dernier segment)
  582. rr2 = (xx-x22)*(xx-x22)+(yy-y22)*(yy-y22)
  583. if (psca2.ge.0.) then
  584. mpova1.vpocha(ipt,1)=psca2
  585. else
  586. yd = (rr2 - dis*dis)
  587. if(yd.lt.0.D0) yd=0.D0
  588. mpova1.vpocha(ipt,1)= -1.*(yd**0.5)
  589. endif
  590. * ---...relatif a la 2eme pointe de fissure (=le 1er segment)
  591. if (ip2.ne.0) then
  592. if (psca1b.le.0.) then
  593. psi2= -1.*psca1b
  594. else
  595. psi2= (rr22 - dis*dis)
  596. if (psi2.lt.0.D0) then
  597. psi2= 0.D0
  598. else
  599. psi2= -1.*(psi2**0.5)
  600. endif
  601. endif
  602. * et on choisit le psi le plus proche
  603. cbp2013 if(psi2.gt.(mpova1.vpocha(ipt,1)))
  604. cbp2013 & mpova1.vpocha(ipt,1) = psi2
  605. if(rr2.gt.rr22) mpova1.vpocha(ipt,1) = psi2
  606. endif
  607. * ---fin du cas 2 pointes
  608. endif
  609. * --fin du cas ideux.ge.1 (calcul de psi)
  610. endif
  611. * -fin du cas ou on stocke phi (ex-dernier segment)
  612.  
  613. 100 continue
  614. * fin de boucle sur point ******************************************
  615.  
  616. c endif
  617. c ***** fin de la distinction icle= 1 / 0 *******************************
  618.  
  619.  
  620. ***********************************************************************
  621. * FIN DU TRAVAIL
  622. ***********************************************************************
  623.  
  624.  
  625. * faut-il ajuster les champs?
  626. iaju=0
  627. do ia=1,isup(/1)
  628. * on inverse isup
  629. if (isup(ia).eq.0) then
  630. iaju=iaju+1
  631. isup(ia)=1
  632. else
  633. isup(ia)=0
  634. endif
  635. enddo
  636. if (iaju+igard.ne.isup(/1)) then
  637. write(ioimp,*) 'psiphi : tout va mal '
  638. $ ,iaju,'+',igard,'ne',isup(/1)
  639. call erreur(26)
  640. return
  641. endif
  642.  
  643. * ajustement des champs
  644. if (iaju.ne.0) then
  645. nbelem=ipt1.num(/2)- iaju
  646. nbnn=1
  647. nbref=0
  648. nbsous=0
  649. segini ipt2
  650. ipt2.itypel=1
  651. nc=1
  652. n=nbelem
  653. segini mpova5
  654. *>>>bp&yt: ajout ligne suivante
  655. if(ideux.ge.1) segini mpova4
  656. if(ideux.ge.2) segini mpova6
  657. iel=0
  658. do i=1,ipt1.num(/2)
  659. if (isup(i).eq.0) then
  660. iel=iel+1
  661. ipt2.num(1,iel)=ipt1.num(1,i)
  662. mpova5.vpocha(iel,1)=mpova2.vpocha(i,1)
  663. *>>>bp&yt: ajout ligne suivante
  664. * on met le meme maillage pour les 2 level set (car toujours utilisées ensemble)
  665. if(ideux.ge.1) mpova4.vpocha(iel,1)=mpova1.vpocha(i,1)
  666. if(ideux.ge.2) mpova6.vpocha(iel,1)=mpova3.vpocha(i,1)
  667. endif
  668. enddo
  669. if (ideux.ge.2) then
  670. c segact msoup3*mod
  671. msoup3.ipoval=mpova6
  672. msoup3.igeoc=ipt2
  673. segsup,mpova3
  674. endif
  675. if (ideux.ge.1) then
  676. c segact msoup1*mod
  677. msoup1.ipoval=mpova4
  678. msoup1.igeoc=ipt2
  679. segsup,mpova1
  680. endif
  681. c segact msoup2*mod
  682. msoup2.ipoval=mpova5
  683. msoup2.igeoc=ipt2
  684. segsup,mpova2
  685. endif
  686. segsup ttrav
  687. segsup isup
  688. c if ((idim.eq.2).and.(icle.eq.1)) then
  689. c if (icle.eq.1) then
  690. c segsup xdis,xdmin,xdmin1,rrlim
  691. c if(ideux.ge.1) segsup xdmin2,xdmin3
  692. c if(ip2.ne.0) segsup xrr2b,xpsc1b
  693. c endif
  694.  
  695.  
  696. end
  697.  
  698.  
  699.  

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