Télécharger la plaquette Cast3M


fr en

Télécharger intgeo.eso

Retour à la liste

Numérotation des lignes : 

intgeo

  1. C INTGEO    SOURCE    PV        20/03/24    21:18:17     10554          
  2. C
  3.       SUBROUTINE INTGEO(IOB1,IOB2,IOB3,IOB4,ICLE)
  4. C-----------------------------------------------------------------------
  5. C
  6. C     INTERSECTION GEOMETRIQUE DE MAILLAGE
  7. C
  8. C-----------------------------------------------------------------------
  9. C
  10. C ENTREE :
  11. C     IOB1   MAILLAGE 1 INITIAL
  12. C    (IOB2   MAILLAGE 2 INITIAL)
  13. C
  14. C SORTIE :
  15. C     IOB1   MAILLAGE 1 MODIFIE
  16. C    (IOB2   MAILLAGE 2 MODIFIE)
  17. C     IOB3   MAILLAGE 3 = SELF-INTERSECTION DU MAILLAGE 1
  18. C                         ou INTERSECTION DU MAILLAGE 1 AVEC 2
  19. C    (IOB4   MAILLAGE 3 = INTERSECTION DU MAILLAGE 2 AVEC 1)
  20. C
  21. C-----------------------------------------------------------------------
  22. C
  23. C  BP208322, 28/08/2012 : Creation
  24. C
  25. C  REM   : - self-intersection abandonnée car compliquée : on préfère
  26. C            chercher à créer directement une bonne surface ...
  27. C
  28. C  TO DO : - cas coplanaire à traiter
  29. c          - optimiser
  30. C
  31. C-----------------------------------------------------------------------
  32. c
  33.       IMPLICIT INTEGER(I-N)
  34.       IMPLICIT REAL*8 (A-H,O-Z)
  35. -INC SMCOORD
  36.  
  37. -INC PPARAM
  38. -INC CCOPTIO
  39. -INC CCGEOME
  40. -INC SMELEME
  41. -INC CCREEL
  42.  
  43.       SEGMENT WORKEL
  44.         INTEGER IDEJVU(NBELEM)
  45.         REAL*8  V(3,NBELEM)
  46.       ENDSEGMENT
  47.       POINTEUR WORK1.WORKEL,WORK2.WORKEL
  48.       INTEGER NODE(2)
  49.  
  50. c     NBNMAX = NomBre de Noeuds MAX
  51. c     ITA(i) = numero du i^eme noeuds du triangle A courant
  52. c     XTA(i,k) = k^ieme coordonnee du i^eme noeuds du triangle A courant
  53. c     XTB(i,k) = k^ieme coordonnee du i^eme noeuds du triangle B courant
  54.       PARAMETER (NBNMAX=3)
  55.       INTEGER ITA(NBNMAX),ITB(NBNMAX)
  56.       REAL*8 XTA(NBNMAX,3),XTB(NBNMAX,3)
  57.       REAL*8 XSEG(1*2,3)
  58.       REAL*8 dA1(3),dA2(3),dA3(3)
  59.       REAL*8 dB1(3),dB2(3),dB3(3)
  60.       REAL*8 xP2(3),xP3(3),xQ2(3),xQ3(3)
  61.       LOGICAL SELF
  62.  
  63.       if(iimpi.GE.10) write(ioimp,*)'====== ENTREE dans INTGEO ======'
  64.  
  65. C-----------------------------------------------------------------------
  66. C     RECUP + VERIFICATION DES DONNEES D ENTREE
  67.  
  68. c     SELF intersection ?
  69.       SELF = (IOB2.eq.0)
  70.       if(SELF) then
  71.         write(ioimp,*) 'Self-intersection non traitée'
  72.         MOTERR(1:8)='MAILLAGE'
  73.         call erreur(37)
  74.       endif
  75. c     on ne traite pas aujourd'hui... abandonnée car compliquée !
  76.  
  77. c     option NOVE => IVERI
  78.       IVERI = ICLE
  79.  
  80. c     Critere absolu de proximite entre 2 noeuds
  81.       CALL LIRREE(XELIM,0,IELIM)
  82. c     Critere absolu de proximite d'un plan (detection du cas coplanaire)
  83.       CALL LIRREE(XPLAN,0,IPLAN)
  84.       if(IPLAN.eq.0) XPLAN=XZPREC**(0.66)
  85.       if(IELIM.ne.0) then
  86.         XELIM2=2.d0*XELIM
  87.         XPLAN2=0.01d0*XELIM
  88.         if(XPLAN2.lt.XPLAN) then
  89.           XPLAN=XPLAN2
  90.         endif
  91.       endif
  92.       XPLANn=-1.d0*XPLAN
  93. c     si XELIM non fourni on calcule avec la taille du 1er element tq :
  94. c     XZPREC << XPLAN << XELIM << taille des elements
  95.  
  96. C     Operateur disponible en dimension 3
  97.       IF (IDIM.NE.3) THEN
  98.         INTERR(1)=IDIM
  99.         CALL ERREUR(709)
  100.         RETURN
  101.       ENDIF
  102.       idim1 = idim+1
  103.       segact,MCOORD*mod
  104.  
  105. C     Recup
  106.       IPT1 = IOB1
  107.       IPT2 = IOB2
  108.       IPT3 = 0
  109.       IPT4 = 0
  110.  
  111. C    +VERIFICATION DE IPT1
  112. C     un maillage non complexe
  113.       segact,IPT1
  114.       NBSOU1=IPT1.LISOUS(/1)
  115.       if (NBSOU1.gt.0) then
  116.         call erreur(25)
  117.         return
  118.       endif
  119. C     maillage de TRI3 uniquement
  120.       ITYPE1=IPT1.ITYPEL
  121.       if(ITYPE1.ne.4) then
  122.         call erreur(16)
  123.         return
  124.       endif
  125.       NBNN1  = IPT1.NUM(/1)
  126.       NBELE1 = IPT1.NUM(/2)
  127.  
  128. C    +VERIFICATION DE IPT2
  129.       IF(.not.SELF) THEN
  130. C     un maillage non complexe
  131.       segact,IPT2
  132.       NBSOU2=IPT2.LISOUS(/1)
  133.       if (NBSOU2.gt.0) then
  134.         call erreur(25)
  135.         return
  136.       endif
  137. C     maillage de TRI3 uniquement
  138.       ITYPE2=IPT2.ITYPEL
  139.       if(ITYPE2.ne.4) then
  140.         call erreur(16)
  141.         return
  142.       endif
  143.       NBNN2  = IPT2.NUM(/1)
  144.       NBELE2 = IPT2.NUM(/2)
  145.       ENDIF
  146.  
  147.  
  148.  
  149. C-----------------------------------------------------------------------
  150. C     PREPARATION DES DONNEES DE SORTIE
  151.  
  152.       NBSOUS = 0
  153.       NBREF  = 0
  154.  
  155. C    +NOUVEAU IPT1 (dimensionné + grand)
  156.       NBNN   = NBNN1
  157.       NBELEM = NBELE1
  158.       segini,MELEME=IPT1
  159.       segdes,IPT1
  160.       IPT1=MELEME
  161.       NBELEM = 4*NBELE1
  162.       segadj,IPT1
  163.       segini,WORK1
  164.  
  165. C    +NOUVEAU IPT2
  166.       IF(.not.SELF) THEN
  167.       NBNN   = NBNN2
  168.       NBELEM = NBELE2
  169.       segini,MELEME=IPT2
  170.       segdes,IPT2
  171.       IPT2=MELEME
  172.       NBELEM = 4*NBELE2
  173.       segadj,IPT2
  174.       segini,WORK2
  175.       ENDIF
  176.  
  177. C    +IPT3 et IPT4 = lignes de seg2 = intersection des tri3
  178.       NBNN3 = 2
  179.       NBNN = NBNN3
  180.       NBELE3 = 0
  181.       NBELE4 = 0
  182.       NBELEM = 3*NBELE1
  183.       segini,IPT3
  184.       IPT3.ITYPEL=2
  185.       IF(.not.SELF) THEN
  186.         segini,IPT4
  187.         IPT4.ITYPEL=2
  188.       ENDIF
  189.  
  190.  
  191. C-----------------------------------------------------------------------
  192. C     POUR LE CAS DE LA SELF-INTERSECTION
  193.       IF(SELF) THEN
  194.         IPT2=IPT1
  195.         NBNN2=NBNN1
  196.         NBELE2=NBELE1
  197.         WORK2=WORK1
  198.         IPT4=IPT3
  199. c       attention : on va remplir IPT2 qui ne refere pas au meme element
  200. c       que IPT1 mais pas IPT4 (car on remplit deja IPT3)
  201.       ENDIF
  202.  
  203. c     pour ajuster la taille par bloc
  204.       nbloc = max(4,max(NBELE1,NBELE2))
  205.  
  206. C=======================================================================
  207. C     BOUCLES POUR LE CALCUL D INTERSECTION
  208.  
  209.  
  210. C==== BOUCLE SUR LES ELEMENTS A ========================================
  211. c     on met une grande valeur mais on s'arrete des qu on a fini :
  212. c     plus de triangles (initiaux + coupés) a traiter (=NBELE1 atteint)
  213. c       IAFIN = 1000*NBELE1
  214.       IAFIN = 100*NBELE1+100*NBELE2
  215.  
  216.       DO 1000 IA=1,IAFIN
  217.  
  218. C-----------------------------------------------------------------------
  219. C       CALCULS RELATIFS AU TRIANGLE A
  220.  
  221. c       1er coup : on prend le triangle "pere" (=avant decoupe)
  222.         IBDEB0=0
  223.  1100   CONTINUE
  224.         ICOUPA = 0
  225.         ICOUPB = 0
  226.  
  227.  
  228. c       coordonnees du Triangle A
  229.         xAmax=XPETIT
  230.         yAmax=XPETIT
  231.         zAmax=XPETIT
  232.         xAmin=XGRAND
  233.         yAmin=XGRAND
  234.         zAmin=XGRAND
  235.         DO INOD1=1,NBNN1
  236.            IPA1 = IPT1.NUM(INOD1,IA)
  237.            if(INOD1.gt.1) then
  238.            do iprec=1,INOD1-1
  239.               if(IPA1.eq.ITA(iprec)) then
  240. c                 write(ioimp,*) '2 noeuds identiques => on saute'
  241.                 goto 900
  242.               endif
  243.            enddo
  244.            endif
  245.            ITA(INOD1)=IPA1
  246.            IREF1= idim1*(IPA1-1)
  247.            XTA(INOD1,1)=XCOOR(IREF1+1)
  248.            XTA(INOD1,2)=XCOOR(IREF1+2)
  249.            XTA(INOD1,3)=XCOOR(IREF1+3)
  250.            xAmax = max(xAmax,XTA(INOD1,1))
  251.            yAmax = max(yAmax,XTA(INOD1,2))
  252.            zAmax = max(zAmax,XTA(INOD1,3))
  253.            xAmin = min(xAmin,XTA(INOD1,1))
  254.            yAmin = min(yAmin,XTA(INOD1,2))
  255.            zAmin = min(zAmin,XTA(INOD1,3))
  256.         ENDDO
  257.  
  258.         if(iimpi.ge.10) then
  259.           write(ioimp,*) '____________________________________________'
  260.           write(ioimp,*) 'TRIANGLE A # ',IA,' / ',NBELE1
  261.           write(ioimp,FMT="(' noeuds ',I5,I5,I5)")ITA(1),ITA(2),ITA(3)
  262. c           write(ioimp,*) '     coordx=',(XTA(iou,1),iou=1,3)
  263. c           write(ioimp,*) '     coordy=',(XTA(iou,2),iou=1,3)
  264. c           write(ioimp,*) '     coordz=',(XTA(iou,3),iou=1,3)
  265.         endif
  266.  
  267.         IF(WORK1.IDEJVU(IA).eq.0) THEN
  268.  
  269. c         vecteurs des aretes
  270.           va12x = XTA(2,1) - XTA(1,1)
  271.           va12y = XTA(2,2) - XTA(1,2)
  272.           va12z = XTA(2,3) - XTA(1,3)
  273.           va13x = XTA(3,1) - XTA(1,1)
  274.           va13y = XTA(3,2) - XTA(1,2)
  275.           va13z = XTA(3,3) - XTA(1,3)
  276.  
  277. c         vecteur normal au Triangle A + aire du triangle
  278.           vax = va12y*va13z - va12z*va13y
  279.           vay = va12z*va13x - va12x*va13z
  280.           vaz = va12x*va13y - va12y*va13x
  281.           xlenA = sqrt(vax*vax + vay*vay + vaz*vaz)
  282.           if (xlenA.le.XPETIT) then
  283.               write(ioimp,*) 'TRI3 A de dimension nulle! ',xlenA,XPETIT
  284. c               write(ioimp,*) 'va12=',va12x,va12y,va12z
  285. c               write(ioimp,*) 'va13=',va13x,va13y,va13z
  286. c               write(ioimp,*) 'va12 PVEC va13 = va=',vax,vay,vaz
  287.               if(IVERI.ge.1) then
  288. c           write(ioimp,*) 'NOVERIF: ON retourne le maillage en l etat'
  289. c                 goto 9000
  290. c         bp : une autre solution est de "sauter" cet element
  291.           write(ioimp,*) 'NOVERIF: l execution continue malgre tout'
  292.                 goto 900
  293.               else
  294.                 call erreur(26)
  295.                 return
  296.               endif
  297.           endif
  298.           vax = vax / xlenA
  299.           vay = vay / xlenA
  300.           vaz = vaz / xlenA
  301.           WORK1.V(1,IA) = vax
  302.           WORK1.V(2,IA) = vay
  303.           WORK1.V(3,IA) = vaz
  304.           WORK1.IDEJVU(IA)=1
  305.           if(iimpi.ge.12) write(ioimp,*) ' normale vA =',vax,vay,vaz
  306.           if(IELIM.eq.0) then
  307.             XELIM=1.D-2*xlenA
  308.             XELIM2=2.d0*XELIM
  309.             IELIM=1
  310.             XPLAN2=1.D-2*XELIM
  311.             if(XPLAN2.lt.XPLAN) then
  312.               XPLAN=XPLAN2
  313.               XPLANn=-1.d0*XPLAN
  314.             endif
  315.           endif
  316.  
  317.         ELSE
  318.           vax = WORK1.V(1,IA)
  319.           vay = WORK1.V(2,IA)
  320.           vaz = WORK1.V(3,IA)
  321.  
  322.         ENDIF
  323.  
  324.         if(IELIM.eq.0) write(ioimp,*)'Valeur par défaut : XELIM=',XELIM
  325.         if(IPLAN.eq.0) write(ioimp,*)'Valeur par défaut : XPLAN=',XPLAN
  326.         IPLAN=1
  327.         if(iimpi.ge.10) write(ioimp,*)'XELIM,XPLAN=',XELIM,XPLAN
  328.         if(iimpi.ge.12) write(ioimp,*)'XELIM2,XPLANn=',XELIM2,XPLANn
  329.  
  330. c       on definit une boite limitant le domaine de recherche
  331.         xAmax = xAmax + (XELIM2)
  332.         yAmax = yAmax + (XELIM2)
  333.         zAmax = zAmax + (XELIM2)
  334.         xAmin = xAmin - (XELIM2)
  335.         yAmin = yAmin - (XELIM2)
  336.         zAmin = zAmin - (XELIM2)
  337.  
  338.  
  339. C====== BOUCLE SUR LES ELEMENTS B ======================================
  340.  
  341. c       dans le cas d'une decoupe, on reprend avec l element IB suivant
  342.         IF(IBDEB0.ne.0) THEN
  343.           IBDEB=IBDEB0
  344.         ELSE
  345. c         dans le cas de la self-intersection on teste les triangles d'apres
  346.           IF(SELF) THEN
  347.             IBDEB=IA+1
  348.           ELSE
  349.             IBDEB=1
  350.           ENDIF
  351.         ENDIF
  352.         IBFIN=NBELE2
  353.  
  354.         DO 2000 IB=IBDEB,IBFIN
  355.  
  356. c         if(iimpi.ge.10) write(ioimp,*) '+ TRIANGLE B # ',IB,' / ',NBELE2
  357.  
  358. C-----------------------------------------------------------------------
  359. C         CALCUL DE LA DISTANCE AU PLAN TA DES NOEUDS DU Triangle B
  360.  
  361. c         coordonnees du triangle B
  362.           ibxmax=0
  363.           ibymax=0
  364.           ibzmax=0
  365.           ibxmin=0
  366.           ibymin=0
  367.           ibzmin=0
  368.           DO INOD2=1,NBNN2
  369.              IPB2 = IPT2.NUM(INOD2,IB)
  370.              if(INOD2.gt.1) then
  371.              do iprec=1,INOD2-1
  372.                 if(IPB2.eq.ITB(iprec)) then
  373.                   write(ioimp,*) '2 noeuds identiques de TB => on saute'
  374.                   goto 1900
  375.                 endif
  376.              enddo
  377.              endif
  378.              ITB(INOD2)=IPB2
  379.              IREF2= idim1*(IPB2-1)
  380.              XTB(INOD2,1)=XCOOR(IREF2+1)
  381.              XTB(INOD2,2)=XCOOR(IREF2+2)
  382.              XTB(INOD2,3)=XCOOR(IREF2+3)
  383.              if(XTB(INOD2,1).gt.xAmax) ibxmax=ibxmax+1
  384.              if(XTB(INOD2,2).gt.yAmax) ibymax=ibymax+1
  385.              if(XTB(INOD2,3).gt.zAmax) ibzmax=ibzmax+1
  386.              if(XTB(INOD2,1).lt.xAmin) ibxmin=ibxmin+1
  387.              if(XTB(INOD2,2).lt.yAmin) ibymin=ibymin+1
  388.              if(XTB(INOD2,3).lt.zAmin) ibzmin=ibzmin+1
  389.           ENDDO
  390.  
  391. c         si on est hors boite selon x y ou z on saute
  392.           if(ibxmax.eq.NBNN2.or.ibymax.eq.NBNN2.or.ibzmax.eq.NBNN2)
  393.      &    goto 1900
  394.           if(ibxmin.eq.NBNN2.or.ibymin.eq.NBNN2.or.ibzmin.eq.NBNN2)
  395.      &    goto 1900
  396.  
  397. c         distance au plan TA
  398.           do iidim=1,idim
  399.             dB1(iidim) = XTB(1,iidim) - XTA(1,iidim)
  400.             dB2(iidim) = XTB(2,iidim) - XTA(1,iidim)
  401.             dB3(iidim) = XTB(3,iidim) - XTA(1,iidim)
  402.           enddo
  403.           xddB1 = vax*dB1(1) + vay*dB1(2) + vaz*dB1(3)
  404.           xddB2 = vax*dB2(1) + vay*dB2(2) + vaz*dB2(3)
  405.           xddB3 = vax*dB3(1) + vay*dB3(2) + vaz*dB3(3)
  406.  
  407. c         si tous les points de TB sont en dessous (ou dessus) de TA
  408. c         + le cas affleurant
  409. c         => il n'y a pas d'intersection : on saute
  410.           if(xddB1.lt.XPLANn.and.xddB2.lt.XPLANn.and.xddB3.lt.XPLANn)
  411.      &    goto 1900
  412.           if(xddB1.gt.XPLAN.and.xddB2.gt.XPLAN.and.xddB3.gt.XPLAN)
  413.      &    goto 1900
  414.  
  415. c         combien de points de TB dans le plan de TA ? => iBinA
  416.           iBinA = 0
  417.           if(abs(xddB1).le.XPLAN) iBinA=iBinA+1
  418.           if(abs(xddB2).le.XPLAN) iBinA=iBinA+2
  419.           if(abs(xddB3).le.XPLAN) iBinA=iBinA+4
  420.  
  421. c         si tous les points de TB sont dans le plan de TA,
  422. c         => on doit traiter le cas coplanaire
  423. c Il n'est pas traité pour l'instant car assez rare en pratique en 3D.
  424.           if (iBinA.eq.7) then
  425. c L'intersection générée pourrait etre consitutée de polygones [Gander]
  426. c etape 1: calcul des intersection segment de A - segments de B
  427. c       2: recherche des points de A (/B) a l'interieur de TB (/TA)
  428. c       3: ordonner dans le sens trigo les points des étapes 1 et 2
  429. c       4: eliminer les doublons
  430.             goto 1900
  431.           endif
  432.  
  433. c         si 1 seul point de TB est dans le plan de TA,
  434. c            et les autres dont du meme coté,
  435. c         => on ne traite pas cette intersection ponctuelle : on saute
  436.           proB12 = xddB1*xddB2
  437.           proB23 = xddB2*xddB3
  438.           proB31 = xddB3*xddB1
  439.           if(iBinA.eq.1.and.proB23.gt.0.d0) goto 1900
  440.           if(iBinA.eq.2.and.proB31.gt.0.d0) goto 1900
  441.           if(iBinA.eq.4.and.proB12.gt.0.d0) goto 1900
  442.  
  443. c         si 2 points de TB sont dans le plan de TA
  444. c         ou si on a une traversée franche => on traite
  445.           if(iimpi.ge.10) then
  446.             write(ioimp,*) '+ TRIANGLE B # ',IB,' / ',NBELE2
  447.             write(ioimp,FMT="('noeuds ',I5,I5,I5)")ITB(1),ITB(2),ITB(3)
  448. c             write(ioimp,*) '+     coordx ',(XTB(iou,1),iou=1,3)
  449. c             write(ioimp,*) '+     coordy ',(XTB(iou,2),iou=1,3)
  450. c             write(ioimp,*) '+     coordz ',(XTB(iou,3),iou=1,3)
  451.             if(iimpi.ge.12)
  452.      &      write(ioimp,*) ' xddB1,2,3=',xddB1,xddB2,xddB3
  453.           endif
  454.  
  455. C-----------------------------------------------------------------------
  456. C         CALCULS RELATIFS AU TRIANGLE B
  457.  
  458.           IF(WORK2.IDEJVU(IB).eq.0) THEN
  459.  
  460. c           vecteurs des aretes
  461.             vb12x = XTB(2,1) - XTB(1,1)
  462.             vb12y = XTB(2,2) - XTB(1,2)
  463.             vb12z = XTB(2,3) - XTB(1,3)
  464.             vb13x = XTB(3,1) - XTB(1,1)
  465.             vb13y = XTB(3,2) - XTB(1,2)
  466.             vb13z = XTB(3,3) - XTB(1,3)
  467.  
  468. c           vecteur normal au Triangle B + aire du triangle
  469.             vbx = vb12y*vb13z - vb12z*vb13y
  470.             vby = vb12z*vb13x - vb12x*vb13z
  471.             vbz = vb12x*vb13y - vb12y*vb13x
  472.             xlenB = sqrt(vbx*vbx + vby*vby + vbz*vbz)
  473.             if (xlenB.le.XPETIT) then
  474.               write(ioimp,*) 'TRI3 B de dimension nulle! ',xlenB,XPETIT
  475.               if(IVERI.ge.1) then
  476. c             write(ioimp,*) 'NOVERIF: ON retourne le maillage en l etat'
  477. c                 goto 9000
  478. c         bp : une autre solution est de "sauter" cet element
  479.           write(ioimp,*) 'NOVERIF: l execution continue malgre tout'
  480.                 goto 1900
  481.               else
  482.                 call erreur(26)
  483.                 return
  484.               endif
  485.             endif
  486.             vbx = vbx / xlenB
  487.             vby = vby / xlenB
  488.             vbz = vbz / xlenB
  489.             WORK2.V(1,IB) = vbx
  490.             WORK2.V(2,IB) = vby
  491.             WORK2.V(3,IB) = vbz
  492.             WORK2.IDEJVU(IB) = 1
  493.             if(iimpi.ge.12) write(ioimp,*) ' normale vB =',vbx,vby,vbz
  494.  
  495.           ELSE
  496.             vbx = WORK2.V(1,IB)
  497.             vby = WORK2.V(2,IB)
  498.             vbz = WORK2.V(3,IB)
  499.  
  500.           ENDIF
  501.  
  502. C-----------------------------------------------------------------------
  503. C         CALCUL DE LA DISTANCE AU PLAN TB DES NOEUDS DU Triangle A
  504.  
  505.           do iidim=1,idim
  506.             dA1(iidim) = XTA(1,iidim) - XTB(1,iidim)
  507.             dA2(iidim) = XTA(2,iidim) - XTB(1,iidim)
  508.             dA3(iidim) = XTA(3,iidim) - XTB(1,iidim)
  509.           enddo
  510.           xddA1 = vbx*dA1(1) + vby*dA1(2) + vbz*dA1(3)
  511.           xddA2 = vbx*dA2(1) + vby*dA2(2) + vbz*dA2(3)
  512.           xddA3 = vbx*dA3(1) + vby*dA3(2) + vbz*dA3(3)
  513.         if(iimpi.ge.12) write(ioimp,*) ' xddA1,2,3=',xddA1,xddA2,xddA3
  514.  
  515. c         si tous les points de TA sont en dessous (ou dessus) de TB,
  516. c         il n'y a pas d'intersection : on quitte
  517. c           if(xddA1.lt.xpetA.and.xddA2.lt.xpetA.and.xddA3.lt.xpetA)
  518. c      &    goto 1900
  519. c           if(xddA1.gt.xpetAn.and.xddA2.gt.xpetAn.and.xddA3.gt.xpetAn)
  520. c      &    goto 1900
  521. c         idem, mais on traite desormais le cas affleurant
  522.           if(xddA1.lt.XPLANn.and.xddA2.lt.XPLANn.and.xddA3.lt.XPLANn)
  523.      &    goto 1900
  524.           if(xddA1.gt.XPLAN.and.xddA2.gt.XPLAN.and.xddA3.gt.XPLAN)
  525.      &    goto 1900
  526.  
  527. c         combien de points de TA dans le plan de TB ? => iAinB
  528.           iAinB = 0
  529.           if(abs(xddA1).le.XPLAN) iAinB=iAinB+1
  530.           if(abs(xddA2).le.XPLAN) iAinB=iAinB+2
  531.           if(abs(xddA3).le.XPLAN) iAinB=iAinB+4
  532.  
  533. c         si tous les points de TA sont dans le plan de TB,
  534. c-------- => on doit traiter le cas coplanaire ------------------------
  535.           if (iAinB.eq.7.or.iBinA.eq.7) then
  536.             write(ioimp,*) 'CAS COPLANAIRE EN COURS DE DEVELOPPEMENT'
  537.             goto 1900
  538.  
  539.           elseif(SELF) then
  540. c           si self et 2 points en commun et pas coplanaire => on saute
  541.             ichev = 0
  542.             do INOD2=1,NBNN2
  543.               IPB2 = IPT2.NUM(INOD2,IB)
  544.               if(IPB2.eq.ITA(1)) ichev=ichev+1
  545.               if(IPB2.eq.ITA(2)) ichev=ichev+1
  546.               if(IPB2.eq.ITA(3)) ichev=ichev+1
  547.             enddo
  548.             if(ichev.eq.2) goto 1900
  549.           endif
  550. c-------- fin du cas coplanaire---------------------------------------
  551.  
  552.  
  553. c         si 1 seul point de TA est dans le plan de TB
  554. c            et les autres dont du meme coté,
  555. c         => on ne traite pas cette intersection ponctuelle : on saute
  556.           proA12 = xddA1*xddA2
  557.           proA23 = xddA2*xddA3
  558.           proA31 = xddA3*xddA1
  559.           if(iAinB.eq.1.and.proA23.gt.0.d0) goto 1900
  560.           if(iAinB.eq.2.and.proA31.gt.0.d0) goto 1900
  561.           if(iAinB.eq.4.and.proA12.gt.0.d0) goto 1900
  562.  
  563. c         si 2 points de TA sont dans le plan de TB
  564. c         ou si on a une traversée franche => on traite
  565.  
  566. C-----------------------------------------------------------------------
  567. C         VECTEUR DIRECTEUR DE LA DROITE INTERSECTANT LES 2 PLANS
  568.           dx = vay*vbz - vaz*vby
  569.           dy = vaz*vbx - vax*vbz
  570.           dz = vax*vby - vay*vbx
  571.         if(iimpi.ge.12) write(ioimp,*) ' directrice D=',dx,dy,dz
  572.  
  573.  
  574. C-----------------------------------------------------------------------
  575. C         EVENTUELLE PERMUTATION DES NOEUDS DE TA
  576. c         afin que xddA1 soit de signe opposé à celui de xddA2 et xddA3
  577. c         on fait attention aux cas ou on a des noeuds de A dans plan TB
  578. c         noeuds A1 et A2 dans plan TB
  579.           if(iAinB.eq.3) then
  580.             proA12 =  1.d0
  581.             proA23 = -1.d0
  582. c         noeuds A1 et A3 dans plan TB
  583.           elseif(iAinB.eq.5) then
  584.             proA12 = -1.d0
  585.             proA23 = -1.d0
  586. c         noeuds A2 et A3 dans plan TB
  587.           elseif(iAinB.eq.6) then
  588.             proA12 = 0.d0
  589.             proA23 = 1.d0
  590.           else
  591.             proA12 = xddA1*xddA2
  592.             proA23 = xddA2*xddA3
  593.           endif
  594.           IF(proA23.ge.0.d0) THEN
  595. c         si xddA2 et xddA3 de meme signe, rien a faire
  596. c rem: ge possible (au lieu de gt) car on est sur de traverser franchement
  597.             iA1=1
  598.             iA2=2
  599.             iA3=3
  600.           ELSE
  601. c           xddA2 et xddA3 de signe oppose : il faut permuter
  602.             IF(proA12.le.0.d0) THEN
  603. c             xddA2 different des 2 autres : 1<-2 2<-3 3<-1
  604.               iA1=2
  605.               iA2=3
  606.               iA3=1
  607.               xddtmp= xddA1
  608.               xddA1 = xddA2
  609.               xddA2 = xddA3
  610.               xddA3 = xddtmp
  611.             ELSE
  612. c             xddA3 different des 2 autres : 1<-3 2<-1 3<-2
  613.               iA1=3
  614.               iA2=1
  615.               iA3=2
  616.               xddtmp= xddA2
  617.               xddA2 = xddA1
  618.               xddA1 = xddA3
  619.               xddA3 = xddtmp
  620.             ENDIF
  621.           ENDIF
  622.           va12x = XTA(iA2,1) - XTA(iA1,1)
  623.           va12y = XTA(iA2,2) - XTA(iA1,2)
  624.           va12z = XTA(iA2,3) - XTA(iA1,3)
  625.           va13x = XTA(iA3,1) - XTA(iA1,1)
  626.           va13y = XTA(iA3,2) - XTA(iA1,2)
  627.           va13z = XTA(iA3,3) - XTA(iA1,3)
  628. c         on recalcule  iAinB
  629.           iAinB = 0
  630.           if(abs(xddA1).le.XPLAN) iAinB=iAinB+1
  631.           if(abs(xddA2).le.XPLAN) iAinB=iAinB+2
  632.           if(abs(xddA3).le.XPLAN) iAinB=iAinB+4
  633. c ici A1 (=1er noeud du Triangle A) est situé de l'autre coté du plan TB
  634.         if(iimpi.ge.12) then
  635.         write(ioimp,*) ' noeuds A permutés=',ITA(iA1),ITA(iA2),ITA(iA3)
  636.         write(ioimp,*) ' xddA1,2,3=',xddA1,xddA2,xddA3
  637.         endif
  638.  
  639. C-----------------------------------------------------------------------
  640. C         EVENTUELLE PERMUTATION DES NOEUDS DE TB
  641. c         afin que xddB1 soit de signe opposé à celui de xddB2 et xddB3
  642. c         on fait attention aux cas ou on a des noeuds de A dans plan TB
  643.           if(iBinA.eq.3) then
  644.             proB12 =  1.d0
  645.             proB23 = -1.d0
  646.           elseif(iBinA.eq.5) then
  647.             proB12 = -1.d0
  648.             proB23 = -1.d0
  649.           elseif(iBinA.eq.6) then
  650.             proB12 = 0.d0
  651.             proB23 = 1.d0
  652.           else
  653.             proB12 = xddB1*xddB2
  654.             proB23 = xddB2*xddB3
  655.           endif
  656.           IF(proB23.ge.0.d0) THEN
  657. c         si xddB2 et xddB3 de meme signe, rien a faire
  658. c rem: ge possible (au lieu de gt) car on est sur de traverser franchement
  659.             iB1=1
  660.             iB2=2
  661.             iB3=3
  662.           ELSE
  663. c           xddB2 et xddB3 de signe oppose : il faut permuter
  664.             IF(proB12.le.0.d0) THEN
  665. c             xddB2 different des 2 autres : 1<-2 2<-3 3<-1
  666.               iB1=2
  667.               iB2=3
  668.               iB3=1
  669.               xddtmp= xddB1
  670.               xddB1 = xddB2
  671.               xddB2 = xddB3
  672.               xddB3 = xddtmp
  673.             ELSE
  674. c             xddB3 different des 2 autres : 1<-3 2<-1 3<-2
  675.               iB1=3
  676.               iB2=1
  677.               iB3=2
  678.               xddtmp= xddB2
  679.               xddB2 = xddB1
  680.               xddB1 = xddB3
  681.               xddB3 = xddtmp
  682.             ENDIF
  683.           ENDIF
  684.           vb12x = XTB(iB2,1) - XTB(iB1,1)
  685.           vb12y = XTB(iB2,2) - XTB(iB1,2)
  686.           vb12z = XTB(iB2,3) - XTB(iB1,3)
  687.           vb13x = XTB(iB3,1) - XTB(iB1,1)
  688.           vb13y = XTB(iB3,2) - XTB(iB1,2)
  689.           vb13z = XTB(iB3,3) - XTB(iB1,3)
  690. c         on recalcule iBinA
  691.           iBinA = 0
  692.           if(abs(xddB1).le.XPLAN) iBinA=iBinA+1
  693.           if(abs(xddB2).le.XPLAN) iBinA=iBinA+2
  694.           if(abs(xddB3).le.XPLAN) iBinA=iBinA+4
  695. c ici B1 (=1er noeud du Triangle B) est situé de l'autre coté du plan TB
  696.         if(iimpi.ge.12) then
  697.         write(ioimp,*) ' noeuds B permutés=',ITB(iB1),ITB(iB2),ITB(iB3)
  698.         write(ioimp,*) ' xddB1,2,3=',xddB1,xddB2,xddB3
  699.         endif
  700.  
  701.  
  702. C-----------------------------------------------------------------------
  703. C         CALCUL DES POINTS P2 et P3
  704.  
  705. c         P2 = intersection arete 1-2 de TA avec plan TB
  706.           rap2 = abs(xddA1/(xddA1-xddA2))
  707.           xP2(1) = XTA(iA1,1) + rap2*va12x
  708.           xP2(2) = XTA(iA1,2) + rap2*va12y
  709.           xP2(3) = XTA(iA1,3) + rap2*va12z
  710. c         P3 = intersection arete 1-3 de TA avec plan TB
  711.           rap3 = abs(xddA1/(xddA1-xddA3))
  712.           xP3(1) = XTA(iA1,1) + rap3*va13x
  713.           xP3(2) = XTA(iA1,2) + rap3*va13y
  714.           xP3(3) = XTA(iA1,3) + rap3*va13z
  715.  
  716. c         pour y repondre il faut calculer les points Q2 et Q3
  717.           if(iimpi.ge.12) write(ioimp,*) ' P2=',(xP2(iou),iou=1,3)
  718.           if(iimpi.ge.12) write(ioimp,*) ' P3=',(xP3(iou),iou=1,3)
  719.  
  720. c         le segment [P1-P2] est il inclut dans TB ?
  721. C-----------------------------------------------------------------------
  722. C         CALCUL DES POINTS Q2 et Q3
  723.  
  724. c         Q2 = intersection arete 1-2 de TB avec plan TA
  725.           rap2 = abs(xddB1/(xddB1-xddB2))
  726.           xQ2(1) = XTB(iB1,1) + rap2*vb12x
  727.           xQ2(2) = XTB(iB1,2) + rap2*vb12y
  728.           xQ2(3) = XTB(iB1,3) + rap2*vb12z
  729. c         Q3 = intersection arete 1-3 de TB avec plan TB
  730.           rap3 = abs(xddB1/(xddB1-xddB3))
  731.           xQ3(1) = XTB(iB1,1) + rap3*vb13x
  732.           xQ3(2) = XTB(iB1,2) + rap3*vb13y
  733.           xQ3(3) = XTB(iB1,3) + rap3*vb13z
  734.  
  735.           if(iimpi.ge.12) write(ioimp,*) ' Q2=',(xQ2(iou),iou=1,3)
  736.           if(iimpi.ge.12) write(ioimp,*) ' Q3=',(xQ3(iou),iou=1,3)
  737.  
  738.  
  739. C-----------------------------------------------------------------------
  740. C         IDENTIFICATION EVENTUELLE DES POINTS P2 P3 Q2 ET Q3
  741. c         calcul de la distance de P2 aux noeuds A1 et A2
  742.           IP2=0
  743.           dP2A2 = sqrt((xP2(1)-XTA(IA2,1))**2
  744.      &               + (xP2(2)-XTA(IA2,2))**2
  745.      &               + (xP2(3)-XTA(IA2,3))**2)
  746.           if (dP2A2.le.XELIM2) then
  747.              IP2   =ITA(IA2)
  748.              xP2(1)=XTA(IA2,1)
  749.              xP2(2)=XTA(IA2,2)
  750.              xP2(3)=XTA(IA2,3)
  751.           else
  752.              dP2A1 = sqrt((xP2(1)-XTA(IA1,1))**2
  753.      &                  + (xP2(2)-XTA(IA1,2))**2
  754.      &                  + (xP2(3)-XTA(IA1,3))**2)
  755.              if (dP2A1.le.XELIM2) then
  756.                 IP2   =ITA(IA1)
  757.                 xP2(1)=XTA(IA1,1)
  758.                 xP2(2)=XTA(IA1,2)
  759.                 xP2(3)=XTA(IA1,3)
  760.              endif
  761.           endif
  762. c         calcul de la distance de P3 aux noeuds A1 et A3
  763.           IP3=0
  764.           dP3A3 = sqrt((xP3(1)-XTA(IA3,1))**2
  765.      &               + (xP3(2)-XTA(IA3,2))**2
  766.      &               + (xP3(3)-XTA(IA3,3))**2)
  767.           if (dP3A3.le.XELIM2) then
  768.              IP3   =ITA(IA3)
  769.              xP3(1)=XTA(IA3,1)
  770.              xP3(2)=XTA(IA3,2)
  771.              xP3(3)=XTA(IA3,3)
  772.           else
  773.              dP3A1 = sqrt((xP3(1)-XTA(IA1,1))**2
  774.      &                  + (xP3(2)-XTA(IA1,2))**2
  775.      &                  + (xP3(3)-XTA(IA1,3))**2)
  776.              if (dP3A1.le.XELIM2) then
  777.                 IP3   =ITA(IA1)
  778.                 xP3(1)=XTA(IA1,1)
  779.                 xP3(2)=XTA(IA1,2)
  780.                 xP3(3)=XTA(IA1,3)
  781.              endif
  782.           endif
  783.  
  784. c         calcul de la distance de Q2 aux noeuds B1 et B2
  785.           IQ2=0
  786.           dQ2B2 = sqrt((xQ2(1)-XTB(IB2,1))**2
  787.      &               + (xQ2(2)-XTB(IB2,2))**2
  788.      &               + (xQ2(3)-XTB(IB2,3))**2)
  789.           if (dQ2B2.le.XELIM2) then
  790.              IQ2  = ITB(IB2)
  791.              xQ2(1)=XTB(IB2,1)
  792.              xQ2(2)=XTB(IB2,2)
  793.              xQ2(3)=XTB(IB2,3)
  794.           else
  795.              dQ2B1 = sqrt((xQ2(1)-XTB(IB1,1))**2
  796.      &                  + (xQ2(2)-XTB(IB1,2))**2
  797.      &                  + (xQ2(3)-XTB(IB1,3))**2)
  798.              if (dQ2B1.le.XELIM2) then
  799.                 IQ2  = ITB(IB1)
  800.                 xQ2(1)=XTB(IB1,1)
  801.                 xQ2(2)=XTB(IB1,2)
  802.                 xQ2(3)=XTB(IB1,3)
  803.              endif
  804.           endif
  805. c         calcul de la distance de Q3 aux noeuds B1 et B3
  806.           IQ3=0
  807.           dQ3B3 = sqrt((xQ3(1)-XTB(IB3,1))**2
  808.      &               + (xQ3(2)-XTB(IB3,2))**2
  809.      &               + (xQ3(3)-XTB(IB3,3))**2)
  810.           if (dQ3B3.le.XELIM2) then
  811.              IQ3   =ITB(IB3)
  812.              xQ3(1)=XTB(IB3,1)
  813.              xQ3(2)=XTB(IB3,2)
  814.              xQ3(3)=XTB(IB3,3)
  815.           else
  816.              dQ3B1 = sqrt((xQ3(1)-XTB(IB1,1))**2
  817.      &                  + (xQ3(2)-XTB(IB1,2))**2
  818.      &                  + (xQ3(3)-XTB(IB1,3))**2)
  819.              if (dQ3B1.le.XELIM2) then
  820.                 IQ3   =ITB(IB1)
  821.                 xQ3(1)=XTB(IB1,1)
  822.                 xQ3(2)=XTB(IB1,2)
  823.                 xQ3(3)=XTB(IB1,3)
  824.              endif
  825.           endif
  826. c rem bp : cette identification déplace la position de P2 P3 Q2 et Q3
  827. c          vers les A1 A2 etc...
  828.  
  829.  
  830. C-----------------------------------------------------------------------
  831. C         ABCISSE s SUR D DES POINTS P2, P3, Q2 et Q3
  832.           sP2 = xP2(1)*dx + xP2(2)*dy + xP2(3)*dz
  833.           sP3 = xP3(1)*dx + xP3(2)*dy + xP3(3)*dz
  834.           sQ2 = xQ2(1)*dx + xQ2(2)*dy + xQ2(3)*dz
  835.           sQ3 = xQ3(1)*dx + xQ3(2)*dy + xQ3(3)*dz
  836.           if(iimpi.ge.10) then
  837.             WRITE(IOIMP,*) 'identif   #P2,3=',IP2,IP3,' Q2,3=',IQ2,IQ3
  838.             WRITE(IOIMP,*) 'abscisses  P2,3=',sP2,sP3,' Q2,3=',sQ2,sQ3
  839.           endif
  840.  
  841. c         si les segments [P2-P3] et [Q2-Q3] ne se chevauchent pas, on quitte
  842.           if(sP2.lt.sP3) then
  843.            sminP = sP2
  844.            smaxP = sP3
  845.           else
  846.            sminP = sP3
  847.            smaxP = sP2
  848.           endif
  849.           if(sQ2.lt.sQ3) then
  850.            sminQ = sQ2
  851.            smaxQ = sQ3
  852.           else
  853.            sminQ = sQ3
  854.            smaxQ = sQ2
  855.           endif
  856.  
  857. c         si pas de chevauchement (exclusion), on saute
  858.           if (sminP.ge.smaxQ.or.smaxP.le.sminQ) goto 1900
  859.  
  860.  
  861. C-----------------------------------------------------------------------
  862. C         PREPARATION DES MAILLAGES DE TA et TB + SEGMENT D'INTERSECTION
  863. c            write(6,*) 'place dispo dan ipt1=',IPT1.NUM(/2),NBELE1
  864. c            write(6,*) 'place dispo dan ipt2=',IPT2.NUM(/2),NBELE2
  865. c         Avant tout on vérifie qu'il reste assez de place
  866.           if(NBELE1+4.gt.IPT1.NUM(/2)) then
  867.             NBNN = NBNN1
  868.             NBELEM = NBELE1+nbloc
  869. c             segadj,IPT1
  870.             segadj,IPT1,WORK1
  871.           endif
  872.           if(NBELE2+4.gt.IPT2.NUM(/2)) then
  873.             NBNN = NBNN2
  874.             NBELEM = NBELE2+nbloc
  875.             segadj,IPT2,WORK2
  876.           endif
  877.           if(NBELE3+1.gt.IPT3.NUM(/2)) then
  878.             NBNN = NBNN3
  879.             NBELEM = NBELE3+nbloc
  880.             segadj,IPT3
  881.           endif
  882.           if(NBELE4+1.gt.IPT4.NUM(/2)) then
  883.             NBNN = NBNN3
  884.             NBELEM = NBELE4+nbloc
  885.             segadj,IPT4
  886.           endif
  887.  
  888.  
  889. c         si chevauchement exact et 2 points de B dans plan A
  890. c         et 2 points de A dans plan B => on saute
  891.           if (   (abs(sminP-sminQ).le.XELIM)
  892.      &      .and.(abs(smaxP-smaxQ).le.XELIM) ) then
  893. c           on ajoute les segments a IPT3 et IPT4
  894.             if (iAinB.eq.6.and.iBinA.eq.6) then
  895. c            +on crée les segments intersection
  896.               NBTMP=NBELE3+1
  897.               NODE(1)=ITA(IA2)
  898.               NODE(2)=ITA(IA3)
  899.               call ajouel(NODE,2,2,IPT3,NBTMP)
  900.               NBELE3=max(NBELE3,NBTMP)
  901.               if(.not.SELF) then
  902.                 NBTMP=NBELE4+1
  903.                 NODE(1)=ITB(IB2)
  904.                 NODE(2)=ITB(IB3)
  905.                 call ajouel(NODE,2,2,IPT4,NBTMP)
  906.                 NBELE4=max(NBELE4,NBTMP)
  907.               endif
  908.               if(iimpi.ge.10) then
  909.                WRITE(IOIMP,*) 'chevauchement exact IPT3 elem',NBTMP
  910.                WRITE(IOIMP,*) IPT3.NUM(1,NBTMP),IPT3.NUM(2,NBTMP)
  911.               endif
  912.               goto 1900
  913.             endif
  914.           endif
  915.  
  916. c         arrivé ici, on a detecté une intersection entre TA et TB
  917.           if(iimpi.ge.10) then
  918.            WRITE(IOIMP,*) 'INTERSECTION DETECTÉE entre ',IA,' et ',IB
  919.            WRITE(IOIMP,*) 'P=[',sminP,smaxP,']  Q=[',sminQ,smaxQ,']'
  920.           endif
  921.  
  922.  
  923.  
  924.  
  925. C-----------------------------------------------------------------------
  926. C         PREPARATION DE LA CONSTRUCTION DES FILS DE TA
  927. C         + CONSTRUCTION DU SEGMENT D'INTERSECTION
  928.  
  929.  
  930. c     --- Cas [P2-P3] inclus dans [Q2-Q3] => segment P2-P3 ---
  931.           IF (sminQ.lt.(sminP+XELIM).and.(smaxP-XELIM).lt.smaxQ) THEN
  932.  
  933. c          -on est sur le bord du triangle => on ne decoupe pas A
  934.             if(IP2.ne.0.and.IP3.ne.0) then
  935.               ICOUPA=0
  936. c          -on a 2 points confondus => TA decoupe en 2
  937.             elseif(IP2.ne.0) then
  938.               call AJOUPO(xP3,IPT3,XELIM,IP3)
  939.               ICOUPA=23
  940.             elseif(IP3.ne.0) then
  941.               call AJOUPO(xP2,IPT3,XELIM,IP2)
  942.               ICOUPA=22
  943. c          -on est dans le cas general => TA decoupe en 3
  944.             else
  945. c             on ajoute P2 et P3 a XCOOR
  946.               call AJOUPO(xP2,IPT3,XELIM,IP2)
  947.               call AJOUPO(xP3,IPT3,XELIM,IP3)
  948.               ICOUPA=30
  949.             endif
  950. c          +on crée le segment intersection
  951.             NBTMP=NBELE3+1
  952.             NODE(1)=IP2
  953.             NODE(2)=IP3
  954.             call ajouel(NODE,2,2,IPT3,NBTMP)
  955.             NBELE3=max(NBELE3,NBTMP)
  956.  
  957. c     --- Cas [P2-P3] exclus de [Q2-Q3] => segment Q2-Q3 ---
  958.           ELSEIF(sminQ.gt.(sminP+XELIM).and.(smaxP-XELIM).gt.smaxQ)THEN
  959.  
  960. c          -on est sur le bord du triangle => TA decoupe en 3
  961.             if(IP2.ne.0.and.IP3.ne.0) then
  962.               ICOUPA=31
  963. c          -on a 2 points confondus => TA decoupe en 5
  964.             else
  965.               if(abs(xddA2).ge.abs(xddA3)) then
  966.                 ICOUPA=52
  967.               else
  968.                 ICOUPA=53
  969.               endif
  970.             endif
  971.               call AJOUPO(xQ2,IPT3,XELIM,IQQ2)
  972. c             cas particulier IQ2 = IQ3
  973.               if(IQ2.eq.IQ3) then
  974.                 IQQ3=IQQ2
  975.               else
  976.                 call AJOUPO(xQ3,IPT3,XELIM,IQQ3)
  977.               endif
  978. c           petite astuce pour garder l orientation
  979.             if ( (sP2.lt.sP3.and.sQ2.lt.sQ3)
  980.      &       .or.(sP2.gt.sP3.and.sQ2.gt.sQ3) ) then
  981.               IQQ2=IQQ2
  982.               IQQ3=IQQ3
  983.             else
  984.               IQQtmp=IQQ3
  985.               IQQ3=IQQ2
  986.               IQQ2=IQQtmp
  987.             endif
  988. c          +on crée le segment intersection
  989.               NBTMP=NBELE3+1
  990.               NODE(1)=IQQ2
  991.               NODE(2)=IQQ3
  992.               call ajouel(NODE,2,2,IPT3,NBTMP)
  993.               NBELE3=max(NBELE3,NBTMP)
  994. c             endif
  995.  
  996. c     --- Cas chevauchement selon [Q2 P3] ---
  997.           ELSEIF (((sP2+XELIM).le.sQ2.and.sQ2.le.(sP3-XELIM)
  998.      &              .and.sQ2.le.sQ3)
  999.      &        .OR.((sP2-XELIM).ge.sQ2.and.sQ2.ge.(sP3+XELIM)
  1000.      &              .and.sQ2.ge.sQ3)) THEN
  1001.  
  1002. c          -on a a minima 2 points confondus
  1003.             if(IP3.ne.0) then
  1004. c             -on est sur le bord du triangle => TA decoupe en 2
  1005.               if(IP2.ne.0) then
  1006.                 ICOUPA=21
  1007. c             -on a 2 points confondus => TA decoupe en 3
  1008.               else
  1009.                 ICOUPA=39
  1010.               endif
  1011.               call AJOUPO(xQ2,IPT3,XELIM,IQQ)
  1012. c          -TA decoupe en 4
  1013.             else
  1014.               call AJOUPO(xP3,IPT3,XELIM,IP3)
  1015.               call AJOUPO(xQ2,IPT3,XELIM,IQQ)
  1016.               ICOUPA=43
  1017.             endif
  1018. c          +on crée le segment intersection
  1019.               NBTMP=NBELE3+1
  1020.               NODE(1)=IQQ
  1021.               NODE(2)=IP3
  1022.               call ajouel(NODE,2,2,IPT3,NBTMP)
  1023.               NBELE3=max(NBELE3,NBTMP)
  1024. c             endif
  1025.  
  1026. c     --- Cas chevauchement selon [Q3 P3] ---
  1027.           ELSEIF (((sP2+XELIM).le.sQ3.and.sQ3.le.(sP3-XELIM)
  1028.      &              .and.sQ3.le.sQ2)
  1029.      &        .OR.((sP2-XELIM).ge.sQ3.and.sQ3.ge.(sP3+XELIM)
  1030.      &              .and.sQ3.ge.sQ2)) THEN
  1031.  
  1032. c          -on a a minima 2 points confondus
  1033.             if(IP3.ne.0) then
  1034. c             -on est sur le bord du triangle => TA decoupe en 2
  1035.               if(IP2.ne.0) then
  1036.                 ICOUPA=21
  1037. c             -on a 2 points confondus => TA decoupe en 3
  1038.               else
  1039.                 ICOUPA=39
  1040.               endif
  1041.               call AJOUPO(xQ3,IPT3,XELIM,IQQ)
  1042. c          -TA decoupe en 4
  1043.             else
  1044.               call AJOUPO(xP3,IPT3,XELIM,IP3)
  1045.               call AJOUPO(xQ3,IPT3,XELIM,IQQ)
  1046.               ICOUPA=43
  1047.             endif
  1048. c          +on crée le segment intersection
  1049.               NBTMP=NBELE3+1
  1050.               NODE(1)=IQQ
  1051.               NODE(2)=IP3
  1052.               call ajouel(NODE,2,2,IPT3,NBTMP)
  1053.               NBELE3=max(NBELE3,NBTMP)
  1054.  
  1055. c     --- Cas chevauchement selon [Q2 P2] ---
  1056.           ELSEIF (((sP3+XELIM).le.sQ2.and.sQ2.le.(sP2-XELIM)
  1057.      &              .and.sQ2.le.sQ3)
  1058.      &        .OR.((sP3-XELIM).ge.sQ2.and.sQ2.ge.(sP2+XELIM)
  1059.      &              .and.sQ2.ge.sQ3)) THEN
  1060.  
  1061. c          -on a a minima 2 points confondus
  1062.             if(IP2.ne.0) then
  1063. c             -on est sur le bord du triangle => TA decoupe en 2
  1064.               if(IP3.ne.0) then
  1065.                 ICOUPA=21
  1066. c             -on a 2 points confondus => TA decoupe en 3
  1067.               else
  1068.                 ICOUPA=39
  1069.               endif
  1070.               call AJOUPO(xQ2,IPT3,XELIM,IQQ)
  1071. c          -TA decoupe en 4
  1072.             else
  1073.               call AJOUPO(xP2,IPT3,XELIM,IP2)
  1074.               call AJOUPO(xQ2,IPT3,XELIM,IQQ)
  1075.               ICOUPA=42
  1076.             endif
  1077. c          +on crée le segment intersection
  1078.               NBTMP=NBELE3+1
  1079.               NODE(1)=IQQ
  1080.               NODE(2)=IP2
  1081.               call ajouel(NODE,2,2,IPT3,NBTMP)
  1082.               NBELE3=max(NBELE3,NBTMP)
  1083. c             endif
  1084.  
  1085. c     --- Cas chevauchement selon [Q3 P2] ---
  1086.           ELSEIF (((sP3+XELIM).le.sQ3.and.sQ3.le.(sP2-XELIM)
  1087.      &              .and.sQ3.le.sQ2)
  1088.      &        .OR.((sP3-XELIM).ge.sQ3.and.sQ3.ge.(sP2+XELIM)
  1089.      &              .and.sQ3.ge.sQ2)) THEN
  1090.  
  1091. c          -on a a minima 2 points confondus
  1092.             if(IP2.ne.0) then
  1093. c             -on est sur le bord du triangle => TA decoupe en 2
  1094.               if(IP3.ne.0) then
  1095.                 ICOUPA=21
  1096. c             -on a 2 points confondus => TA decoupe en 3
  1097.               else
  1098.                 ICOUPA=39
  1099.               endif
  1100.               call AJOUPO(xQ3,IPT3,XELIM,IQQ)
  1101. c          -TA decoupe en 4
  1102.             else
  1103.               call AJOUPO(xP2,IPT3,XELIM,IP2)
  1104.               call AJOUPO(xQ3,IPT3,XELIM,IQQ)
  1105.               ICOUPA=42
  1106.             endif
  1107. c          +on crée le segment intersection
  1108.               NBTMP=NBELE3+1
  1109.               NODE(1)=IP2
  1110.               NODE(2)=IQQ
  1111.               call ajouel(NODE,2,2,IPT3,NBTMP)
  1112.               NBELE3=max(NBELE3,NBTMP)
  1113. c             endif
  1114.  
  1115.           ENDIF
  1116.           if(iimpi.ge.10) then
  1117.             write(ioimp,fmt="('IPT3(1et2',I5,')=',I5,I5)")
  1118.      &      NBELE3,IPT3.NUM(1,NBELE3),IPT3.NUM(2,NBELE3)
  1119.           endif
  1120.  
  1121.  
  1122. C-----------------------------------------------------------------------
  1123. C         PREPARATION DE LA CONSTRUCTION DES FILS DE TB
  1124. C         + CONSTRUCTION DU SEGMENT D'INTERSECTION
  1125.  
  1126. c     --- Cas [Q2-Q3] inclus dans [P2-P3] => segment Q2-Q3 ---
  1127.           IF (sminP.lt.(sminQ+XELIM).and.(smaxQ-XELIM).lt.smaxP) THEN
  1128.  
  1129. c          -on est sur le bord du triangle => on ne decoupe pas B
  1130.             if(IQ2.ne.0.and.IQ3.ne.0) then
  1131.               ICOUPB=0
  1132. c          -on a 2 points confondus => TB decoupe en 2
  1133.             elseif(IQ2.ne.0) then
  1134.               call AJOUPO(xQ3,IPT4,XELIM,IQ3)
  1135.               ICOUPB=23
  1136.             elseif(IQ3.ne.0) then
  1137.               call AJOUPO(xQ2,IPT4,XELIM,IQ2)
  1138.               IQ3 = ITB(IB3)
  1139.               ICOUPB=22
  1140. c          -on est dans le cas general => TB decoupe en 3
  1141.             else
  1142. c             on ajoute Q2 et Q3 a XCOOR
  1143.               call AJOUPO(xQ2,IPT4,XELIM,IQ2)
  1144.               call AJOUPO(xQ3,IPT4,XELIM,IQ3)
  1145.               ICOUPB=30
  1146.             endif
  1147. c          +on crée le segment intersection
  1148.             if(.not.SELF) then
  1149.               NBTMP=NBELE4+1
  1150.               NODE(1)=IQ2
  1151.               NODE(2)=IQ3
  1152.               call ajouel(NODE,2,2,IPT4,NBTMP)
  1153.               NBELE4=max(NBELE4,NBTMP)
  1154.             endif
  1155.  
  1156. c     --- Cas [Q2-Q3] exclus de [P2-P3] => segment P2-P3 ---
  1157.           ELSEIF(sminP.gt.(sminQ+XELIM).and.(smaxQ-XELIM).gt.smaxP)THEN
  1158.  
  1159. c          -on est sur le bord du triangle => TB decoupe en 3
  1160.             if(IQ2.ne.0.and.IQ3.ne.0) then
  1161.               ICOUPB=31
  1162. c          -on a 2 points confondus => TB decoupe en 5
  1163.             else
  1164.               if(abs(xddB2).ge.abs(xddB3)) then
  1165.                 ICOUPB=52
  1166.               else
  1167.                 ICOUPB=53
  1168.               endif
  1169.             endif
  1170.               call AJOUPO(xP2,IPT4,XELIM,IPP2)
  1171. c             cas particulier IP2 = IP3
  1172.               if(IP2.eq.IP3) then
  1173.                 IPP3=IPP2
  1174.               else
  1175.                 call AJOUPO(xP3,IPT4,XELIM,IPP3)
  1176.               endif
  1177. c           petite astuce pour garder l orientation
  1178.             if ( (sQ2.lt.sQ3.and.sP2.lt.sP3)
  1179.      &       .or.(sQ2.gt.sQ3.and.sP2.gt.sP3) ) then
  1180.               IPP2=IPP2
  1181.               IPP3=IPP3
  1182.             else
  1183.               IPPtmp=IPP3
  1184.               IPP3=IPP2
  1185.               IPP2=IPPtmp
  1186.             endif
  1187. c          +on crée le segment intersection
  1188.             if(.not.SELF) then
  1189.               NBTMP=NBELE4+1
  1190.               NODE(1)=IPP2
  1191.               NODE(2)=IPP3
  1192.               call ajouel(NODE,2,2,IPT4,NBTMP)
  1193.               NBELE4=max(NBELE4,NBTMP)
  1194.             endif
  1195.  
  1196. c     --- Cas chevauchement selon [P2 Q3] ---
  1197.           ELSEIF (((sQ2+XELIM).le.sP2.and.sP2.le.(sQ3-XELIM)
  1198.      &              .and.sP2.le.sP3)
  1199.      &        .OR.((sQ2-XELIM).ge.sP2.and.sP2.ge.(sQ3+XELIM)
  1200.      &              .and.sP2.ge.sP3)) THEN
  1201.  
  1202. c          -on a a minima 2 points confondus
  1203.             if(IQ3.ne.0) then
  1204.               IQ3 = ITB(IB3)
  1205. c             -on est sur le bord du triangle => TB decoupe en 2
  1206.               if(IQ2.ne.0) then
  1207.                 ICOUPB=21
  1208. c             -on a 2 points confondus => TB decoupe en 3
  1209.               else
  1210.                 ICOUPB=39
  1211.               endif
  1212.               call AJOUPO(xP2,IPT4,XELIM,IPP)
  1213. c          -TB decoupe en 4
  1214.             else
  1215.               call AJOUPO(xQ3,IPT4,XELIM,IQ3)
  1216.               call AJOUPO(xP2,IPT4,XELIM,IPP)
  1217.              ICOUPB=43
  1218.             endif
  1219. c          +on crée le segment intersection
  1220.             if(.not.SELF) then
  1221.               NBTMP=NBELE4+1
  1222.               NODE(1)=IPP
  1223.               NODE(2)=IQ3
  1224.               call ajouel(NODE,2,2,IPT4,NBTMP)
  1225.               NBELE4=max(NBELE4,NBTMP)
  1226.             endif
  1227.  
  1228. c     --- Cas chevauchement selon [P3 Q3] ---
  1229.           ELSEIF (((sQ2+XELIM).le.sP3.and.sP3.le.(sQ3-XELIM)
  1230.      &              .and.sP3.le.sP2)
  1231.      &        .OR.((sQ2-XELIM).ge.sP3.and.sP3.ge.(sQ3+XELIM)
  1232.      &              .and.sP3.ge.sP2)) THEN
  1233.  
  1234. c          -on a a minima 2 points confondus
  1235.             if(IQ3.ne.0) then
  1236.               IQ3 = ITB(IB3)
  1237. c             -on est sur le bord du triangle => TB decoupe en 2
  1238.               if(IQ2.ne.0) then
  1239.                 ICOUPB=21
  1240. c             -on a 2 points confondus => TB decoupe en 3
  1241.               else
  1242.                 ICOUPB=39
  1243.               endif
  1244.               call AJOUPO(xP3,IPT4,XELIM,IPP)
  1245. c          -TB decoupe en 4
  1246.             else
  1247.               call AJOUPO(xQ3,IPT4,XELIM,IQ3)
  1248.               call AJOUPO(xP3,IPT4,XELIM,IPP)
  1249.               ICOUPB=43
  1250.             endif
  1251. c          +on crée le segment intersection
  1252.             if(.not.SELF) then
  1253.               NBTMP=NBELE4+1
  1254.               NODE(1)=IPP
  1255.               NODE(2)=IQ3
  1256.               call ajouel(NODE,2,2,IPT4,NBTMP)
  1257.               NBELE4=max(NBELE4,NBTMP)
  1258.             endif
  1259.  
  1260. c     --- Cas chevauchement selon [P2 Q2] ---
  1261.           ELSEIF (((sQ3+XELIM).le.sP2.and.sP2.le.(sQ2-XELIM)
  1262.      &              .and.sP2.le.sP3)
  1263.      &        .OR.((sQ3-XELIM).ge.sP2.and.sP2.ge.(sQ2+XELIM)
  1264.      &              .and.sP2.ge.sP3)) THEN
  1265.  
  1266. c          -on a a minima 2 points confondus
  1267.             if(IQ2.ne.0) then
  1268. c             -on est sur le bord du triangle => TB decoupe en 2
  1269.               if(IQ3.ne.0) then
  1270.                 ICOUPB=21
  1271. c             -on a 2 points confondus => TB decoupe en 3
  1272.               else
  1273.                 ICOUPB=39
  1274.               endif
  1275.               call AJOUPO(xP2,IPT4,XELIM,IPP)
  1276. c          -TB decoupe en 4
  1277.             else
  1278.               call AJOUPO(xQ2,IPT4,XELIM,IQ2)
  1279.               call AJOUPO(xP2,IPT4,XELIM,IPP)
  1280.               ICOUPB=42
  1281.             endif
  1282. c          +on crée le segment intersection
  1283.             if(.not.SELF) then
  1284.               NBTMP=NBELE4+1
  1285.               NODE(1)=IQ2
  1286.               NODE(2)=IPP
  1287.               call ajouel(NODE,2,2,IPT4,NBTMP)
  1288.               NBELE4=max(NBELE4,NBTMP)
  1289.             endif
  1290.  
  1291. c     --- Cas chevauchement selon [P3 Q2] ---
  1292.           ELSEIF (((sQ3+XELIM).le.sP3.and.sP3.le.(sQ2-XELIM)
  1293.      &              .and.sP3.le.sP2)
  1294.      &        .OR.((sQ3-XELIM).ge.sP3.and.sP3.ge.(sQ2+XELIM)
  1295.      &              .and.sP3.ge.sP2)) THEN
  1296.  
  1297. c          -on a a minima 2 points confondus
  1298.             if(IQ2.ne.0) then
  1299. c             -on est sur le bord du triangle => TB decoupe en 2
  1300.               if(IQ3.ne.0) then
  1301.                 ICOUPB=21
  1302. c             -on a 2 points confondus => TB decoupe en 3
  1303.               else
  1304.                 ICOUPB=39
  1305.               endif
  1306.               call AJOUPO(xP3,IPT4,XELIM,IPP)
  1307. c          -TB decoupe en 4
  1308.             else
  1309.               call AJOUPO(xQ2,IPT4,XELIM,IQ2)
  1310.               call AJOUPO(xP3,IPT4,XELIM,IPP)
  1311.               ICOUPB=42
  1312.             endif
  1313. c          +on crée le segment intersection
  1314.             if(.not.SELF) then
  1315.               NBTMP=NBELE4+1
  1316.               NODE(1)=IQ2
  1317.               NODE(2)=IPP
  1318.               call ajouel(NODE,2,2,IPT4,NBTMP)
  1319.               NBELE4=max(NBELE4,NBTMP)
  1320.             endif
  1321.  
  1322.           ENDIF
  1323.  
  1324.  
  1325.           if(ICOUPA.eq.0.and.ICOUPB.eq.0) goto 1900
  1326. C-----------------------------------------------------------------------
  1327. C         CREATION DU SEGMENT INTERSECTION (fait + haut)
  1328. c         CREATION DES TRIANGLES DECOUPES ISSUS DE TA
  1329.           IF    (ICOUPA.EQ.21) THEN
  1330. c           on decoupe TA en 2 triangles (rem: ITA(IA3)=IP3)
  1331. c          -1er sous triangle de TA
  1332.             IPT1.NUM(1,IA)=ITA(IA1)
  1333.             IPT1.NUM(2,IA)=ITA(IA2)
  1334.             IPT1.NUM(3,IA)=IQQ
  1335. c          -2eme sous triangle de TA
  1336.             IPT1.NUM(1, NBELE1+1)=IQQ
  1337.             IPT1.NUM(2, NBELE1+1)=ITA(IA3)
  1338.             IPT1.NUM(3, NBELE1+1)=ITA(IA1)
  1339.             IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1340.             NBELE1  = NBELE1+1
  1341.           ELSEIF(ICOUPA.EQ.22) THEN
  1342. c           on decoupe TA en 2 triangles (rem: ITA(IA3)=IP3)
  1343. c          -1er sous triangle de TA
  1344.             IPT1.NUM(1,IA)=ITA(IA1)
  1345.             IPT1.NUM(2,IA)=IP2
  1346.             IPT1.NUM(3,IA)=ITA(IA3)
  1347. c          -2eme sous triangle de TA
  1348.             IPT1.NUM(1, NBELE1+1)=IP2
  1349.             IPT1.NUM(2, NBELE1+1)=ITA(IA2)
  1350.             IPT1.NUM(3, NBELE1+1)=ITA(IA3)
  1351.             IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1352.             NBELE1  = NBELE1+1
  1353.           ELSEIF(ICOUPA.EQ.23) THEN
  1354. c           on decoupe TA en 2 triangles (rem: ITA(IA2)=IP2)
  1355. c          -1er sous triangle de TA
  1356.             IPT1.NUM(1,IA)=ITA(IA1)
  1357.             IPT1.NUM(2,IA)=ITA(IA2)
  1358.             IPT1.NUM(3,IA)=IP3
  1359. c          -2eme sous triangle de TA
  1360.             IPT1.NUM(1, NBELE1+1)=ITA(IA2)
  1361.             IPT1.NUM(2, NBELE1+1)=ITA(IA3)
  1362.             IPT1.NUM(3, NBELE1+1)=IP3
  1363.             IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1364.             NBELE1  = NBELE1+1
  1365.           ELSEIF(ICOUPA.EQ.30) THEN
  1366. c           on decoupe TA en 3 triangles
  1367. c          -1er sous triangle de TA
  1368.             IPT1.NUM(1,IA)=ITA(IA1)
  1369.             IPT1.NUM(2,IA)=IP2
  1370.             IPT1.NUM(3,IA)=IP3
  1371. c          -2eme sous triangle de TA
  1372.             IPT1.NUM(1, NBELE1+1)=IP2
  1373.             IPT1.NUM(2, NBELE1+1)=ITA(IA2)
  1374.             IPT1.NUM(3, NBELE1+1)=ITA(IA3)
  1375.             IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1376. c          -3eme sous triangle de TA
  1377.             IPT1.NUM(1, NBELE1+2)=ITA(IA3)
  1378.             IPT1.NUM(2, NBELE1+2)=IP3
  1379.             IPT1.NUM(3, NBELE1+2)=IP2
  1380.             IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1381.             NBELE1  = NBELE1+2
  1382.           ELSEIF(ICOUPA.EQ.31) THEN
  1383. c           on decoupe TA en 3 triangles
  1384. c          -1er sous triangle de TA
  1385.             IPT1.NUM(1,IA)=ITA(IA1)
  1386.             IPT1.NUM(2,IA)=ITA(IA2)
  1387.             IPT1.NUM(3,IA)=IQQ2
  1388. c          -2eme sous triangle de TA
  1389.             IPT1.NUM(1, NBELE1+1)=ITA(IA1)
  1390.             IPT1.NUM(2, NBELE1+1)=IQQ2
  1391.             IPT1.NUM(3, NBELE1+1)=IQQ3
  1392.             IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1393. c          -3eme sous triangle de TA
  1394.             IPT1.NUM(1, NBELE1+2)=ITA(IA1)
  1395.             IPT1.NUM(2, NBELE1+2)=IQQ3
  1396.             IPT1.NUM(3, NBELE1+2)=ITA(IA3)
  1397.             IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1398.             NBELE1  = NBELE1+2
  1399.           ELSEIF(ICOUPA.EQ.39) THEN
  1400. c           on decoupe TA en 3 triangles
  1401. c          -1er sous triangle de TA
  1402.             IPT1.NUM(1,IA)=ITA(IA1)
  1403.             IPT1.NUM(2,IA)=ITA(IA2)
  1404.             IPT1.NUM(3,IA)=IQQ
  1405. c          -2eme sous triangle de TA
  1406.             IPT1.NUM(1, NBELE1+1)=ITA(IA2)
  1407.             IPT1.NUM(2, NBELE1+1)=IQQ
  1408.             IPT1.NUM(3, NBELE1+1)=ITA(IA3)
  1409.             IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1410. c          -3eme sous triangle de TA
  1411.             IPT1.NUM(1, NBELE1+2)=ITA(IA3)
  1412.             IPT1.NUM(2, NBELE1+2)=IQQ
  1413.             IPT1.NUM(3, NBELE1+2)=ITA(IA1)
  1414.             IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1415.             NBELE1  = NBELE1+2
  1416.           ELSEIF(ICOUPA.EQ.42) THEN
  1417. c           on decoupe TA en 4 triangles (arete 1-2 coupee)
  1418. c          -1er sous triangle de TA
  1419.             IPT1.NUM(1,IA)=ITA(IA1)
  1420.             IPT1.NUM(2,IA)=IP2
  1421.             IPT1.NUM(3,IA)=IQQ
  1422. c          -2eme sous triangle de TA
  1423.             IPT1.NUM(1, NBELE1+1)=IP2
  1424.             IPT1.NUM(2, NBELE1+1)=ITA(IA2)
  1425.             IPT1.NUM(3, NBELE1+1)=IQQ
  1426.             IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1427. c          -3eme sous triangle de TA
  1428.             IPT1.NUM(1, NBELE1+2)=ITA(IA2)
  1429.             IPT1.NUM(2, NBELE1+2)=ITA(IA3)
  1430.             IPT1.NUM(3, NBELE1+2)=IQQ
  1431.             IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1432. c          -4eme sous triangle de TA
  1433.             IPT1.NUM(1, NBELE1+3)=ITA(IA3)
  1434.             IPT1.NUM(2, NBELE1+3)=ITA(IA1)
  1435.             IPT1.NUM(3, NBELE1+3)=IQQ
  1436.             IPT1.ICOLOR(NBELE1+3)=IPT1.ICOLOR(IA)
  1437.             NBELE1  = NBELE1+3
  1438.           ELSEIF(ICOUPA.EQ.43) THEN
  1439. c           on decoupe TA en 4 triangles (arete 1-3 coupee)
  1440. c          -1er sous triangle de TA
  1441.             IPT1.NUM(1,IA)=ITA(IA1)
  1442.             IPT1.NUM(2,IA)=ITA(IA2)
  1443.             IPT1.NUM(3,IA)=IQQ
  1444. c          -2eme sous triangle de TA
  1445.             IPT1.NUM(1, NBELE1+1)=ITA(IA2)
  1446.             IPT1.NUM(2, NBELE1+1)=ITA(IA3)
  1447.             IPT1.NUM(3, NBELE1+1)=IQQ
  1448.             IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1449. c          -3eme sous triangle de TA
  1450.             IPT1.NUM(1, NBELE1+2)=ITA(IA3)
  1451.             IPT1.NUM(2, NBELE1+2)=IP3
  1452.             IPT1.NUM(3, NBELE1+2)=IQQ
  1453.             IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1454. c          -4eme sous triangle de TA
  1455.             IPT1.NUM(1, NBELE1+3)=IP3
  1456.             IPT1.NUM(2, NBELE1+3)=ITA(IA1)
  1457.             IPT1.NUM(3, NBELE1+3)=IQQ
  1458.             IPT1.ICOLOR(NBELE1+3)=IPT1.ICOLOR(IA)
  1459.             NBELE1  = NBELE1+3
  1460.           ELSEIF(ICOUPA.EQ.52) THEN
  1461. c           on decoupe TA en 5 triangles
  1462. c          -1er sous triangle de TA
  1463.             IPT1.NUM(1,IA)=ITA(IA1)
  1464.             IPT1.NUM(2,IA)=IQQ2
  1465.             IPT1.NUM(3,IA)=IQQ3
  1466. c          -2eme sous triangle de TA
  1467.             IPT1.NUM(1,NBELE1+1)=ITA(IA1)
  1468.             IPT1.NUM(2,NBELE1+1)=ITA(IA2)
  1469.             IPT1.NUM(3,NBELE1+1)=IQQ2
  1470.             IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1471. c          -3eme sous triangle de TA
  1472.             IPT1.NUM(1,NBELE1+2)=ITA(IA2)
  1473.             IPT1.NUM(2,NBELE1+2)=IQQ3
  1474.             IPT1.NUM(3,NBELE1+2)=IQQ2
  1475.             IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1476. c          -4eme sous triangle de TA
  1477.             IPT1.NUM(1,NBELE1+3)=ITA(IA2)
  1478.             IPT1.NUM(2,NBELE1+3)=ITA(IA3)
  1479.             IPT1.NUM(3,NBELE1+3)=IQQ3
  1480.             IPT1.ICOLOR(NBELE1+3)=IPT1.ICOLOR(IA)
  1481. c          -5eme sous triangle de TA
  1482.             IPT1.NUM(1,NBELE1+4)=ITA(IA3)
  1483.             IPT1.NUM(2,NBELE1+4)=ITA(IA1)
  1484.             IPT1.NUM(3,NBELE1+4)=IQQ3
  1485.             IPT1.ICOLOR(NBELE1+4)=IPT1.ICOLOR(IA)
  1486.             NBELE1  = NBELE1+4
  1487.           ELSEIF(ICOUPA.EQ.53) THEN
  1488. c           on decoupe TA en 5 triangles
  1489. c          -1er sous triangle de TA
  1490.             IPT1.NUM(1,IA)=ITA(IA1)
  1491.             IPT1.NUM(2,IA)=IQQ2
  1492.             IPT1.NUM(3,IA)=IQQ3
  1493. c          -2eme sous triangle de TA
  1494.             IPT1.NUM(1,NBELE1+1)=ITA(IA1)
  1495.             IPT1.NUM(2,NBELE1+1)=ITA(IA2)
  1496.             IPT1.NUM(3,NBELE1+1)=IQQ2
  1497.             IPT1.ICOLOR(NBELE1+1)=IPT1.ICOLOR(IA)
  1498. c          -3eme sous triangle de TA
  1499.             IPT1.NUM(1,NBELE1+2)=ITA(IA2)
  1500.             IPT1.NUM(2,NBELE1+2)=ITA(IA3)
  1501.             IPT1.NUM(3,NBELE1+2)=IQQ2
  1502.             IPT1.ICOLOR(NBELE1+2)=IPT1.ICOLOR(IA)
  1503. c          -4eme sous triangle de TA
  1504.             IPT1.NUM(1,NBELE1+3)=ITA(IA3)
  1505.             IPT1.NUM(2,NBELE1+3)=IQQ3
  1506.             IPT1.NUM(3,NBELE1+3)=IQQ2
  1507.             IPT1.ICOLOR(NBELE1+3)=IPT1.ICOLOR(IA)
  1508. c          -5eme sous triangle de TA
  1509.             IPT1.NUM(1,NBELE1+4)=ITA(IA3)
  1510.             IPT1.NUM(2,NBELE1+4)=ITA(IA1)
  1511.             IPT1.NUM(3,NBELE1+4)=IQQ3
  1512.             IPT1.ICOLOR(NBELE1+4)=IPT1.ICOLOR(IA)
  1513.             NBELE1  = NBELE1+4
  1514.           ENDIF
  1515.  
  1516. c         CREATION DES TRIANGLES DECOUPES ISSUS DE TB
  1517.           if(SELF) NBELE2=NBELE1
  1518.           IF    (ICOUPB.EQ.21) THEN
  1519. c           on decoupe TB en 2 triangles
  1520. c          -1er sous triangle de TB
  1521.             IPT2.NUM(1,IB)=ITB(IB1)
  1522.             IPT2.NUM(2,IB)=ITB(IB2)
  1523.             IPT2.NUM(3,IB)=IPP
  1524. c          -2eme sous triangle de TB
  1525.             IPT2.NUM(1, NBELE2+1)=IPP
  1526.             IPT2.NUM(2, NBELE2+1)=ITB(IB3)
  1527.             IPT2.NUM(3, NBELE2+1)=ITB(IB1)
  1528.             IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1529.             NBELE2  = NBELE2+1
  1530.           ELSEIF(ICOUPB.EQ.22) THEN
  1531. c           on decoupe TB en 2 triangles (rem: ITB(IB3)=IQ3)
  1532. c          -1er sous triangle de TB
  1533.             IPT2.NUM(1,IB)=ITB(IB1)
  1534.             IPT2.NUM(2,IB)=IQ2
  1535.             IPT2.NUM(3,IB)=ITB(IB3)
  1536. c          -2eme sous triangle de TB
  1537.             IPT2.NUM(1, NBELE2+1)=IQ2
  1538.             IPT2.NUM(2, NBELE2+1)=ITB(IB2)
  1539.             IPT2.NUM(3, NBELE2+1)=ITB(IB3)
  1540.             IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1541.             NBELE2  = NBELE2+1
  1542.           ELSEIF(ICOUPB.EQ.23) THEN
  1543. c           on decoupe TB en 2 triangles (rem: ITB(IB2)=IQ2)
  1544. c          -1er sous triangle de TB
  1545.             IPT2.NUM(1,IB)=ITB(IB1)
  1546.             IPT2.NUM(2,IB)=ITB(IB2)
  1547.             IPT2.NUM(3,IB)=IQ3
  1548. c          -2eme sous triangle de TB
  1549.             IPT2.NUM(1, NBELE2+1)=ITB(IB2)
  1550.             IPT2.NUM(2, NBELE2+1)=ITB(IB3)
  1551.             IPT2.NUM(3, NBELE2+1)=IQ3
  1552.             IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1553.             NBELE2  = NBELE2+1
  1554.           ELSEIF(ICOUPB.EQ.30) THEN
  1555. c           on decoupe TB en 3 triangles
  1556. c          -1er sous triangle de TB
  1557.             IPT2.NUM(1,IB)=ITB(IB1)
  1558.             IPT2.NUM(2,IB)=IQ2
  1559.             IPT2.NUM(3,IB)=IQ3
  1560. c          -2eme sous triangle de TB
  1561.             IPT2.NUM(1,NBELE2+1)=IQ2
  1562.             IPT2.NUM(2,NBELE2+1)=ITB(IB2)
  1563.             IPT2.NUM(3,NBELE2+1)=ITB(IB3)
  1564.             IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1565. c          -3eme sous triangle de TB
  1566.             IPT2.NUM(1,NBELE2+2)=IQ2
  1567.             IPT2.NUM(2,NBELE2+2)=ITB(IB3)
  1568.             IPT2.NUM(3,NBELE2+2)=IQ3
  1569.             IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1570.             NBELE2=NBELE2+2
  1571.           ELSEIF(ICOUPB.EQ.31) THEN
  1572. c           on decoupe TB en 3 triangles
  1573. c          -1er sous triangle de TB
  1574.             IPT2.NUM(1,IB)=ITB(IB1)
  1575.             IPT2.NUM(2,IB)=ITB(IB2)
  1576.             IPT2.NUM(3,IB)=IPP2
  1577. c          -2eme sous triangle de TB
  1578.             IPT2.NUM(1, NBELE2+1)=ITB(IB1)
  1579.             IPT2.NUM(2, NBELE2+1)=IPP2
  1580.             IPT2.NUM(3, NBELE2+1)=IPP3
  1581.             IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1582. c          -3eme sous triangle de TB
  1583.             IPT2.NUM(1, NBELE2+2)=ITB(IB1)
  1584.             IPT2.NUM(2, NBELE2+2)=IPP3
  1585.             IPT2.NUM(3, NBELE2+2)=ITB(IB3)
  1586.             IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1587.             NBELE2=NBELE2+2
  1588.           ELSEIF(ICOUPB.EQ.39) THEN
  1589. c           on decoupe TB en 3 triangles
  1590. c          -1er sous triangle de TB
  1591.             IPT2.NUM(1,IB)=ITB(IB1)
  1592.             IPT2.NUM(2,IB)=ITB(IB2)
  1593.             IPT2.NUM(3,IB)=IPP
  1594. c          -2eme sous triangle de TB
  1595.             IPT2.NUM(1, NBELE2+1)=ITB(IB2)
  1596.             IPT2.NUM(2, NBELE2+1)=IPP
  1597.             IPT2.NUM(3, NBELE2+1)=ITB(IB3)
  1598.             IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1599. c          -3eme sous triangle de TB
  1600.             IPT2.NUM(1, NBELE2+2)=ITB(IB3)
  1601.             IPT2.NUM(2, NBELE2+2)=IPP
  1602.             IPT2.NUM(3, NBELE2+2)=ITB(IB1)
  1603.             IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1604.             NBELE2=NBELE2+2
  1605.           ELSEIF(ICOUPB.EQ.42) THEN
  1606. c           on decoupe TB en 4 triangles (arete 1-2 coupee)
  1607. c          -1er sous triangle de TB
  1608.             IPT2.NUM(1,IB)=ITB(IB1)
  1609.             IPT2.NUM(2,IB)=IQ2
  1610.             IPT2.NUM(3,IB)=IPP
  1611. c          -2eme sous triangle de TB
  1612.             IPT2.NUM(1, NBELE2+1)=IQ2
  1613.             IPT2.NUM(2, NBELE2+1)=ITB(IB2)
  1614.             IPT2.NUM(3, NBELE2+1)=IPP
  1615.             IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1616. c          -3eme sous triangle de TB
  1617.             IPT2.NUM(1, NBELE2+2)=ITB(IB2)
  1618.             IPT2.NUM(2, NBELE2+2)=ITB(IB3)
  1619.             IPT2.NUM(3, NBELE2+2)=IPP
  1620.             IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1621. c          -4eme sous triangle de TB
  1622.             IPT2.NUM(1, NBELE2+3)=ITB(IB3)
  1623.             IPT2.NUM(2, NBELE2+3)=ITB(IB1)
  1624.             IPT2.NUM(3, NBELE2+3)=IPP
  1625.             IPT2.ICOLOR(NBELE2+3)=IPT2.ICOLOR(IB)
  1626.             NBELE2=NBELE2+3
  1627.           ELSEIF(ICOUPB.EQ.43) THEN
  1628. c           on decoupe TB en 4 triangles (arete 1-3 coupee)
  1629. c          -1er sous triangle de TB
  1630.             IPT2.NUM(1,IB)=ITB(IB1)
  1631.             IPT2.NUM(2,IB)=ITB(IB2)
  1632.             IPT2.NUM(3,IB)=IPP
  1633. c          -2eme sous triangle de TB
  1634.             IPT2.NUM(1, NBELE2+1)=ITB(IB2)
  1635.             IPT2.NUM(2, NBELE2+1)=ITB(IB3)
  1636.             IPT2.NUM(3, NBELE2+1)=IPP
  1637.             IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1638. c          -3eme sous triangle de TB
  1639.             IPT2.NUM(1, NBELE2+2)=ITB(IB3)
  1640.             IPT2.NUM(2, NBELE2+2)=IQ3
  1641.             IPT2.NUM(3, NBELE2+2)=IPP
  1642.             IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1643. c          -4eme sous triangle de TB
  1644.             IPT2.NUM(1, NBELE2+3)=IQ3
  1645.             IPT2.NUM(2, NBELE2+3)=ITB(IB1)
  1646.             IPT2.NUM(3, NBELE2+3)=IPP
  1647.             IPT2.ICOLOR(NBELE2+3)=IPT2.ICOLOR(IB)
  1648.             NBELE2=NBELE2+3
  1649.           ELSEIF(ICOUPB.EQ.52) THEN
  1650. c           on decoupe TB en 5 triangles
  1651. c          -1er sous triangle de TB
  1652.             IPT2.NUM(1,IB)=ITB(IB1)
  1653.             IPT2.NUM(2,IB)=IPP2
  1654.             IPT2.NUM(3,IB)=IPP3
  1655. c          -2eme sous triangle de TB
  1656.             IPT2.NUM(1,NBELE2+1)=ITB(IB1)
  1657.             IPT2.NUM(2,NBELE2+1)=ITB(IB2)
  1658.             IPT2.NUM(3,NBELE2+1)=IPP2
  1659.             IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1660. c          -3eme sous triangle de TB
  1661.             IPT2.NUM(1,NBELE2+2)=ITB(IB2)
  1662.             IPT2.NUM(2,NBELE2+2)=IPP3
  1663.             IPT2.NUM(3,NBELE2+2)=IPP2
  1664.             IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1665. c          -4eme sous triangle de TB
  1666.             IPT2.NUM(1,NBELE2+3)=ITB(IB2)
  1667.             IPT2.NUM(2,NBELE2+3)=ITB(IB3)
  1668.             IPT2.NUM(3,NBELE2+3)=IPP3
  1669.             IPT2.ICOLOR(NBELE2+3)=IPT2.ICOLOR(IB)
  1670. c          -5eme sous triangle de TB
  1671.             IPT2.NUM(1,NBELE2+4)=ITB(IB3)
  1672.             IPT2.NUM(2,NBELE2+4)=ITB(IB1)
  1673.             IPT2.NUM(3,NBELE2+4)=IPP3
  1674.             IPT2.ICOLOR(NBELE2+4)=IPT2.ICOLOR(IB)
  1675.             NBELE2=NBELE2+4
  1676.           ELSEIF(ICOUPB.EQ.53) THEN
  1677. c           on decoupe TB en 5 triangles
  1678. c          -1er sous triangle de TB
  1679.             IPT2.NUM(1,IB)=ITB(IB1)
  1680.             IPT2.NUM(2,IB)=IPP2
  1681.             IPT2.NUM(3,IB)=IPP3
  1682. c          -2eme sous triangle de TB
  1683.             IPT2.NUM(1,NBELE2+1)=ITB(IB1)
  1684.             IPT2.NUM(2,NBELE2+1)=ITB(IB2)
  1685.             IPT2.NUM(3,NBELE2+1)=IPP2
  1686.             IPT2.ICOLOR(NBELE2+1)=IPT2.ICOLOR(IB)
  1687. c          -3eme sous triangle de TB
  1688.             IPT2.NUM(1,NBELE2+2)=ITB(IB2)
  1689.             IPT2.NUM(2,NBELE2+2)=ITB(IB3)
  1690.             IPT2.NUM(3,NBELE2+2)=IPP2
  1691.             IPT2.ICOLOR(NBELE2+2)=IPT2.ICOLOR(IB)
  1692. c          -4eme sous triangle de TB
  1693.             IPT2.NUM(1,NBELE2+3)=ITB(IB3)
  1694.             IPT2.NUM(2,NBELE2+3)=IPP3
  1695.             IPT2.NUM(3,NBELE2+3)=IPP2
  1696.             IPT2.ICOLOR(NBELE2+3)=IPT2.ICOLOR(IB)
  1697. c          -5eme sous triangle de TB
  1698.             IPT2.NUM(1,NBELE2+4)=ITB(IB3)
  1699.             IPT2.NUM(2,NBELE2+4)=ITB(IB1)
  1700.             IPT2.NUM(3,NBELE2+4)=IPP3
  1701.             IPT2.ICOLOR(NBELE2+4)=IPT2.ICOLOR(IB)
  1702.             NBELE2=NBELE2+4
  1703.           ENDIF
  1704.           if(SELF) NBELE1=NBELE2
  1705.  
  1706.           if(iimpi.ge.10) then
  1707.           write(ioimp,*) '------------------------------->',ICOUPA
  1708.           write(ioimp,*) '------------------------------->',ICOUPB
  1709.           write(ioimp,*) '-----------------------------------------'
  1710.           endif
  1711.  
  1712. c         Tant qu'il y a une decoupe, on reprend le traitement relatif
  1713. c         au 1er element "fils" du triangle decoupe ...
  1714.           IBDEB0 = IB+1
  1715.           if(IBDEB0.le.NBELE2) GOTO 1100
  1716.  
  1717.  
  1718. c         1900 : label ou l on arrive s il n y a pas eu de decoupe
  1719.  1900     CONTINUE
  1720.           IBDEB0=0
  1721.           if(IA.ge.NBELE1.and.IB.ge.NBELE2) goto 9000
  1722.  
  1723.  
  1724.  2000   CONTINUE
  1725. C====== fin de BOUCLE SUR LES ELEMENTS B ===============================
  1726.  
  1727.  
  1728.   900     CONTINUE
  1729.           if(IA.ge.NBELE1) goto 9000
  1730.  
  1731.  1000 CONTINUE
  1732. C==== fin de BOUCLE SUR LES ELEMENTS A =================================
  1733.  
  1734.  
  1735.  
  1736.  9000 CONTINUE
  1737. C-----------------------------------------------------------------------
  1738. C     FIN NORMALE
  1739.  
  1740. c     on supprime les segments de travaux locaux
  1741.       segsup,WORK1
  1742.  
  1743. c     on ajuste IPT1
  1744.       NBNN  =NBNN1
  1745.       NBELEM=NBELE1
  1746.       segadj,IPT1
  1747. c     on ajuste IPT3
  1748.       NBNN = 2
  1749.       NBELEM = NBELE3
  1750.       segadj,IPT3
  1751. c     on desactive
  1752.       segdes,IPT1,IPT3
  1753.  
  1754.       IF(.not.SELF) THEN
  1755.         segsup,WORK2
  1756. c       on ajuste IPT2
  1757.         NBNN  =NBNN2
  1758.         NBELEM=NBELE2
  1759.         segadj,IPT2
  1760. c       on ajuste IPT4
  1761.         NBNN = 2
  1762.         NBELEM = NBELE4
  1763.         segadj,IPT4
  1764. c       on desactive
  1765.         segdes,IPT2,IPT4
  1766.       ENDIF
  1767.  
  1768.  
  1769. C     ON REGENERE et on ne garde que les triangles (ITYPEL=4)
  1770.       CALL REGENE(IPT1)
  1771.       IPT1ok = IPT1
  1772.       segact,IPT1
  1773.       NBSOU1=IPT1.LISOUS(/1)
  1774.       if(NBSOU1.gt.0) then
  1775.       do isou1=1,NBSOU1
  1776.         MELEME=IPT1.LISOUS(isou1)
  1777.         segact,MELEME
  1778.         if(ITYPEL.eq.4) then
  1779.           IPT1ok=MELEME
  1780.           segdes,MELEME
  1781.           segsup,IPT1
  1782.           goto 9010
  1783.         else
  1784.           segsup,MELEME
  1785.         endif
  1786.         segdes,MELEME
  1787.       enddo
  1788.       endif
  1789.  9010 CONTINUE
  1790.  
  1791.       IF(SELF) GOTO 9020
  1792. C     ON REGENERE et on ne garde que les triangles (ITYPEL=4)
  1793.       CALL REGENE(IPT2)
  1794.       IPT2ok = IPT2
  1795.       segact,IPT2
  1796.       NBSOU1=IPT2.LISOUS(/1)
  1797.       if(NBSOU1.gt.0) then
  1798.       do isou1=1,NBSOU1
  1799.         MELEME=IPT2.LISOUS(isou1)
  1800.         segact,MELEME
  1801.         if(ITYPEL.eq.4) then
  1802.           IPT2ok=MELEME
  1803.           segdes,MELEME
  1804.           segsup,IPT2
  1805.           goto 9020
  1806.         else
  1807.           segsup,MELEME
  1808.         endif
  1809.         segdes,MELEME
  1810.       enddo
  1811.       endif
  1812.  9020 CONTINUE
  1813.  
  1814.  
  1815. c     on renvoie
  1816.       IOB1 = IPT1ok
  1817.       IF(.not.SELF) IOB2 = IPT2ok
  1818.       IOB3 = IPT3
  1819.       IF(.not.SELF) IOB4 = IPT4
  1820.  
  1821.  
  1822.       END
  1823.  
  1824.  
  1825.  
  1826.  
  1827.  
  1828.  
  1829.  

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