Télécharger la plaquette Cast3M


fr en

Télécharger psip3d.eso

Retour à la liste

Numérotation des lignes : 

psip3d

  1. C PSIP3D    SOURCE    CB215821  20/11/25    13:37:56     10792          
  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         BP, 2018-06 correction d'une petite erreur + //isation
  14. c
  15. ***********************************************************************
  16.       IMPLICIT REAL*8(A-H,O-Z)
  17.       IMPLICIT INTEGER(I-N)
  18.  
  19. -INC PPARAM
  20. -INC CCOPTIO
  21. -INC SMCOORD
  22. -INC SMELEME
  23. -INC SMCHPOI
  24. -INC CCREEL
  25. -INC CCASSIS
  26. C     Declaration du COMMON pour le travail en parallele
  27.       COMMON/psiabc/IPARAL
  28.  
  29.       segment icpr(nbpts)
  30.       segment mxicpr
  31.         real*8 xicpr(nbpts)
  32.       endsegment 
  33.       segment isup(npt)
  34.       segment sfis
  35.          integer ifis(nbelfi,2)
  36.       endsegment
  37.  
  38. c     MTRI3 : toutes les infos pour decrire un ensemble de TRI3  
  39. c     noeuds, coordonnees, vecteur unitaire des aretes, normale, aire
  40.       SEGMENT MTRI3
  41.         INTEGER KNODE(3,NTRI3)
  42.         REAL*8  XNODE(3,NTRI3),YNODE(3,NTRI3),ZNODE(3,NTRI3)
  43. c         REAL*8  XMIMA(NTRI3),YMIMA(NTRI3),ZMIMA(NTRI3)
  44.         REAL*8  VARE12(4,NTRI3),VARE23(3,NTRI3),VARE31(3,NTRI3)
  45.         REAL*8  VNORM(3,NTRI3),SURFAC(NTRI3)
  46.       ENDSEGMENT
  47.  
  48. c     SPARAL : pour la parallelisation
  49. C     NBTHRD  : nombre de threads demandes
  50. c     KPT3    : ...
  51.       SEGMENT SPARAL
  52.         INTEGER NBTHRD
  53.         INTEGER IERROR(NBTHR)
  54.         INTEGER KPT3,KPOVA1,KPOVA2,KPOVA3
  55.         INTEGER KTRI3,KXICPR,KSUP,KFIS
  56.         INTEGER KCLE,KDEUX,KMIL5
  57.         REAL*8  XCRIT1
  58.       ENDSEGMENT
  59.  
  60.       DATA EPSI/1.D-5/
  61.       INTEGER ooperm(7,2)
  62.       DATA ooperm/ 2, 4, 6, 1, 3, 5, 7,
  63.      $             4, 1, 5, 2, 6, 3, 7/
  64.  
  65.       EXTERNAL PSI3Di
  66.       LOGICAL BTHRD
  67.  
  68.  
  69. ***********************************************************************
  70. *     Initialisation, Recup, Activation
  71. ***********************************************************************
  72.  
  73.       idim1=idim+1
  74. c     idbug=numero du noeud pour lequel on souhaite afficher le calcul
  75.       idbug = 0
  76. c       idbug = 16346
  77.  
  78. c     des chpoints (les msoup sont deja actifs)
  79.       mpova1=0
  80.       mpova2=0
  81.       mpova3=0
  82.       if(msoup2.gt.0) mpova2=msoup2.ipoval
  83.       if(msoup1.gt.0) mpova1=msoup1.ipoval
  84.       if(msoup3.gt.0) mpova3=msoup3.ipoval
  85. c     activation des maillages             / desactivation
  86. c     melmai : fait dans psiphi            / ligne 1325 ou 1329
  87. c     melfis : fait ligne 64               / ligne 1274
  88. c     melfro : fait (si ideux>=1) ligne 78 / ligne 356
  89.       ipt3=melmai
  90. c        segact,ipt3
  91.       npt=ipt3.num(/2)
  92.  
  93. c     variables parfois testées
  94.       imil5=0
  95.  
  96.  
  97. ***********************************************************************
  98. *     Preparation du maillage de la fissure
  99. ***********************************************************************
  100.         meleme=melfis
  101.         segact meleme
  102. *       transfo du front en TRI3
  103.         if(itypel.ne.4)  call change(melfis,4)
  104.         meleme=melfis
  105.         segact meleme
  106. c rem: inutile car change laisse melfis actif
  107.         nseg=num(/2)
  108.  
  109. *------ Cas ou l'on calcule phi ET psi --------------------------------
  110.         if (ideux.ge.1) then
  111. * on rearrange le melfis de facon à prevenir dans ifis(iel,1) quels
  112. * element touchent le front de fissure. Ceux ci sont ordonnées de facon
  113. * que le front de fissure soit du noeud 1 vers le noeud 2 de l'element
  114.           ipt2=melfro
  115.           segact ipt2
  116.           segini icpr
  117.           nbelem2 = ipt2.num(/2)
  118.           nbnode2 = ipt2.num(/1)
  119. * petite numerotation locale
  120. * icpr(noeud)>0 si le noeud appartient au front de fissure
  121. *            =numero du 1er element de melfro ou on a vu ce noeud
  122. * + on en profite pour calculer l'abscisse curviligne du front si ideux.eq.2
  123. * icpr(noeud)<0 si le noeud appartient au reste du contour de la fissure
  124. *            =-numero du 1er element de melcon ou on a vu ce noeud
  125.           if (ideux.ge.2) then
  126.             segini,mxicpr
  127.             xtau=0.d0
  128.             ic=ipt2.num(1,1)
  129.             x11= xcoor((ic-1)*idim1 +1)
  130.             y11= xcoor((ic-1)*idim1 +2)
  131.             z11= xcoor((ic-1)*idim1 +3)
  132.             do ia=1,nbelem2
  133.               do ib=1,nbnode2
  134.                 ic=ipt2.num(ib,ia)
  135.                 if (icpr(ic).eq.0) then
  136.                   icpr(ic)=ia
  137.                   x22= xcoor((ic-1)*idim1 +1)
  138.                   y22= xcoor((ic-1)*idim1 +2)
  139.                   z22= xcoor((ic-1)*idim1 +3)
  140.                   dxtau= sqrt( (x22-x11)*(x22-x11)+(y22-y11)*(y22-y11)
  141.      $           +(z22-z11)*(z22-z11) )
  142.                   xtau=xtau+dxtau
  143.                   xicpr(ic)=xtau
  144.                   x11=x22
  145.                   y11=y22
  146.                   z11=z22
  147. c           if(idebug.ge.1)
  148. c      $    write(6,*) 'ia,ib,ic,xicpr(ic)=',ia,ib,ic,xicpr(ic)
  149.                 endif
  150.               enddo
  151.             enddo
  152.           else
  153.             do ia=1,nbelem2
  154.               do ib=1,nbnode2
  155.                 ic=ipt2.num(ib,ia)
  156.                 if (icpr(ic).eq.0) then
  157.                   icpr(ic)=ia
  158.                 endif
  159.               enddo
  160.             enddo
  161.           endif
  162. * bp (23/02/2012) : calcul du contour et fin de remplissage de icpr
  163.          CALL ECROBJ('MAILLAGE',melfis)
  164.          CALL PRCONT
  165.          CALL LIROBJ('MAILLAGE',melcon,1,IRETOU)
  166.          IF (IERR.NE.0) RETURN
  167.          meleme = melfis
  168.          IPT6   = melcon
  169.          SEGACT,meleme,IPT6
  170.           nbelCON = IPT6.num(/2)
  171.           nbnoCON = IPT6.num(/1)
  172.             do ia=1,nbelCON
  173.               do ib=1,nbnoCON
  174.                 ic=IPT6.num(ib,ia)
  175. c               noeud qui n'a pas ete vu dans melfro et qui appartient au contour
  176.                 if(icpr(ic).eq.0)  icpr(ic)=-ia
  177.               enddo
  178.             enddo
  179. * nombre de segments appartenant au contour sans le front
  180.          nbelem6 = nbelCON - nbelem2
  181. * remplissage effectif de ifis contenu dans le segment sfis
  182.          i2firs = ipt2.num(1,1)
  183.          i2last = ipt2.num(nbnode2,nbelem2)
  184.          ivui2f = 0
  185.          ivui2l = 0
  186. c          write(6,*) 'i2firs,i2last=',i2firs,i2last
  187.          nbelfi=nseg
  188.          nbelem=nseg
  189.          nbnn=3
  190.          nbsous=0
  191.          nbref=0
  192. c          segini,ipt5=meleme
  193. * bp 20/12/2011 : on propose de remplir ipt5 avec :
  194. c + des elements du front (et du contour) d abord,
  195. c + des elements avec 1 point appartenant au front (et au contour) ensuite,
  196. c + les autres elements enfin.
  197. c bp 06/03/2012 : on propose de remplir ifis avec :
  198. c ifis(iseg,1) =2 si le segment [1-2] du iseg^ieme triangle appartient
  199. c au front de fissure, =3 si en plus c'est l extrémité du front,
  200. c =1 si l un des segments appartient au contour
  201. c =-1 si l'un des noeuds appartient au front
  202. c ifis(iseg,2) =+1+2+4 si les segments [1-2] [2-3] [3-1] du iseg^ieme triangle appartiennent au contour
  203.          segini,ipt5
  204.          ideb5 = 0
  205.          imil5 = nbelem2
  206.          ifin5 = nseg+1
  207.          segini sfis
  208.  
  209. *--------boucle sur les elements de la fissure ----------------------
  210.          do iel=1,nseg
  211. c             write(*,*) 'boucle sur elements de fissure',iel,'/',nseg
  212.  
  213. *     io est incrémenté de +1,+2,+4 si node 1,2,3 appartient au front
  214. *     ioo est incrémenté de +1,+2,+4 si node 1,2,3 appartient au contour
  215.             io=0
  216.             ioo=0
  217. *        ---boucle sur les noeuds des elements de la fissure ----------
  218.             do ib=1,num(/1)
  219.               ia=ib
  220.               if(ib.eq.3) ia = 4
  221.               ic = num(ib,iel)
  222.               if(icpr(ic).gt.0) io=io+ia
  223.               if(icpr(ic).lt.0) ioo=ioo+ia
  224. *             cas particulier d un segment avec 1 point appartenant au contour et un autre
  225.               if(ic.eq.i2firs.or.ic.eq.i2last) ioo=ioo+ia
  226.             enddo
  227. *        ---fin de boucle sur les noeuds des elements de la fissure ---
  228.  
  229. *     remplissage de ipt5(i5,...) et ifis(i5,1)
  230.             i5=0
  231.             iperm=0
  232. *        ---cas ou 1 segment appartient au front
  233.             if (io.eq.3.or.io.gt.4) then
  234.               ideb5=ideb5+1
  235.               i5=ideb5
  236. *             (eventuelle) permutation des noeuds pour avoir segment 1-2 = front
  237.               if (io.eq.3) then
  238.                   ipt5.num(1,i5)=num(1,iel)
  239.                   ipt5.num(2,i5)=num(2,iel)
  240.                   ipt5.num(3,i5)=num(3,iel)
  241.               elseif(io.eq.5) then
  242.                   ipt5.num(1,i5)=num(3,iel)
  243.                   ipt5.num(2,i5)=num(1,iel)
  244.                   ipt5.num(3,i5)=num(2,iel)
  245.                   iperm=1
  246.               elseif(io.eq.6) then
  247.                   ipt5.num(1,i5)=num(2,iel)
  248.                   ipt5.num(2,i5)=num(3,iel)
  249.                   ipt5.num(3,i5)=num(1,iel)
  250.                   iperm=2
  251.               elseif(io.eq.7) then
  252.                   write(ioimp,*) 'ERREUR: '
  253.      $           ,'Plusieurs cotes d un meme triangle du maillage de '
  254.      $           ,'la fissure converti en TRI3 appartiennent au front!'
  255.                   return
  256.               endif
  257. *             ifis( ,1)=3 si c'est une extremité
  258.               if ((ipt5.num(1,i5)).eq.i2firs) then
  259.                 ifis(i5,1)=3
  260.                 ivui2f = ivui2f+1
  261.               elseif((ipt5.num(2,i5)).eq.i2last) then
  262.                 ifis(i5,1)=3
  263.                 ivui2l = ivui2l+1
  264.               else
  265. *             ifis( ,1)=2 si c'est un tri3 du front avec des voisins
  266.                 ifis(i5,1)=2
  267.               endif
  268. *        ---cas ou 1 (ou 2) segment(s) appartient au contour
  269.             elseif(ioo.eq.3.or.ioo.gt.4) then
  270.               imil5=imil5+1
  271.               i5=imil5
  272. *             pas de permutation des noeuds
  273.               ipt5.num(1,i5)=num(1,iel)
  274.               ipt5.num(2,i5)=num(2,iel)
  275.               ipt5.num(3,i5)=num(3,iel)
  276. *             ifis( ,1)=1 si un segment appartient au contour
  277.               ifis(i5,1)=1
  278. *        ---cas ou seulement 1 point appartient au front
  279.             elseif(io.eq.1.or.io.eq.2.or.io.eq.4) then
  280.               ifin5=ifin5-1
  281.               i5=ifin5
  282. *             pas de permutation des noeuds
  283.               ipt5.num(1,i5)=num(1,iel)
  284.               ipt5.num(2,i5)=num(2,iel)
  285.               ipt5.num(3,i5)=num(3,iel)
  286. *             ifis( ,1)=-1 si 1 seul point du tri3 appartient au front
  287.               ifis(i5,1)=-1
  288. *        ---autres cas
  289.             else
  290.               ifin5=ifin5-1
  291.               i5=ifin5
  292. *             pas de permutation des noeuds
  293.               ipt5.num(1,i5)=num(1,iel)
  294.               ipt5.num(2,i5)=num(2,iel)
  295.               ipt5.num(3,i5)=num(3,iel)
  296.             endif
  297.  
  298. *     remplissage de ifis(i5,2)
  299. *        ---uniquement si 1 (ou 2) point(s) appartient au contour
  300.             if(ioo.eq.3.or.ioo.gt.4) then
  301. *             attention a l eventuelle permutation effectuee
  302.               if(iperm.ne.0)  ioo=ooperm(ioo,iperm)
  303. *            -ioo=3 => noeuds 1-2 = segment 1 appartient au contour
  304.               if (ioo.eq.3) then
  305.                  ifis(i5,2)=1
  306. *            -ioo=5 => noeuds 3-1 = segment 3 appartient au contour
  307.               elseif(ioo.eq.5) then
  308.                  ifis(i5,2)=4
  309. *            -ioo=6 => noeuds 2-3 = segment 2 appartient au contour
  310.               elseif(ioo.eq.6) then
  311.                  ifis(i5,2)=2
  312. *            -cas ou 3 points appartiennent au contour => travail en +
  313. *             pour determiner les 2 segments dans le contour 
  314.               elseif(ioo.eq.7) then
  315. c             write(6,*)'tri3 avec 2 segments dans le contour!',iel,i5
  316. * rem: si les 3 points appartiennent au contour, on n'a pas de segments
  317. *      appartenant au front
  318.                  ii1 = ipt5.num(1,i5)
  319.                  ii2 = ipt5.num(2,i5)
  320.                  ii3 = ipt5.num(3,i5)
  321. c il faut trouver quels sont les 2 segments appartenant au contour !
  322. c                  ia1 = -icpr(ii1)
  323. c                  ia2 = -icpr(ii2)
  324. c                  ia3 = -icpr(ii3)
  325. c dans de tres rares cas, on peut avoir des points appartenant aussi au front (i2firs et i2last)
  326.                  ia1 = abs(icpr(ii1))
  327.                  ia2 = abs(icpr(ii2))
  328.                  ia3 = abs(icpr(ii3))
  329.                  ia1min = min(ia1,ia2)
  330.                  ia1min = min(ia1min,ia3)
  331.                  ia1seg = 0
  332.                  do 11 ia1=ia1min,nbelCON
  333.                    ifis1 = 0
  334.                    do ib=1,nbnoCON
  335.                      if(ii1.eq.IPT6.num(ib,ia1))  ifis1=ifis1+1
  336.                      if(ii2.eq.IPT6.num(ib,ia1))  ifis1=ifis1+2
  337.                      if(ii3.eq.IPT6.num(ib,ia1))  ifis1=ifis1+4
  338.                    enddo
  339. c                  si on n'a trouve que 1 noeud dans le ia1 ieme SEG2 du contour, on continue a chercher
  340.                    if(ifis1.le.2) goto 11
  341.                    if(ifis1.eq.4) goto 11
  342. c                  on a trouve les 2 noeuds dans un SEG2 du contour
  343.                    if(ifis1.eq.3) ifis(i5,2)=ifis(i5,2)+1
  344.                    if(ifis1.eq.6) ifis(i5,2)=ifis(i5,2)+2
  345.                    if(ifis1.eq.5) ifis(i5,2)=ifis(i5,2)+4
  346.                    if(ifis1.ge.7) then
  347.              write(ioimp,*) 'tri3 avec 3 segments dans le contour !'
  348.                      call erreur(26)
  349.                    endif
  350.                    if(ifis1.eq.3.or.ifis1.ge.5)  ia1seg=ia1seg+1
  351.                    if(ia1seg.eq.2) goto 12
  352.   11               continue
  353. c on n a pas trouvé les 2 segments appartenant au contour
  354.                    write(ioimp,*)'on n a pas trouvé les 2 segments',
  355.      $                   ' du tri3 devant appartenant au contour !'
  356.                    call erreur(26)
  357.                    return
  358. c on a trouvé les 2 segments appartenant au contour
  359.   12               continue  
  360.                  if(ifis(i5,2).lt.3.or.ifis(i5,2).ge.7) then
  361.              write(ioimp,*) 'impossible de detecter les 2 segments ',
  362.      $       'appartenant au contour avec les noeuds ',ii1,ii2,ii3
  363.              write(ioimp,*) 'ifis(ideb5,2)=',ifis(ideb5,2)
  364.              write(ioimp,*) 'IPT6 pour les segments :'
  365.              write(ioimp,*) ' ia1',ia1,(IPT6.num(iou,ia1),iou=1,2)
  366.              write(ioimp,*) ' ia2',ia2,(IPT6.num(iou,ia2),iou=1,2)
  367.              write(ioimp,*) ' ia3',ia3,(IPT6.num(iou,ia3),iou=1,2)
  368.                    call erreur(26)
  369.                    return
  370.                  endif
  371.               endif
  372. c            -fin du cas elseif(ioo.eq.7) then  
  373.             endif
  374. c         --fin du cas if(ioo.eq.3.or.ioo.gt.4) then    
  375.  
  376.  
  377.  
  378.          enddo
  379. *        fin de boucle sur les elements de la fissure -----------------
  380.  
  381. c          if ((ifin5-ideb5).ne.1) then
  382.          if (ideb5.ne.nbelem2.or.(ifin5-imil5).ne.1) then
  383.        write(ioimp,*) 'pb avec le rearrangement de elements de fissure'
  384.        write(ioimp,*) ideb5,nbelem2,imil5,ifin5
  385.            call erreur(26)
  386.            return
  387.          endif
  388. *        si le front n est pas ferme, il faut avoir detecte ses 2 extremites
  389. *        sinon il est probable que la ligne du front melfro ne soit pas
  390. *        orientée comme le bord du maillage melfis
  391.          if ((ivui2f.ne.1).and.(ivui2l.ne.1)) then
  392.             write(ioimp,*) 'ERREUR: Les extremites du front (noeuds',
  393.      $      i2firs,i2last,') ont été detectees ',ivui2f, 'et ',ivui2l
  394.      $      ,' fois au lieu de 1 fois attendue'
  395.             write(ioimp,*) 'Verifiez la coherence avec l orientation'
  396.      $      ,'du maillage de la fissure fourni'
  397.      $      ,' et INVErser le sens du front le cas échéant.'
  398. c             write(6,*) 'front de fissure='
  399. c             CALL ECMAIL(ipt2,JENT2)
  400. c             write(6,*) 'surface de fissure='
  401. c             CALL ECMAIL(ipt5,JENT5)
  402.             call erreur(26)
  403.             return
  404.          endif
  405. c          segdes,meleme,IPT6
  406.          segdes,meleme,IPT6,ipt2
  407.          meleme=ipt5
  408.          melfis=meleme
  409.          segsup icpr
  410.  
  411.         endif
  412. *------ fin du Cas ou l'on calcule phi ET psi --------------------------
  413.  
  414. * sont actifs: melmai(=struct changée en poi1), ipt5=meleme=melfis,
  415. * sont inactifs: ipt2=melfro
  416.  
  417.  
  418.  
  419. ***********************************************************************
  420. * debut du travail suivant valeur de icle
  421.  
  422.       if (xcrit.eq.0.d0) then
  423.         icle=0
  424.         xcrit=1.D10
  425.       else
  426.         icle=1
  427. *       ajout d'un epsilon 1.D-5 de rattrapage
  428.         xcrit = 1.00001D0 * xcrit
  429.       endif
  430.  
  431.       melfis=meleme
  432.       melpoi=melmai
  433.  
  434. ***********************************************************************
  435. *     IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT
  436. ***********************************************************************
  437.  
  438. *  pour chaque segment  de la fissure (melfis 2eme maillage lu)
  439. *  on calcule le champ phi et psi des noeuds autour de l'element
  440.         if(iimpi.ge.1) then
  441. c         call gibtem(xkt)
  442.         write(ioimp,*)'debut travail'
  443.         endif
  444. *
  445.        meleme=melfis
  446. c        segact,meleme
  447.        ipt3=melpoi
  448. c        segact,ipt3
  449.        npt=ipt3.num(/2)
  450. * attention on initialise isup avec les valeurs inverses,il faudra
  451. * le retourner
  452. c c        segini isup,isupfi,xdis,xdmin,xdmin1,rrlim
  453. c        segini isup,isupfi,xdis,xdmin,xdmin1
  454.        segini isup
  455. c        xdis,xdmin,xdmin1
  456. c        if(ideux.ge.1) segini xdmin2,xdmin3
  457.  
  458.        NTRI3=nseg
  459.        SEGINI,MTRI3
  460.        call cpu_time(start_time)
  461. *=====================================================================
  462. *      Boucle sur les elements de la fissure =========================
  463.  
  464.        do 444 iseg=1,nseg
  465.  
  466. c          coordonnees du tri3
  467.            ia=num(1,iseg)
  468.            x11= xcoor((ia-1)*idim1 +1)
  469.            y11= xcoor((ia-1)*idim1 +2)
  470.            z11= xcoor((ia-1)*idim1 +3)
  471.            ib=num(2,iseg)
  472.            x22= xcoor((ib-1)*idim1 +1)
  473.            y22= xcoor((ib-1)*idim1 +2)
  474.            z22= xcoor((ib-1)*idim1 +3)
  475.            ic=num(3,iseg)
  476.            x33= xcoor((ic-1)*idim1 +1)
  477.            y33= xcoor((ic-1)*idim1 +2)
  478.            z33= xcoor((ic-1)*idim1 +3)
  479. c            if (x11.ge.x22) then
  480. c             xxmi=min(x22,x33) - xcrit
  481. c             xxma=max(x11,x33) + xcrit
  482. c            else
  483. c             xxmi=min(x11,x33) - xcrit
  484. c             xxma=max(x22,x33) + xcrit
  485. c            endif
  486. c            if (y11.ge.y22) then
  487. c             yymi=min(y22,y33) - xcrit
  488. c             yyma=max(y11,y33) + xcrit
  489. c            else
  490. c             yymi=min(y11,y33) - xcrit
  491. c             yyma=max(y22,y33) + xcrit
  492. c            endif
  493. c            if (z11.ge.z22) then
  494. c             zzmi=min(z22,z33) - xcrit
  495. c             zzma=max(z11,z33) + xcrit
  496. c            else
  497. c             zzmi=min(z11,z33) - xcrit
  498. c             zzma=max(z22,z33) + xcrit
  499. c            endif
  500.  
  501. c          vecteurs des aretes
  502.            v12x=x22-x11
  503.            v12y=y22-y11
  504.            v12z=z22-z11
  505.            v23x=x33-x22
  506.            v23y=y33-y22
  507.            v23z=z33-z22
  508.            v31x=x11-x33
  509.            v31y=y11-y33
  510.            v31z=z11-z33
  511.  
  512. c          vecteur normal
  513.            vx=v12y*(z33-z11)-v12z*(y33-y11)
  514.            vy=v12z*(x33-x11)-v12x*(z33-z11)
  515.            vz=v12x*(y33-y11)-v12y*(x33-x11)
  516.            bb=vx*vx + vy*vy+  vz*vz
  517.            aa= sqrt(bb)
  518.            if (aa.le.XPETIT) then
  519.              write(ioimp,*) 'element fissure de longueur nulle ',aa
  520.              write(ioimp,*) 'tri3.',iseg,'.',ia,'=',x11,y11,z11
  521.              write(ioimp,*) 'tri3.',iseg,'.',ib,'=',x22,y22,z22
  522.              write(ioimp,*) 'tri3.',iseg,'.',ic,'=',x33,y33,z33
  523.              call erreur(26)
  524.              return
  525.            endif
  526.            sur = bb/2.d0
  527.            vx=vx/aa
  528.            vy=vy/aa
  529.            vz=vz/aa
  530.            vxyz = max(abs(vx),abs(vy))
  531.            vxyz = max(vxyz,abs(vz))
  532.            if (vxyz.le.XPETIT) then
  533.            write(ioimp,*)'element fissure de longueur nulle',aa,vx,vy,vz
  534.              call erreur(26)
  535.              return
  536.            endif
  537. c            xsur = 1.D-8*sur
  538.  
  539. c          vecteurs des aretes (suite)
  540.            v12lo= sqrt(v12x*v12x+v12y*v12y+v12z*v12z)
  541.            v23lo= sqrt(v23x*v23x+v23y*v23y+v23z*v23z)
  542.            v31lo= sqrt(v31x*v31x+v31y*v31y+v31z*v31z)
  543.            v12x= v12x/v12lo
  544.            v12y= v12y/v12lo
  545.            v12z= v12z/v12lo
  546.            v23x= v23x/v23lo
  547.            v23y= v23y/v23lo
  548.            v23z= v23z/v23lo
  549.            v31x= v31x/v31lo
  550.            v31y= v31y/v31lo
  551.            v31z= v31z/v31lo
  552.  
  553. c          REMPLISSAGE DU SEGMENT MTRI3  
  554.            MTRI3.KNODE(1,iseg)=ia
  555.            MTRI3.KNODE(2,iseg)=ib
  556.            MTRI3.KNODE(3,iseg)=ic
  557.            MTRI3.XNODE(1,iseg)=x11
  558.            MTRI3.XNODE(2,iseg)=x22
  559.            MTRI3.XNODE(3,iseg)=x33
  560.            MTRI3.YNODE(1,iseg)=y11
  561.            MTRI3.YNODE(2,iseg)=y22
  562.            MTRI3.YNODE(3,iseg)=y33
  563.            MTRI3.ZNODE(1,iseg)=z11
  564.            MTRI3.ZNODE(2,iseg)=z22
  565.            MTRI3.ZNODE(3,iseg)=z33
  566.            MTRI3.VARE12(1,iseg)=v12x
  567.            MTRI3.VARE12(2,iseg)=v12y
  568.            MTRI3.VARE12(3,iseg)=v12z
  569.            MTRI3.VARE12(4,iseg)=v12lo
  570.            MTRI3.VARE23(1,iseg)=v23x
  571.            MTRI3.VARE23(2,iseg)=v23y
  572.            MTRI3.VARE23(3,iseg)=v23z
  573.            MTRI3.VARE31(1,iseg)=v31x
  574.            MTRI3.VARE31(2,iseg)=v31y
  575.            MTRI3.VARE31(3,iseg)=v31z
  576.            MTRI3.VNORM(1,iseg)=vx
  577.            MTRI3.VNORM(2,iseg)=vy
  578.            MTRI3.VNORM(3,iseg)=vz
  579.            MTRI3.SURFAC(iseg)=sur
  580.  
  581.  444   continue
  582. *=====================================================================
  583.  
  584.  
  585. *=====================================================================
  586. *     2EME PHASE DU TRAVAIL (DOUBLE BOUCLE)
  587. *=====================================================================
  588.  
  589. C     FAUT-IL PASSER EN // ?
  590. c     estimation grossiere 
  591. c     car on n'a trouve que des cas ou la // est benefique
  592.       if (nseg*npt.lt.100) then
  593.         NBTHR = 1
  594.         BTHRD=.false.
  595.       else
  596.         NBTHR = NBTHRS      
  597.         BTHRD = .TRUE.
  598.         CALL THREADII
  599.       endif
  600.  
  601. C     CREATION ET REMPLISSAGE DU SEGMENT POUR LA //iSATION
  602.       SEGINI,SPARAL
  603.       SPARAL.NBTHRD = NBTHR
  604.       SPARAL.KPT3   = IPT3
  605.       SPARAL.KTRI3  = MTRI3
  606.       SPARAL.KCLE   = icle
  607.       SPARAL.KDEUX  = ideux
  608.       SPARAL.XCRIT1 = xcrit
  609.       SPARAL.KSUP   = isup
  610.       SPARAL.KPOVA1 = mpova1
  611.       SPARAL.KPOVA2 = mpova2
  612.       SPARAL.KPOVA3 = mpova3
  613.       SPARAL.KFIS   = sfis
  614.       SPARAL.KMIL5  = imil5
  615.       SPARAL.Kxicpr = mxicpr
  616.  
  617.       IF (BTHRD) THEN
  618. C       Remplissage du 'COMMON/psiabc'
  619.         IPARAL=SPARAL
  620.         DO ith=2,NBTHR
  621.           CALL THREADID(ith,PSI3Di)
  622.         ENDDO
  623.         CALL PSI3Di(1)
  624.  
  625. C       Attente de la fin de tous les threads en cours de travail
  626.         DO ith=2,NBTHR
  627.           CALL THREADIF(ith)
  628.         ENDDO
  629.  
  630. C       On libere les Threads
  631.         CALL THREADIS
  632.  
  633. C       Verification de l'erreur (Apres liberation des THREADS)
  634.         DO ith=1,NBTHR
  635.           IRETOU=SPARAL.IERROR(ith)
  636.           IF (IRETOU .GT. 0) THEN
  637.             CALL ERREUR(IRETOU)
  638.             RETURN
  639.           ENDIF
  640.         ENDDO
  641.  
  642.       ELSE
  643. C       Appel de la SUBROUTINE qui fait le travail
  644.         CALL PSI3D0(1,SPARAL)
  645. C       Verification de l'erreur
  646.         IRETOU=SPARAL.IERROR(1)
  647.         IF (IRETOU .GT. 0) THEN
  648.           CALL ERREUR(IRETOU)
  649.           RETURN
  650.         ENDIF
  651.  
  652.       ENDIF      
  653.  
  654. ***********************************************************************
  655. *     FIN DU TRAVAIL
  656. ***********************************************************************
  657.  
  658. * on inverse isup
  659. * AVANT : isup(i) =0 veut dire que le point na pas été traité
  660. * APRES : isup(i) =1 veut dire que le point peut etre enleve (=0 sinon)
  661.         iaju=0
  662.         do ia=1,isup(/1)
  663.          if (isup(ia).eq.0) then
  664.             iaju=iaju+1
  665.             isup(ia)=1
  666.          else
  667.             isup(ia)=0
  668.          endif
  669.         enddo
  670. c         if (iaju+igard.ne.isup(/1)) then
  671. c             write(ioimp,*) 'psiphi : tout va mal '
  672. c      $                 ,iaju,'+',igard,'ne',isup(/1)
  673. c             call erreur(26)
  674. c             return
  675. c         endif
  676. c bp : test ci dessus un peu chiant a garder en // et pas tre utile
  677.         ipt1=ipt3
  678.         ipt4=melfis
  679.         segdes ipt4
  680. * iaju est le nombre de points a éliminer
  681. * isup(i) =1 veut dire que le point peut etre enleve (=0 sinon)
  682.  
  683. * faut-il ajuster les champs?
  684.       if (iaju.ne.0) then
  685. *        si oui on recree les bons chpoints
  686.          nbelem=ipt1.num(/2)- iaju
  687.          nbnn=1
  688.          nbref=0
  689.          nbsous=0
  690.          segini ipt2
  691.          ipt2.itypel=1
  692.          nc=1
  693.          n=nbelem
  694.          segini mpova5
  695. *>>>bp&yt: ajout ligne suivante
  696.          if(ideux.ge.1) segini mpova4
  697.          if(ideux.ge.2) segini mpova6
  698.          iel=0
  699.          do i=1,ipt1.num(/2)
  700.           if (isup(i).eq.0) then
  701.             iel=iel+1
  702.             ipt2.num(1,iel)=ipt1.num(1,i)
  703.             mpova5.vpocha(iel,1)=mpova2.vpocha(i,1)
  704. *>>>bp&yt:  ajout ligne suivante
  705. * on met le meme maillage pour les 2 level set (car toujours utilisées ensemble)
  706.            if(ideux.ge.1)  mpova4.vpocha(iel,1)=mpova1.vpocha(i,1)
  707.            if(ideux.ge.2)  mpova6.vpocha(iel,1)=mpova3.vpocha(i,1)
  708.           endif
  709.          enddo
  710.          if (ideux.ge.2) then
  711. c            segact msoup3*mod
  712.            msoup3.ipoval=mpova6
  713.            msoup3.igeoc=ipt2
  714.            segdes msoup3,mpova6
  715.            segsup,mpova3
  716.          endif
  717.          if (ideux.ge.1) then
  718. c            segact msoup1*mod
  719.            msoup1.ipoval=mpova4
  720.            msoup1.igeoc=ipt2
  721.            segdes msoup1,mpova4
  722.            segsup,mpova1
  723.          endif
  724. c          segact msoup2*mod
  725.          msoup2.ipoval=mpova5
  726.          msoup2.igeoc=ipt2
  727.          segdes msoup2,mpova5,ipt2
  728.          segsup,mpova2
  729.          segdes ipt1
  730.       else
  731.          if(ideux.ge.2) segdes,msoup3,mpova3
  732.          if(ideux.ge.1) segdes,msoup1,mpova1
  733.          segdes msoup2,mpova2,ipt1
  734.       endif
  735.       if (ideux.ge.1) then
  736.           segsup,icpr,sfis
  737.           if(ideux.ge.2) segsup,mxicpr
  738.       endif
  739.       if (icle.eq.1) then
  740. c        segsup isup,isupfi,xdis,xdmin,xdmin1,rrlim
  741. c        segsup isup,isupfi,xdis,xdmin,xdmin1
  742.        segsup isup
  743. c        if(ideux.ge.1) segsup xdmin2,xdmin3
  744.       endif
  745.       SEGSUP,MTRI3
  746.  
  747.       RETURN
  748.       END
  749.  
  750.  
  751.  
  752.  
  753.  
  754.  
  755.  
  756.  
  757.  
  758.  
  759.  
  760.  
  761.  

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