Télécharger psip2d.eso

Retour à la liste

Numérotation des lignes :

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

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