Télécharger la plaquette Cast3M


fr en

Télécharger psip2d.eso

Retour à la liste

Numérotation des lignes : 

psip2d

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

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