Télécharger la plaquette Cast3M


fr en

Télécharger crcoup.eso

Retour à la liste

Numérotation des lignes : 

crcoup

  1. C CRCOUP    SOURCE    PV        20/03/24    21:16:27     10554          
  2. C  CALCUL D'UN OBJET COUPE VIRTUEL
  3. C
  4. C  CRCOUP REMPLACE MELEME PAR L'OBJET COUPE
  5. C  MCOUP ET LE DERNIER COMPOSANT DE MELEME REMPLACE DECRIVENT
  6. C  L'INTERSECTION
  7. C  IL FAUDRA METTRE UNE INDICATION AUX TRIANGLES SE TROUVANT DANS LE PLA
  8. C  DE COUPE POUR REUSSIR LES PARTIES VUES PARTIES CACHEES
  9. C
  10. C SG 2016/07/18 : traitement du cas des faces tri7/qua9
  11. C
  12.       SUBROUTINE CRCOUP(JOEIL,ICOUP1,ICOUP2,ICOUP3,MELEME,MCOUP,VCHCA,
  13.      #  MELEM2,MCOU2,mchamo,isect)
  14.       IMPLICIT INTEGER(I-N)
  15. -INC CCGEOME
  16. -INC CCREEL
  17. -INC PPARAM
  18. -INC CCOPTIO
  19. -INC SMCOORD
  20. -INC SMELEME
  21. -INC SMCHAML
  22.       REAL*8 XN,YN,ZN,SIGA,SIGB,SIGC,SIGRAV
  23. * sg
  24.       REAL*8 xo,yo,zo,xct,yct,zct,xc2,yc2,zc2,xc3,yc3,zc3
  25.       REAL*8 xv,yv,zv,xw,yw,zw,xp,yp,zp,xgrav,ygrav,zgrav
  26.       REAL*8 xa,ya,za,xb,yb,zb,xc,yc,zc
  27.       REAL*8 rn,sig,chgrav
  28. * sg
  29.       SEGMENT VCHCA(0)
  30.       SEGMENT ITR(0)
  31.       SEGMENT ITRT(0)
  32.       SEGMENT XTR(0)
  33.       SEGMENT XTRT(0)
  34.       SEGMENT IQUAT(0)
  35.       SEGMENT XQUAT(0)
  36.       SEGMENT ISEGM(0)
  37.       SEGMENT XSEGM(NBPOIN)
  38.       SEGMENT MCOUP(0)
  39.       SEGMENT IJGRAV(2,100)
  40.       DIMENSION JTT(9),ISEG(9)
  41.       LOGICAL VOLUM
  42.       logical lquaf
  43. *
  44. *dbg      write(ioimp,*) 'coucou crcoup isect=',isect
  45.       n2ptel=0
  46.       mcham=mchamo
  47.       MEDSAU=0
  48.       SEGINI ITR,ISEGM,MCOUP,IJGRAV,ITRT,IQUAT
  49.       SEGACT MCOORD*mod
  50.       NBPOIN=nbpts
  51.       if (mcham.ne.0) segini xtr,xtrt,xquat,xsegm
  52.       IREF=(JOEIL-1)*4
  53.       XO=XCOOR(IREF+1)
  54.       YO=XCOOR(IREF+2)
  55.       ZO=XCOOR(IREF+3)
  56.       IREF=(ICOUP1-1)*4
  57.       XCT=XCOOR(IREF+1)
  58.       YCT=XCOOR(IREF+2)
  59.       ZCT=XCOOR(IREF+3)
  60.       IREF=(ICOUP2-1)*4
  61.       XC2=XCOOR(IREF+1)
  62.       YC2=XCOOR(IREF+2)
  63.       ZC2=XCOOR(IREF+3)
  64.       IREF=(ICOUP3-1)*4
  65.       XC3=XCOOR(IREF+1)
  66.       YC3=XCOOR(IREF+2)
  67.       ZC3=XCOOR(IREF+3)
  68. C  NORMALE AU PLAN :
  69.       XV=XC2-XCT
  70.       YV=YC2-YCT
  71.       ZV=ZC2-ZCT
  72.       XW=XC3-XCT
  73.       YW=YC3-YCT
  74.       ZW=ZC3-ZCT
  75.       XN=YV*ZW-ZV*YW
  76.       YN=ZV*XW-XV*ZW
  77.       ZN=XV*YW-YV*XW
  78.       RN=XN**2+YN**2+ZN**2
  79. *      IF (RN.EQ.0.) CALL ERREUR(21)
  80.       IF (ABS(RN).LT.(XPETIT**0.25D0)) CALL ERREUR(21)
  81.       IF (IERR.NE.0) RETURN
  82.       RN=SQRT(RN)
  83.       XN=XN/RN
  84.       YN=YN/RN
  85.       ZN=ZN/RN
  86.       SIG=(XO-XCT)*XN+(YO-YCT)*YN+(ZO-ZCT)*ZN
  87.       IF (SIG.LE.0.) THEN
  88.          XN=-XN
  89.          YN=-YN
  90.          ZN=-ZN
  91.       ENDIF
  92.  5000 CONTINUE
  93.       SEGACT MELEME
  94.       NBSOU=LISOUS(/1)
  95.       NBSOUS=MAX(2,NBSOU+1)
  96.       NBREF=0
  97.       if (mcham.ne.0) nbref=nbsous
  98.       NBNN=0
  99.       NBELEM=0
  100.       SEGINI IPT8
  101.       IPT1=MELEME
  102.       ISU=0
  103. *  prob optimiseur
  104. C       melva1=meleme
  105. C       melva2=meleme
  106.       DO 10 ISOUS=1,MAX(1,NBSOU)
  107.          IF (NBSOU.NE.0) THEN
  108.             IPT1=LISOUS(ISOUS)
  109.             if (mcham.ne.0) then
  110.                melva1=lisref(isous)
  111.                segact melva1
  112.                lva1=melva1.velche(/1)
  113.             endif
  114.             SEGACT IPT1
  115.          ENDIF
  116.          NBREF=0
  117.          NBSOUS=0
  118.          NBNN=IPT1.NUM(/1)
  119.          NBELEM=IPT1.NUM(/2)
  120.          SEGINI IPT2
  121.          IPT2.ITYPEL=IPT1.ITYPEL
  122.          if (mcham.ne.0) then
  123.             n1ptel=nbnn
  124.             n1el=nbelem
  125.             n2ptel=0
  126.             n2el=0
  127.             segini melva2
  128.             segact melva1
  129.             idi2= melva1.velche(/2)
  130.          endif
  131.          JELN=0
  132.          DO 20 JEL=1,NBELEM
  133.             DO 25 INOEU=1,NBNN
  134.                IPT=IPT1.NUM(INOEU,JEL)
  135.                IREF=4*(IPT-1)
  136.                XP=XCOOR(IREF+1)
  137.                YP=XCOOR(IREF+2)
  138.                ZP=XCOOR(IREF+3)
  139.                SIG=(XP-XCT)*XN+(YP-YCT)*YN+(ZP-ZCT)*ZN
  140.                IF (SIG.GT.0.) GOTO 30
  141.  25         CONTINUE
  142. *  tous les points sont du bon cote du plan de coupe
  143.             if (isect.ne.0) goto 20
  144.             JELN=JELN+1
  145.             DO 27 INOEU=1,NBNN
  146.                IPT2.NUM(INOEU,JELN)=IPT1.NUM(INOEU,JEL)
  147.                if (mcham.ne.0) melva2.velche(INOEU,JELN)=
  148.      $              melva1.velche(min(INOEU,lva1),JEL)
  149.  27         CONTINUE
  150.             IPT2.ICOLOR(JELN)=IPT1.ICOLOR(JEL)
  151.             GOTO 20
  152.  30         CONTINUE
  153.             VOLUM=KSURF(IPT1.ITYPEL).NE.IPT1.ITYPEL
  154. *     IF (VOLUM) THEN
  155.             IPT=IPT1.NUM(1,JEL)
  156.             IREF=(IPT-1)*4
  157.             XGRAV=XCOOR(IREF+1)
  158.             YGRAV=XCOOR(IREF+2)
  159.             ZGRAV=XCOOR(IREF+3)
  160.             IF (VCHCA.NE.0) VGRAV=VCHCA(IPT)
  161.             if (mcham.ne.0) chgrav=melva1.velche(1,min(JEL,idi2))
  162.             DO 21 INOEU=2,NBNN
  163.                IPT=IPT1.NUM(INOEU,JEL)
  164.                IREF=(IPT-1)*4
  165.                XGRAV=XGRAV+XCOOR(IREF+1)
  166.                YGRAV=YGRAV+XCOOR(IREF+2)
  167.                ZGRAV=ZGRAV+XCOOR(IREF+3)
  168.                IF (VCHCA.NE.0) VGRAV=VGRAV+VCHCA(IPT)
  169.                if (mcham.ne.0) chgrav=chgrav+
  170.      $              melva1.velche(min(INOEU,lva1),min(idi2,JEL))
  171.  21         CONTINUE
  172.             XGRAV=XGRAV/NBNN
  173.             YGRAV=YGRAV/NBNN
  174.             ZGRAV=ZGRAV/NBNN
  175. C  ON PREND COMME POINT DU PLAN LA PROJECTION DU CENTRE GRAVITE
  176.             SIGRAV=(XGRAV-XCT)*XN+(YGRAV-YCT)*YN+(ZGRAV-ZCT)*ZN
  177.             XCT=XGRAV-SIGRAV*XN
  178.             YCT=YGRAV-SIGRAV*YN
  179.             ZCT=ZGRAV-SIGRAV*ZN
  180.             SIGRAV=(XGRAV-XCT)*XN+(YGRAV-YCT)*YN+(ZGRAV-ZCT)*ZN
  181.             IF (ABS(SIGRAV).LE.1D-30) SIGRAV=1D-30
  182.             NJGRAV=0
  183.             IF (VCHCA.NE.0) VGRAV=VGRAV/NBNN
  184.             if (mcham.ne.0) CHGRAV=CHGRAV/NBNN
  185. *     ENDIF
  186. C  DECOMPOSITION EN FACE PUIS EN TRIANGLE
  187.             NBFAC=LTEL(1,IPT1.ITYPEL)
  188.             IAD=LTEL(2,IPT1.ITYPEL)-1
  189.             IF (NBFAC.EQ.0) GOTO 20
  190.             DO 161 IFAC=1,NBFAC
  191.                ITYP=LDEL(1,IAD+IFAC)
  192. * SG  Pour les TRI7/QUA9, il y a un traitement particulier pour ne pas
  193. *     tracer les segments relies au noeud central. Il prend en compte le
  194. *     fait que le decoupage est barycentrique
  195.          lquaf=(ityp.eq.7.OR.ityp.eq.8) 
  196.          JAD=LDEL(2,IAD+IFAC)-1
  197.          NPFAC=KDFAC(1,ITYP)
  198.          IDEP=KDFAC(2,ITYP)
  199.          IFEP=IDEP+3*(KDFAC(3,ITYP)-1)
  200. *         write(ioimp,*) 'ifac,npfac=',ifac,npfac
  201.          DO 160 ITRIAN=IDEP,IFEP,3
  202.             kfa=kfac(itrian)
  203.             kfb=kfac(itrian+1)
  204.             kfc=kfac(itrian+2)
  205.             JFA=JAD+KFAC(ITRIAN)
  206.             JFB=JAD+KFAC(ITRIAN+1)
  207.             JFC=JAD+KFAC(ITRIAN+2)
  208.             ISEG(1)=0
  209.             ISEG(2)=0
  210.             ISEG(3)=0
  211.             ISEG(4)=0
  212.             ISEG1=0
  213.             ISEG2=0
  214.             ISEG3=0
  215.             if (lquaf) then
  216.                if (KFA.NE.NPFAC.AND.KFB.NE.NPFAC) ISEG1=1
  217.                if (KFB.NE.NPFAC.AND.KFC.NE.NPFAC) ISEG2=1
  218.                if (KFA.NE.NPFAC.AND.KFC.NE.NPFAC) ISEG3=1
  219.             else
  220. *               imod1 = MOD(JFB-JFA+NPFAC,NPFAC)
  221. *               imod2 = MOD(JFC-JFB+NPFAC,NPFAC)
  222. *               imod3 = MOD(JFA-JFC+NPFAC,NPFAC)
  223.                IF (MOD(JFB-JFA+NPFAC,NPFAC).EQ.1) ISEG1=1
  224.                IF (MOD(JFB-JFA+NPFAC,NPFAC).EQ.NPFAC-1) ISEG1=1
  225.                IF (MOD(JFC-JFB+NPFAC,NPFAC).EQ.1) ISEG2=1
  226.                IF (MOD(JFC-JFB+NPFAC,NPFAC).EQ.NPFAC-1) ISEG2=1
  227.                IF (MOD(JFA-JFC+NPFAC,NPFAC).EQ.1) ISEG3=1
  228.                IF (MOD(JFA-JFC+NPFAC,NPFAC).EQ.NPFAC-1) ISEG3=1
  229. *               write(ioimp,*) '   itrian,kfa,kfb,kfc=',itrian
  230. *     $              ,kfac(itrian),kfac(itrian+1),kfac(itrian+2)
  231. *               write(ioimp,*) '   iseg1,iseg2,iseg3=',iseg1,iseg2,iseg3
  232. *               write(ioimp,*) '   imod1,imod2,imod3=',imod1,imod2,imod3
  233. *               write(ioimp,*) '   imod1p,imod2p,imod3p=',imod1p,imod2p
  234. *     $              ,imod3p
  235.             endif
  236.             IAFA=LFAC(JFA)
  237.             IBFA=LFAC(JFB)
  238.             ICFA=LFAC(JFC)
  239.             IA=IPT1.NUM(IAFA,JEL)
  240.             IB=IPT1.NUM(IBFA,JEL)
  241.             IC=IPT1.NUM(ICFA,JEL)
  242.             IREF=4*(IA-1)
  243.             XA=XCOOR(IREF+1)
  244.             YA=XCOOR(IREF+2)
  245.             ZA=XCOOR(IREF+3)
  246.             IREF=4*(IB-1)
  247.             XB=XCOOR(IREF+1)
  248.             YB=XCOOR(IREF+2)
  249.             ZB=XCOOR(IREF+3)
  250.             IREF=4*(IC-1)
  251.             XC=XCOOR(IREF+1)
  252.             YC=XCOOR(IREF+2)
  253.             ZC=XCOOR(IREF+3)
  254. C  TESTONS SI LES POINTS SONT DU BON COTE
  255.             SIGA=(XA-XCT)*XN+(YA-YCT)*YN+(ZA-ZCT)*ZN
  256.             IF (ABS(SIGA).LE.1D-30) SIGA=1D-30
  257.             SIGB=(XB-XCT)*XN+(YB-YCT)*YN+(ZB-ZCT)*ZN
  258.             IF (ABS(SIGB).LE.1D-30) SIGB=1D-30
  259.             SIGC=(XC-XCT)*XN+(YC-YCT)*YN+(ZC-ZCT)*ZN
  260.             IF (ABS(SIGC).LE.1D-30) SIGC=1D-30
  261. *     if (isect.eq.1) goto 350
  262.             IF (SIGA.GT.0..OR.SIGB.GT.0..OR.SIGC.GT.0.) GOTO 103
  263.             IF (VOLUM) THEN
  264.                IA1=MIN(IA,IB,IC)
  265.                IB1=MAX(IA,IB,IC)
  266.                IC1=IA+IB+IC
  267.                DO 102 IDJF=ITRT(/1)-4,1,-5
  268.                   ITR1=MIN(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
  269.                   ITR2=MAX(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
  270.                   ITR3=ITRT(IDJF)+ITRT(IDJF+1)+ITRT(IDJF+2)
  271.                   IF (IA1.NE.ITR1) GOTO 102
  272.                   IF (IB1.NE.ITR2) GOTO 102
  273.                   IF (IC1.NE.ITR3) GOTO 102
  274. *      WRITE (6,*) ' TRIANGLE EN DOUBLE CAS 1 ',IA,IB,IC
  275. C  LE TRIANGLE EST DECLARE NON TRACABLE
  276.                   ITRT(IDJF+4)=MOD(ITRT(IDJF+4),16)+16
  277.                   GOTO 301
  278.  102           CONTINUE
  279.             ENDIF
  280.             if (isect.eq.0) then
  281. *       WRITE (6,*) ' NOUVEAU TRIANGLE ',IA,IB,IC
  282.                ITRT(**)=IA
  283.                ITRT(**)=IB
  284.                ITRT(**)=IC
  285.                ITRT(**)=IPT1.ICOLOR(JEL)
  286.                ITRT(**)=ISEG1+2*ISEG2+4*ISEG3
  287.                if (mcham.ne.0) then
  288.                   xtrt(**)=melva1.velche(min(lva1,iafa),min(idi2,jel))
  289.                   xtrt(**)=melva1.velche(min(lva1,ibfa),min(idi2,jel))
  290.                   xtrt(**)=melva1.velche(min(lva1,icfa),min(idi2,jel))
  291.                endif
  292.             endif
  293.  301        CONTINUE
  294.             GOTO 101
  295.  103        CONTINUE
  296.             IF (SIGA.GT.0..AND.SIGB.GT.0..AND.SIGC.GT.0.) GOTO 101
  297.             IPOS=0
  298.             IF (SIGA.LE.0.) THEN
  299.                IPOS=IPOS+1
  300.                JTT(IPOS)=IA
  301.                ISEG(IPOS)=ISEG1
  302.                if(mcham.ne.0) xsegm(ia)=
  303.      $              melva1.velche(min(lva1,iafa),min(idi2,jel))
  304.             ENDIF
  305.             IF (SIGA*SIGB.LT.0..AND.(.NOT.(SIGA.LE.0..AND.SIGB.LE.0.)))
  306.      $           THEN
  307.                IA1=MIN(IA,IB)
  308.                IB1=MAX(IA,IB)
  309.                DO 104 NSEGM=ISEGM(/1)-2,1,-3
  310.                   IF (IA1.NE.ISEGM(NSEGM)) GOTO 104
  311.                   IF (IB1.NE.ISEGM(NSEGM+1)) GOTO 104
  312.                   ISUP=ISEGM(NSEGM+2)
  313.                   if (mcham.ne.0) xsegm(isup)=
  314.      $               (melva1.velche(min(lva1,iafa),min(idi2,jel))*sigb
  315.      $               -melva1.velche(min(lva1,ibfa),min(idi2,jel))*siga)
  316.      $                 /(sigb-siga)
  317. *       WRITE (6,*) ' POINT EN DOUBLE ',IA1,IB1,ISUP
  318.                   GOTO 105
  319.  104           CONTINUE
  320.                NBPTS=nbpts+1
  321.                SEGADJ MCOORD
  322.                XCOOR((NBPTS-1)*4+1)=(XA*SIGB-XB*SIGA)/(SIGB-SIGA)
  323.                XCOOR((NBPTS-1)*4+2)=(YA*SIGB-YB*SIGA)/(SIGB-SIGA)
  324.                XCOOR((NBPTS-1)*4+3)=(ZA*SIGB-ZB*SIGA)/(SIGB-SIGA)
  325.                XCOOR(NBPTS*(IDIM+1))=0.D0
  326.                IF (VCHCA.NE.0)
  327.      $              VCHCA(**)=(VCHCA(IA)*SIGB-VCHCA(IB)*SIGA)
  328.      $              /(SIGB-SIGA)
  329.                ISUP=NBPTS
  330.                ISEGM(**)=IA1
  331.                ISEGM(**)=IB1
  332.                ISEGM(**)=ISUP
  333.                if (mcham.ne.0)
  334.      >    xsegm(**)=(melva1.velche(min(lva1,iafa),min(idi2,jel))*sigb-
  335.      >    melva1.velche(min(lva1,ibfa),min(idi2,jel))*siga)/(sigb-siga)
  336.  105           CONTINUE
  337.                IPOS=IPOS+1
  338.                JTT(IPOS)=ISUP
  339.                IF (SIGB.LE.0.) THEN
  340.                   ISEG(IPOS)=ISEG1
  341.                ELSE
  342.                   ISEG(IPOS)=1
  343.                ENDIF
  344.             ENDIF
  345.             IF (SIGB.LE.0.) THEN
  346.                IPOS=IPOS+1
  347.                JTT(IPOS)=IB
  348.                ISEG(IPOS)=ISEG2
  349.                if (mcham.ne.0) 
  350.      $             xsegm(ib)=melva1.velche(min(lva1,ibfa),min(idi2,jel))
  351.             ENDIF
  352.             IF (SIGB*SIGC.LT.0..AND.(.NOT.(SIGB.LE.0..AND.SIGC.LE.0.)))
  353.      $           THEN
  354.                IB1=MIN(IB,IC)
  355.                IC1=MAX(IB,IC)
  356.                DO 106 NSEGM=ISEGM(/1)-2,1,-3
  357.                   IF (IB1.NE.ISEGM(NSEGM)) GOTO 106
  358.                   IF (IC1.NE.ISEGM(NSEGM+1)) GOTO 106
  359.                   ISUP=ISEGM(NSEGM+2)
  360.                   if (mcham.ne.0)
  361.      $                 xsegm(isup)=(
  362.      $                 melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigc
  363.      $                -melva1.velche(min(lva1,icfa),min(idi2,jel))*sigb)
  364.      $                 /(sigc-sigb)
  365. *       WRITE (6,*) ' POINT EN DOUBLE ',IB1,IC1,ISUP
  366.                   GOTO 107
  367.  106           CONTINUE
  368.                NBPTS=nbpts+1
  369.                SEGADJ MCOORD
  370.                XCOOR((NBPTS-1)*4+1)=(XB*SIGC-XC*SIGB)/(SIGC-SIGB)
  371.                XCOOR((NBPTS-1)*4+2)=(YB*SIGC-YC*SIGB)/(SIGC-SIGB)
  372.                XCOOR((NBPTS-1)*4+3)=(ZB*SIGC-ZC*SIGB)/(SIGC-SIGB)
  373.                XCOOR(NBPTS*4)=0.D0
  374.                IF (VCHCA.NE.0)
  375.      $              VCHCA(**)=(VCHCA(IB)*SIGC-VCHCA(IC)*SIGB)/(SIGC
  376.      $              -SIGB)
  377.                ISUP=NBPTS
  378.                ISEGM(**)=IB1
  379.                ISEGM(**)=IC1
  380.                ISEGM(**)=ISUP
  381.                if (mcham.ne.0)
  382.      $              xsegm(**)=(
  383.      $              melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigc
  384.      $             -melva1.velche(min(lva1,icfa),min(idi2,jel))*sigb)
  385.      $              /(sigc-sigb)
  386.  107           CONTINUE
  387.                IPOS=IPOS+1
  388.                JTT(IPOS)=ISUP
  389.                IF (SIGC.LE.0.) THEN
  390.                   ISEG(IPOS)=ISEG2
  391.                ELSE
  392.                   ISEG(IPOS)=1
  393.                ENDIF
  394.             ENDIF
  395.             IF (SIGC.LE.0.) THEN
  396.                IPOS=IPOS+1
  397.                JTT(IPOS)=IC
  398.                ISEG(IPOS)=ISEG3
  399.                if (mcham.ne.0)
  400.      $             xsegm(ic)=melva1.velche(min(lva1,icfa),min(idi2,jel))
  401.             ENDIF
  402.             IF (SIGC*SIGA.LT.0..AND.(.NOT.(SIGC.LE.0..AND.SIGA.LE.0.)))
  403.      $           THEN
  404.                IC1=MIN(IC,IA)
  405.                IA1=MAX(IC,IA)
  406.                DO 108 NSEGM=ISEGM(/1)-2,1,-3
  407.                   IF (IC1.NE.ISEGM(NSEGM)) GOTO 108
  408.                   IF (IA1.NE.ISEGM(NSEGM+1)) GOTO 108
  409.                   ISUP=ISEGM(NSEGM+2)
  410.                   if (mcham.ne.0)
  411.      $                 xsegm(isup)=(
  412.      $                 melva1.velche(min(lva1,icfa),min(idi2,jel))*siga
  413.      $                -melva1.velche(min(lva1,iafa),min(idi2,jel))*sigc)
  414.      $                 /(siga-sigc)
  415. *       WRITE (6,*) ' POINT EN DOUBLE ',IC1,IA1,ISUP
  416.                   GOTO 109
  417.  108           CONTINUE
  418.                NBPTS=nbpts+1
  419.                SEGADJ MCOORD
  420.                XCOOR((NBPTS-1)*4+1)=(XC*SIGA-XA*SIGC)/(SIGA-SIGC)
  421.                XCOOR((NBPTS-1)*4+2)=(YC*SIGA-YA*SIGC)/(SIGA-SIGC)
  422.                XCOOR((NBPTS-1)*4+3)=(ZC*SIGA-ZA*SIGC)/(SIGA-SIGC)
  423.                XCOOR(NBPTS*4)=0.D0
  424.                IF (VCHCA.NE.0)
  425.      $              VCHCA(**)=(VCHCA(IC)*SIGA-VCHCA(IA)*SIGC)/(SIGA
  426.      $              -SIGC)
  427.                ISUP=NBPTS
  428.                ISEGM(**)=IC1
  429.                ISEGM(**)=IA1
  430.                ISEGM(**)=ISUP
  431.                if (mcham.ne.0)
  432.      $              xsegm(**)=(
  433.      $              melva1.velche(min(lva1,icfa),min(idi2,jel))*siga
  434.      $             -melva1.velche(min(lva1,iafa),min(idi2,jel))*sigc)
  435.      $              /(siga-sigc)
  436.  109           CONTINUE
  437.                IPOS=IPOS+1
  438.                JTT(IPOS)=ISUP
  439.                IF (SIGA.LE.0.) THEN
  440.                   ISEG(IPOS)=ISEG3
  441.                ELSE
  442.                   ISEG(IPOS)=1
  443.                ENDIF
  444.             ENDIF
  445.             IF (IPOS.GT.4) WRITE (6,*) ' ANOMALIE IPOS SURFACE ',IPOS
  446.             IF (IPOS.LT.3) WRITE (6,*) ' ANOMALIE IPOS SURFACE ',IPOS
  447.             IF (IPOS.LT.3) GOTO 101
  448.             IF (IPOS.EQ.3) THEN
  449.                IF (VOLUM) THEN
  450.                   IA1=MIN(JTT(1),JTT(2),JTT(3))
  451.                   IB1=MAX(JTT(1),JTT(2),JTT(3))
  452.                   IC1=JTT(1)+JTT(2)+JTT(3)
  453.                   DO 110 IDJF=ITRT(/1)-4,1,-5
  454.                      ITR1=MIN(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
  455.                      ITR2=MAX(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
  456.                      ITR3=ITRT(IDJF)+ITRT(IDJF+1)+ITRT(IDJF+2)
  457.                      IF (IA1.NE.ITR1) GOTO 110
  458.                      IF (IB1.NE.ITR2) GOTO 110
  459.                      IF (IC1.NE.ITR3) GOTO 110
  460. *      WRITE (6,*) ' TRIANGLE EN DOUBLE CAS 2 ',JTT(1),JTT(2),JTT(3)
  461. C  LE TRIANGLE EST DECLARE NON TRACABLE
  462.                      ITRT(IDJF+4)=MOD(ITRT(IDJF+4),16)+16
  463.                      GOTO 303
  464.  110              CONTINUE
  465.                ENDIF
  466. *     WRITE (6,*) ' NOUVEAU TRIANGLE ',JTT(1),JTT(2),JTT(3)
  467.                if (isect.eq.1.and.jtt(1).le.nbpoin) jtt(1)=isup
  468.                ITRT(**)=JTT(1)
  469.                if (isect.eq.1.and.jtt(2).le.nbpoin) jtt(2)=isup
  470.                ITRT(**)=JTT(2)
  471.                if (isect.eq.1.and.jtt(3).le.nbpoin) jtt(3)=isup
  472.                ITRT(**)=JTT(3)
  473.                ITRT(**)=IPT1.ICOLOR(JEL)
  474.                ITRT(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)
  475.                if (mcham.ne.0) then
  476.                   xtrt(**)=xsegm(jtt(1))
  477.                   xtrt(**)=xsegm(jtt(2))
  478.                   xtrt(**)=xsegm(jtt(3))
  479.                endif
  480.             ENDIF
  481.  303        CONTINUE
  482.             IF (IPOS.EQ.4) THEN
  483.                IF (VOLUM) THEN
  484.       IA1=MIN(JTT(1),JTT(2),JTT(3),JTT(4))
  485.       IB1=MAX(JTT(1),JTT(2),JTT(3),JTT(4))
  486.       IC1=JTT(1)+JTT(2)+JTT(3)+JTT(4)
  487.       ID1=JTT(1)*JTT(2)*JTT(3)*JTT(4)
  488.       DO 111 IDJF=IQUAT(/1)-5,1,-6
  489.       ITR1=MIN(IQUAT(IDJF),IQUAT(IDJF+1),IQUAT(IDJF+2),IQUAT(IDJF+3))
  490.       ITR2=MAX(IQUAT(IDJF),IQUAT(IDJF+1),IQUAT(IDJF+2),IQUAT(IDJF+3))
  491.       ITR3=IQUAT(IDJF)+IQUAT(IDJF+1)+IQUAT(IDJF+2)+IQUAT(IDJF+3)
  492.       ITR4=IQUAT(IDJF)*IQUAT(IDJF+1)*IQUAT(IDJF+2)*IQUAT(IDJF+3)
  493.        IF (IA1.NE.ITR1) GOTO 111
  494.        IF (IB1.NE.ITR2) GOTO 111
  495.        IF (IC1.NE.ITR3) GOTO 111
  496.        IF (ID1.NE.ITR4) GOTO 111
  497. *      WRITE (6,*) ' QUADRILA EN DOUBLE ',JTT(1),JTT(2),JTT(3),JTT(4)
  498. C  LE QUADRILATERE EST DECLARE NON TRACABLE
  499.        IQUAT(IDJF+5)=MOD(IQUAT(IDJF+5),16)+16
  500.        GOTO 302
  501.  111  CONTINUE
  502.       ENDIF
  503. *     WRITE (6,*) ' NOUVEAU QUADRILATERE ',JTT(1),JTT(2),JTT(3),JTT(4)
  504.        if (isect.eq.1.and.jtt(1).le.nbpoin) jtt(1)=isup
  505.       IQUAT(**)=JTT(1)
  506.        if (isect.eq.1.and.jtt(2).le.nbpoin) jtt(2)=isup
  507.       IQUAT(**)=JTT(2)
  508.        if (isect.eq.1.and.jtt(3).le.nbpoin) jtt(3)=isup
  509.       IQUAT(**)=JTT(3)
  510.        if (isect.eq.1.and.jtt(4).le.nbpoin) jtt(4)=isup
  511.       IQUAT(**)=JTT(4)
  512.       IQUAT(**)=IPT1.ICOLOR(JEL)
  513.       IQUAT(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)+8*ISEG(4)
  514.       if (mcham.ne.0) then
  515.        xquat(**)=xsegm(jtt(1))
  516.        xquat(**)=xsegm(jtt(2))
  517.        xquat(**)=xsegm(jtt(3))
  518.        xquat(**)=xsegm(jtt(4))
  519.       endif
  520.       ENDIF
  521.  302  CONTINUE
  522.       GOTO 101
  523.  101  CONTINUE
  524.  350  CONTINUE
  525.       IF (.NOT.VOLUM) GOTO 160
  526. C  ETUDE POUR LE CAS DES VOLUMES DU TETRAEDRE S'APPUYANT SUR LE
  527. C  TRIANGLE ET LE CENTRE DE GRAVITE
  528. C  TETRAEDRE
  529.       IF (SIGA.LT.0..AND.SIGB.LT.0..AND.SIGC.LT.0..AND.SIGRAV
  530.      #  .LT.0.) GOTO 160
  531.       IF (SIGA.GT.0..AND.SIGB.GT.0..AND.SIGC.GT.0..AND.SIGRAV
  532.      #  .GT.0.) GOTO 160
  533.       IPOS=0
  534.       ISEG(1)=0
  535.       ISEG(2)=0
  536.       ISEG(3)=0
  537.       ISEG(4)=0
  538.       IF (SIGRAV*SIGA.LE.0.) THEN
  539.        IF (SIGA.EQ.0.) THEN
  540.         ISUP=IA
  541.         if (mcham.ne.0)
  542.      >    xsegm(isup)=melva1.velche(min(lva1,iafa),min(idi2,jel))
  543.         GOTO 359
  544.        ENDIF
  545.         DO 358 NSEGM=1,NJGRAV
  546.         IF (IA.NE.IJGRAV(1,NSEGM)) GOTO 358
  547.         ISUP=IJGRAV(2,NSEGM)
  548.         if (mcham.ne.0)
  549.      >  xsegm(isup)=(melva1.velche(min(lva1,iafa),min(idi2,jel))*sigrav-
  550.      >              chgrav*siga)/(sigrav-siga)
  551.         GOTO 359
  552.  358    CONTINUE
  553.         NBPTS=nbpts+1
  554.         SEGADJ MCOORD
  555.         XCOOR((NBPTS-1)*4+1)=(XA*SIGRAV-XGRAV*SIGA)/(SIGRAV-SIGA)
  556.         XCOOR((NBPTS-1)*4+2)=(YA*SIGRAV-YGRAV*SIGA)/(SIGRAV-SIGA)
  557.         XCOOR((NBPTS-1)*4+3)=(ZA*SIGRAV-ZGRAV*SIGA)/(SIGRAV-SIGA)
  558.         XCOOR(NBPTS*4)=0.D0
  559.         IF (VCHCA.NE.0)
  560.      #   VCHCA(**)=(VCHCA(IA)*SIGRAV-VGRAV*SIGA)/(SIGRAV-SIGA)
  561.        if (mcham.ne.0)
  562.      >   xsegm(**)=(melva1.velche(min(lva1,iafa),min(idi2,jel))*sigrav-
  563.      >              chgrav*siga)/(sigrav-siga)
  564.         ISUP=NBPTS
  565.         NJGRAV=NJGRAV+1
  566.         IJGRAV(1,NJGRAV)=IA
  567.         IJGRAV(2,NJGRAV)=ISUP
  568.  359    CONTINUE
  569.         IPOS=IPOS+1
  570.         JTT(IPOS)=ISUP
  571.       ENDIF
  572.       IF (SIGA*SIGB.LT.0..AND.
  573.      # (.NOT.(SIGRAV*SIGA.LE.0..AND.SIGRAV*SIGB.LE.0.))) THEN
  574.         IA1=MIN(IA,IB)
  575.         IB1=MAX(IA,IB)
  576.         DO 404 NSEGM=ISEGM(/1)-2,1,-3
  577.         IF (IA1.NE.ISEGM(NSEGM)) GOTO 404
  578.         IF (IB1.NE.ISEGM(NSEGM+1)) GOTO 404
  579.         ISUP=ISEGM(NSEGM+2)
  580. *       WRITE (6,*) ' POINT EN DOUBLE ',IA1,IB1,ISUP
  581.         if (mcham.ne.0)
  582.      >  xsegm(isup)=(melva1.velche(min(lva1,iafa),min(idi2,jel))*sigb-
  583.      >   melva1.velche(min(lva1,ibfa),min(idi2,jel))*siga)/(sigb-siga)
  584.         GOTO 405
  585.  404    CONTINUE
  586.         NBPTS=nbpts+1
  587.         SEGADJ MCOORD
  588.         XCOOR((NBPTS-1)*4+1)=(XA*SIGB-XB*SIGA)/(SIGB-SIGA)
  589.         XCOOR((NBPTS-1)*4+2)=(YA*SIGB-YB*SIGA)/(SIGB-SIGA)
  590.         XCOOR((NBPTS-1)*4+3)=(ZA*SIGB-ZB*SIGA)/(SIGB-SIGA)
  591.         XCOOR(NBPTS*4)=0.D0
  592.         IF (VCHCA.NE.0)
  593.      #   VCHCA(**)=(VCHCA(IA)*SIGB-VCHCA(IB)*SIGA)/(SIGB-SIGA)
  594.         ISUP=NBPTS
  595.         ISEGM(**)=IA1
  596.         ISEGM(**)=IB1
  597.         ISEGM(**)=ISUP
  598.         if (mcham.ne.0)
  599.      >    xsegm(**)=(melva1.velche(min(lva1,iafa),jel)*sigb-
  600.      >               melva1.velche(min(lva1,ibfa),jel)*siga)/(sigb-siga)
  601.  405    CONTINUE
  602.         IPOS=IPOS+1
  603.         JTT(IPOS)=ISUP
  604.       IF (SIGRAV*SIGB.GT.0.) THEN
  605.         ISEG(IPOS)=1
  606.       ENDIF
  607.       ENDIF
  608.       IF (SIGRAV*SIGB.LE.0.) THEN
  609.        IF (SIGB.EQ.0.) THEN
  610.         ISUP=IB
  611.         if (mcham.ne.0)
  612.      >    xsegm(isup)=melva1.velche(min(lva1,ibfa),min(idi2,jel))
  613.         GOTO 369
  614.        ENDIF
  615.         DO 368 NSEGM=1,NJGRAV
  616.         IF (IB.NE.IJGRAV(1,NSEGM)) GOTO 368
  617.         ISUP=IJGRAV(2,NSEGM)
  618.         GOTO 369
  619.  368    CONTINUE
  620.         NBPTS=nbpts+1
  621.         SEGADJ MCOORD
  622.         XCOOR((NBPTS-1)*4+1)=(XB*SIGRAV-XGRAV*SIGB)/(SIGRAV-SIGB)
  623.         XCOOR((NBPTS-1)*4+2)=(YB*SIGRAV-YGRAV*SIGB)/(SIGRAV-SIGB)
  624.         XCOOR((NBPTS-1)*4+3)=(ZB*SIGRAV-ZGRAV*SIGB)/(SIGRAV-SIGB)
  625.         XCOOR(NBPTS*4)=0.D0
  626.         IF (VCHCA.NE.0)
  627.      #   VCHCA(**)=(VCHCA(IB)*SIGRAV-VGRAV*SIGB)/(SIGRAV-SIGB)
  628.         if (mcham.ne.0)
  629.      >  xsegm(**)=(melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigrav-
  630.      >               chgrav*sigb)/(sigrav-sigb)
  631.         ISUP=NBPTS
  632.         NJGRAV=NJGRAV+1
  633.         IJGRAV(1,NJGRAV)=IB
  634.         IJGRAV(2,NJGRAV)=ISUP
  635.  369    CONTINUE
  636.         IPOS=IPOS+1
  637.         JTT(IPOS)=ISUP
  638.         if (mcham.ne.0)
  639.      >  xsegm(isup)=(melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigrav-
  640.      >               chgrav*sigb)/(sigrav-sigb)
  641.       ENDIF
  642.       IF (SIGB*SIGC.LT.0..AND.
  643.      # (.NOT.(SIGRAV*SIGB.LE.0..AND.SIGRAV*SIGC.LE.0.))) THEN
  644.         IB1=MIN(IB,IC)
  645.         IC1=MAX(IB,IC)
  646.         DO 406 NSEGM=ISEGM(/1)-2,1,-3
  647.         IF (IB1.NE.ISEGM(NSEGM)) GOTO 406
  648.         IF (IC1.NE.ISEGM(NSEGM+1)) GOTO 406
  649.         ISUP=ISEGM(NSEGM+2)
  650. *       WRITE (6,*) ' POINT EN DOUBLE ',IB1,IC1,ISUP
  651.         if (mcham.ne.0)
  652.      >  xsegm(isup)=(melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigc-
  653.      >  melva1.velche(min(lva1,icfa),min(idi2,jel))*sigb)/(sigc-sigb)
  654.         GOTO 407
  655.  406    CONTINUE
  656.         NBPTS=nbpts+1
  657.         SEGADJ MCOORD
  658.         XCOOR((NBPTS-1)*4+1)=(XB*SIGC-XC*SIGB)/(SIGC-SIGB)
  659.         XCOOR((NBPTS-1)*4+2)=(YB*SIGC-YC*SIGB)/(SIGC-SIGB)
  660.         XCOOR((NBPTS-1)*4+3)=(ZB*SIGC-ZC*SIGB)/(SIGC-SIGB)
  661.         XCOOR(NBPTS*4)=0.D0
  662.         IF (VCHCA.NE.0)
  663.      #   VCHCA(**)=(VCHCA(IB)*SIGC-VCHCA(IC)*SIGB)/(SIGC-SIGB)
  664.         ISUP=NBPTS
  665.         ISEGM(**)=IB1
  666.         ISEGM(**)=IC1
  667.         ISEGM(**)=ISUP
  668.         if (mcham.ne.0)
  669.      >  xsegm(**)=(melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigc-
  670.      >  melva1.velche(min(lva1,icfa),min(idi2,jel))*sigb)/(sigc-sigb)
  671.  407    CONTINUE
  672.         IPOS=IPOS+1
  673.         JTT(IPOS)=ISUP
  674.       IF (SIGRAV*SIGC.GT.0.) THEN
  675.         ISEG(IPOS)=1
  676.       ENDIF
  677.       ENDIF
  678.       IF (SIGRAV*SIGC.LE.0.) THEN
  679.        IF (SIGC.EQ.0.) THEN
  680.         ISUP=IC
  681.        if (mcham.ne.0)
  682.      >    xsegm(isup)=melva1.velche(min(lva1,icfa),min(idi2,jel))
  683.         GOTO 379
  684.        ENDIF
  685.         DO 378 NSEGM=1,NJGRAV
  686.         IF (IC.NE.IJGRAV(1,NSEGM)) GOTO 378
  687.         ISUP=IJGRAV(2,NSEGM)
  688.         if (mcham.ne.0)
  689.      >  xsegm(isup)=(melva1.velche(min(lva1,icfa),min(idi2,jel))*sigrav-
  690.      >               chgrav*sigc)/(sigrav-sigc)
  691.         GOTO 379
  692.  378    CONTINUE
  693.         NBPTS=nbpts+1
  694.         SEGADJ MCOORD
  695.         XCOOR((NBPTS-1)*4+1)=(XC*SIGRAV-XGRAV*SIGC)/(SIGRAV-SIGC)
  696.         XCOOR((NBPTS-1)*4+2)=(YC*SIGRAV-YGRAV*SIGC)/(SIGRAV-SIGC)
  697.         XCOOR((NBPTS-1)*4+3)=(ZC*SIGRAV-ZGRAV*SIGC)/(SIGRAV-SIGC)
  698.         XCOOR(NBPTS*4)=0.D0
  699.         IF (VCHCA.NE.0)
  700.      #   VCHCA(**)=(VCHCA(IC)*SIGRAV-VGRAV*SIGC)/(SIGRAV-SIGC)
  701.         if (mcham.ne.0)
  702.      >  xsegm(**)=(melva1.velche(min(lva1,icfa),min(idi2,jel))*sigrav-
  703.      >               chgrav*sigc)/(sigrav-sigc)
  704.         ISUP=NBPTS
  705.         NJGRAV=NJGRAV+1
  706.         IJGRAV(1,NJGRAV)=IC
  707.         IJGRAV(2,NJGRAV)=ISUP
  708.  379    CONTINUE
  709.         IPOS=IPOS+1
  710.         JTT(IPOS)=ISUP
  711.       ENDIF
  712.       IF (SIGC*SIGA.LT.0..AND.
  713.      # (.NOT.(SIGRAV*SIGC.LE.0..AND.SIGRAV*SIGA.LE.0.))) THEN
  714.         IC1=MIN(IC,IA)
  715.         IA1=MAX(IC,IA)
  716.         DO 408 NSEGM=ISEGM(/1)-2,1,-3
  717.         IF (IC1.NE.ISEGM(NSEGM)) GOTO 408
  718.         IF (IA1.NE.ISEGM(NSEGM+1)) GOTO 408
  719.         ISUP=ISEGM(NSEGM+2)
  720. *       WRITE (6,*) ' POINT EN DOUBLE ',IC1,IA1,ISUP
  721.        if (mcham.ne.0)
  722.      > xsegm(isup)=(melva1.velche(min(lva1,icfa),min(idi2,jel))*siga-
  723.      > melva1.velche(min(lva1,iafa),min(idi2,jel))*sigc)/(siga-sigc)
  724.         GOTO 409
  725.  408    CONTINUE
  726.         NBPTS=nbpts+1
  727.         SEGADJ MCOORD
  728.         XCOOR((NBPTS-1)*4+1)=(XC*SIGA-XA*SIGC)/(SIGA-SIGC)
  729.         XCOOR((NBPTS-1)*4+2)=(YC*SIGA-YA*SIGC)/(SIGA-SIGC)
  730.         XCOOR((NBPTS-1)*4+3)=(ZC*SIGA-ZA*SIGC)/(SIGA-SIGC)
  731.         XCOOR(NBPTS*4)=0.D0
  732.         IF (VCHCA.NE.0)
  733.      #   VCHCA(**)=(VCHCA(IC)*SIGA-VCHCA(IA)*SIGC)/(SIGA-SIGC)
  734.         ISUP=NBPTS
  735.         ISEGM(**)=IC1
  736.         ISEGM(**)=IA1
  737.         ISEGM(**)=ISUP
  738.         if (mcham.ne.0)
  739.      >  xsegm(**)=(melva1.velche(min(lva1,icfa),min(idi2,jel))*siga-
  740.      >   melva1.velche(min(lva1,iafa),min(idi2,jel))*sigc)/(siga-sigc)
  741.  409    CONTINUE
  742.         IPOS=IPOS+1
  743.         JTT(IPOS)=ISUP
  744.       IF (SIGRAV*SIGA.GT.0.) THEN
  745.         ISEG(IPOS)=1
  746.       ENDIF
  747.       ENDIF
  748.       IF (IPOS.GT.4) WRITE (6,*) ' ANOMALIE IPOS VOLUME ',IPOS
  749.       IF (IPOS.LT.3) WRITE (6,*) ' ANOMALIE IPOS VOLUME ',IPOS
  750.       IF (IPOS.LT.3) GOTO 160
  751.       ITR(**)=JTT(1)
  752.       ITR(**)=JTT(2)
  753.       ITR(**)=JTT(3)
  754.       ITR(**)=IPT1.ICOLOR(JEL)
  755.       if (mcham.ne.0) then
  756.        xtr(**)=xsegm(jtt(1))
  757.        xtr(**)=xsegm(jtt(2))
  758.        xtr(**)=xsegm(jtt(3))
  759.       endif
  760. *     WRITE (6,*) ' TRIANGLE (CAS 1) DU PLAN DE COUPE '
  761. *     WRITE (6,*) JTT(1),JTT(2),JTT(3)
  762. C  TRIANGLE DU PLAN DE COUPE
  763.       IF (IPOS.EQ.3) THEN
  764.        MCOUP(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)+8
  765.       ELSE
  766.        MCOUP(**)=ISEG(1)+2*ISEG(2)+4*0      +8
  767.       ENDIF
  768.  403  CONTINUE
  769.       IF (IPOS.GE.4) THEN
  770.       ITR(**)=JTT(1)
  771.       ITR(**)=JTT(3)
  772.       ITR(**)=JTT(4)
  773.       ITR(**)=IPT1.ICOLOR(JEL)
  774.       if (mcham.ne.0) then
  775.        xtr(**)=xsegm(jtt(1))
  776.        xtr(**)=xsegm(jtt(3))
  777.        xtr(**)=xsegm(jtt(4))
  778.       endif
  779. *     WRITE (6,*) ' TRIANGLE (CAS 2) DU PLAN DE COUPE '
  780. *     WRITE (6,*) JTT(1),JTT(3),JTT(4)
  781. C  TRIANGLE DU PLAN DE COUPE
  782.        MCOUP(**)=0      +2*ISEG(3)+4*ISEG(4)+8
  783.       ENDIF
  784.  402  CONTINUE
  785.  160  CONTINUE
  786.  161  CONTINUE
  787.   20  CONTINUE
  788.       NBELEM=JELN
  789.       SEGADJ IPT2
  790.       ISU=ISU+1
  791.       IPT8.LISOUS(ISU)=IPT2
  792.       if (mcham.ne.0) then
  793.        n1el=nbelem
  794.        segadj melva2
  795.        ipt8.lisref(isu)=melva2
  796.       endif
  797.       SEGDES IPT1
  798.   10  CONTINUE
  799.       IF (NBSOU.NE.0) SEGDES MELEME
  800.       NBSOUS=ISU+1
  801. *     IF (ITR(/1).EQ.0) NBSOUS=ISU
  802.       NBELEM=0
  803.       NBNN=0
  804.       NBREF=0
  805.       if (mcham.ne.0) nbref=nbsous
  806.       SEGADJ IPT8
  807.       MELEME=IPT8
  808. *     IF (ITR(/1).EQ.0) RETURN
  809. *      WRITE(6,*) ' NOMBRE DE TRIANGLE ENGENDREES DANS LE PLAN DE COUPE',
  810. *     $     ITR(/1)/4
  811. *      WRITE(6,*) ' NOMBRE DE TRIANGLE ENGENDREES HORS LE PLAN DE COUPE',
  812. *     $     ITRT(/1)/5
  813. *      WRITE(6,*) ' NOMBRE DE QUADRILA ENGENDREES HORS LE PLAN DE COUPE',
  814. *     $     IQUAT(/1)/6
  815.       NBELEM=ITR(/1)/4+ITRT(/1)/5+(IQUAT(/1)/6)*2
  816.       NBNN=3
  817.       NBREF=0
  818.       NBSOUS=0
  819.       SEGINI IPT4
  820.       IPT4.ITYPEL=4
  821.       if (mcham.ne.0) then
  822.        n1el=nbelem
  823.        n1ptel=nbnn
  824.        segini melva4
  825.       endif
  826.       ITRL1=ITR(/1)/4
  827.       DO 200 IEL=1,ITRL1
  828.       DO 201 NBNN=1,3
  829.       IPT4.NUM(NBNN,IEL)=ITR((IEL-1)*4+NBNN)
  830.       if (mcham.ne.0) melva4.velche(nbnn,iel)=xtr((IEL-1)*3+nbnn)
  831.  201  CONTINUE
  832.       IPT4.ICOLOR(IEL)=ITR(IEL*4)
  833.  200  CONTINUE
  834.       ITRL2=ITRT(/1)/5
  835.       DO 202 IEL=1,ITRL2
  836.       DO 203 NBNN=1,3
  837.       IPT4.NUM(NBNN,IEL+ITRL1)=ITRT((IEL-1)*5+NBNN)
  838.       if (mcham.ne.0)
  839.      >  melva4.velche(nbnn,iel+itrl1)=xtrt((IEL-1)*3+nbnn)
  840.  203  CONTINUE
  841.       IPT4.ICOLOR(IEL+ITRL1)=ITRT(IEL*5-1)
  842.       MCOUP(**)=ITRT(IEL*5)
  843.  202  CONTINUE
  844.       ITRL3=(IQUAT(/1)/6)
  845.       DO 204 IEL=1,ITRL3
  846.       IEL1=2*IEL-1
  847.       DO 205 NBNN=1,3
  848.       IPT4.NUM(NBNN,IEL1+ITRL1+ITRL2)=IQUAT((IEL-1)*6+NBNN)
  849.       if (mcham.ne.0)
  850.      >  melva4.velche(nbnn,iel1+itrl1+itrl2)=xquat((IEL-1)*4+nbnn)
  851.  205  CONTINUE
  852.       IPT4.ICOLOR(IEL1+ITRL1+ITRL2)=IQUAT(IEL*6-1)
  853.       MCOUP(**)=MOD(IQUAT(IEL*6),4)+16*(IQUAT(IEL*6)/16)
  854.       IEL2=2*IEL
  855.       IPT4.NUM(1,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+1)
  856.       IPT4.NUM(2,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+3)
  857.       IPT4.NUM(3,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+4)
  858.       IPT4.ICOLOR(IEL2+ITRL1+ITRL2)=IQUAT(IEL*6-1)
  859.       MCOUP(**)=0+MOD(IQUAT(IEL*6)/4,2)*2+MOD(IQUAT(IEL*6)/8,2)*4
  860.      #  +16*(IQUAT(IEL*6)/16)
  861.       if (mcham.ne.0) then
  862.        melva4.velche(1,iel2+itrl1+itrl2)=xquat((IEL-1)*4+1)
  863.        melva4.velche(2,iel2+itrl1+itrl2)=xquat((IEL-1)*4+3)
  864.        melva4.velche(3,iel2+itrl1+itrl2)=xquat((IEL-1)*4+4)
  865.       endif
  866.  204  CONTINUE
  867. *     WRITE (6,*) ' LISTE DES ELEMENTS RAJOUTES '
  868. *     DO 9834 J=1,NBELEM
  869. *     WRITE (6,*) (IPT4.NUM(I,J),I=1,3)
  870. *9834 CONTINUE
  871.       LISOUS(ISU+1)=IPT4
  872.       if (mcham.ne.0) lisref(ISU+1)=melva4
  873.       IF (MELEM2.NE.0) THEN
  874.          IF (MEDSAU.EQ.0) THEN
  875. C           ON COUPE MELEM2
  876.             mcham=0
  877.             MEDSAU=MELEME
  878.             MEDCOU=MCOUP
  879.             MELEME=MELEM2
  880.             SEGSUP ITR,ITRT,IQUAT
  881.             SEGINI MCOUP,ITR,ITRT,IQUAT
  882.             GOTO 5000
  883.          ELSE
  884.             MELEM2=MELEME
  885.             MELEME=MEDSAU
  886.             MCOU2=MCOUP
  887.             MCOUP=MEDCOU
  888.          ENDIF
  889.       ENDIF
  890.       SEGSUP ITR,ISEGM,IJGRAV,ITRT,IQUAT
  891.       if (mchamo.ne.0) segsup xtr,xtrt,xquat,xsegm
  892.       RETURN
  893.       END
  894.  
  895.  
  896.  
  897.  

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