Télécharger la plaquette Cast3M


fr en

Télécharger crcou2.eso

Retour à la liste

Numérotation des lignes : 

crcou2

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

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