C CCUBIC    SOURCE    PV        20/03/24    21:15:39     10554          
      SUBROUTINE CCUBIC
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C************************************************************************
C Ce sp transforme des elemnts quadratique d'un type pris
C dans la liste ci-dessous
C SEG3 TRI6 QUA8 CU20 PR15 TE10 PY13
C  3    6    10   15   17   24   26
C
C en les éléments correspondant MACRO (iso p2 ou q2) de la liste
C ci-dessous
C
C SEG3 TRI6 QUA9 CU27 PR18 TE15 PY14
C  3    6    11   33   ??   35   ??
C************************************************************************
-INC SMELEME
-INC SMCOORD

-INC PPARAM
-INC CCOPTIO
-INC SMLENTI
      SEGMENT TRAV
      INTEGER IFAC4(4,NBFTM)
      ENDSEGMENT
      SEGMENT TRAV1
      INTEGER NUMF4(9,NBFTM)
      ENDSEGMENT
      SEGMENT SITAB
      INTEGER ITAB(NPA,NBATM)
      ENDSEGMENT
      SEGMENT SITAC
      INTEGER ITAC(NBMAX,NNMAX)
      ENDSEGMENT
      SEGMENT SITAF
      INTEGER ITAF(NBFAX,MMMAX)
      ENDSEGMENT
      SEGMENT TRAV2
      INTEGER NUF(6,NBELT)
      ENDSEGMENT
      SEGMENT CARA
      INTEGER KM(6,NBSOU1)
      ENDSEGMENT


      DIMENSION IAR(3,35),IFR(4,20)
      DATA IAR /1,2,3,3,4,5,5,6,7,7,8,1,
     &          1,9,13,3,10,15,5,11,17,7,12,19,
     &          13,14,15,15,16,17,17,18,19,19,20,13,
     &          1,2,3,3,4,5,5,6,1,1,7,10,3,8,12,5,9,14,
     &          10,11,12,12,13,14,14,15,10,
     &          1,2,3,3,4,5,5,6,1,1,7,10,3,8,10,5,9,10,
     &          1,2,3,3,4,5,5,6,7,7,8,1,1,9,13,3,10,13,
     &          5,11,13,7,12,13/

      DATA IFR /1,6,9,5, 2,7,10,6, 3,8,11,7, 4,5,12,8,
     &          1,2,3,4, 9,10,11,12,
     &          1,5,7,4, 2,6,8,5, 3,4,9,6, 1,2,3,0, 7,8,9,0,
     &          1,2,3,0, 1,5,4,0, 2,6,5,0, 3,6,4,0,
     &          1,2,3,4, 1,6,5,0, 2,7,6,0, 3,8,7,0, 4,5,8,0/

C Nb de pts par face par type elt
      DIMENSION KNPF(6,7)
      DATA KNPF/1,1,0,0,0,0,
     &          3,3,3,0,0,0,
     &          3,3,3,3,0,0,
     &          8,8,8,8,8,8,
     &          8,8,8,6,6,0,
     &          6,6,6,6,0,0,
     &          8,6,6,6,6,0/

C Nb d'arretes par face par type elt
      DIMENSION KNAF(6,7)
      DATA KNAF/0,0,0,0,0,0,
     &          0,0,0,0,0,0,
     &          0,0,0,0,0,0,
     &          4,4,4,4,4,4,
     &          4,4,4,3,3,0,
     &          3,3,3,3,0,0,
     &          4,3,3,3,3,0/

C?    DIMENSION I12(35),J12(11)
C?    DATA I12  /1,2,3,4,5,6,7,8,10,12,14,16,19,20,21,22,23,24,25,26,
C?   &           1,2,3,4,5,6,8,10,12,15,16,17,18,19,20/
C?    DATA J12  /11,13,15,17,9,27,
C?   &           9 ,11,13,7 ,21/

      DIMENSION INF(8,20)
C CU20
      DATA INF/1,2,3,10,15,14,13,9,  3,4,5,11,17,16,15,10,
     &         5,6,7,12,19,18,17,11, 7,8,1,9,13,20,19,12,
     &         1,2,3,4,5,6,7,8,      13,14,15,16,17,18,19,20,
C PR15
     &         1,2,3,8 ,12,11,10,7,  3,4,5,9 ,14,13,12,8 ,
     &         5,6,1,7,10,15,14,9 ,  1,2,3,4,5,6,0,0,
     &         10,11,12,13,14,15,0,0,
C TE15
     &  1,2,3,4,5,6,0,0,  1,2,3,8,10,7,0,0,  3,4,5,9,10,8,0,0,
     &  1,6,5,9,10,7,0,0,
C PY15
     &  1,2,3,4,5,6,7,8,  1,2,3,10,13,9,0,0,  3,4,5,11,13,10,0,0,
     &  5,6,7,12,13,11,0,0,  7,8,1,9,13,12,0,0/





      DIMENSION NUA(12),NUMA(4),XA(3,21),KTA(7,9)
      DATA KTA/3,6,10,15,17,24,26,
     &         3,6,11,33,00,35,00,
C nb de faces
     &         2,3,4 ,6 ,5 ,4 ,5 ,
C nb d'arretes
     &         0,3,4 ,12,9 ,6 ,8,
C Idec
     &         0,0,0 ,0 ,12,21,27,
C Idc3
     &         0,0,0 ,0 ,6 ,11,15,
C Nbnn
     &         3,6,9 ,27,18,15,14,
C NP
     &         3,6,8 ,20,15,10,13,
C IDF
     &         0,0,0 ,0 ,6 ,11,15/

C?Numero pt centre N18
C?   &         0,0,0 ,18,14,14,18/
C?Idc2
C?   &         0,0,0 ,0 ,20,35,45,

C SEG3 TRI6 QUA8 CU20 PR15 TE10 PY13
C  3    6    10   15   17   24   26
C SEG3 TRI6 QUA9 CU27 PR18 TE15 PY14
C  3    6    11   33   ??   35   ??



C*************************************************************

      CALL LIROBJ('MAILLAGE',MELEME,1,IRET)
      IF(IRET.EQ.0)RETURN

      SEGACT MELEME
      NBSOU1=LISOUS(/1)
      IF(NBSOU1.EQ.0)NBSOU1=1

      SEGINI CARA

      NBELT=0
      DO 11 L=1,NBSOU1
      IPT1=MELEME
      IF(NBSOU1.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
      ITYP=IPT1.ITYPEL
      IK=0
      DO 111 I=1,7
      IF(ITYP.EQ.KTA(I,1))IK=I
 111  CONTINUE

      IF(IK.EQ.0)THEN
      CALL ERREUR(29)
      RETURN
      ENDIF


      NP    =IPT1.NUM(/1)
      NBELEM=IPT1.NUM(/2)
      IF(IK.EQ.3.OR.IK.EQ.4)THEN
      NBELT=NBELT+NBELEM
      ENDIF
      KM(1,L)=NBELEM
      KM(2,L)=IPT1
      KM(3,L)=IK
C nb d'aretes par element
      KM(4,L)=KTA(IK,4)
C nb de faces
      KM(5,L)=KTA(IK,3)
      SEGDES IPT1
 11   CONTINUE
      segact mcoord*mod
      NBPT=NBPTS    
      JG=NBPT
      SEGINI MLENTI,MLENT1

      NPA=3
      NBATM=5*NBELT+500
      SEGINI SITAB

      NBMAX=4+1
      NNMAX=3*NBELT+300
      NBFAX=4+1
      MMMAX=3*NBELT+300
      SEGINI SITAF,SITAC

      NBFTM=5*NBELT+500
      SEGINI TRAV,TRAV1,TRAV2

      SEGACT MELEME

      NBFT4=0
      NBAT=0
      NN=0
      MM=0

      NK=0
      DO 1 L=1,NBSOU1
      IK=KM(3,L)

      IPT1=MELEME
      IF(NBSOU1.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1
      NBELEM=IPT1.NUM(/2)
      NP    =IPT1.NUM(/1)
C     write(6,*)' ITYP,NBELEM=',ITYP,NBELEM

      IF(IK.EQ.1)GO TO 1
      IF(IK.EQ.2.OR.IK.EQ.3)THEN
      NK=NK+NBELEM
      GO TO 1
      ENDIF

      IDEC=KTA(IK,5)
      IDC3=KTA(IK,6)
      IDF =KTA(IK,9)
      NBA=KTA(IK,4)
      NBF=KTA(IK,3)
C     write(6,*)' NBA,NBF,NPF=',NBA,NBF

      DO 2 K=1,NBELEM
      NK=NK+1

      DO 5 NA=1,NBA
      N1=IPT1.NUM(IAR(1,NA+IDEC),K)
      N2=IPT1.NUM(IAR(3,NA+IDEC),K)
      IF(N1.GT.N2)THEN
      NM=N2
      N2=N1
      N1=NM
      ENDIF
      NM=IPT1.NUM(IAR(2,NA+IDEC),K)

      IF(LECT(N1).EQ.0)THEN
C le sommet n'a pas encore été touché
      NBAT=NBAT+1
      IF(NBAT.GT.NBATM)THEN
      write(6,*)' Taille ITAB 2ème dime insuffisante NBATM=',NBATM
      NBATM=NBATM+NBELEM
C     write(6,*)' NBATM=',nbatm,' NPA=',npa,nbat
      SEGADJ SITAB
      ENDIF

      NN=NN+1

      IF(NN.GT.NNMAX)THEN
      write(6,*)' Taille ITAC 2eme dime insuffisante NN=',NN
      NNMAX=NNMAX+NBELEM
      SEGADJ SITAC
      ENDIF

      LECT(N1)=NN
      ITAC(1,NN)=1
      ITAC(1+1,NN)=NBAT
C     write(6,*)' NBAT=',nbat,' nn=',nn,' NB=',ITAC(1,NN),nbat
      NUA(NA)=NBAT
      ITAB(1,NBAT)=N1
      ITAB(2,NBAT)=NM
      ITAB(3,NBAT)=N2
      GO TO 5
      ENDIF

C On cherche si l'arrete existe deja dans la table ITAC
      NN1=LECT(N1)
      NB=ITAC(1,NN1)
      DO 31 II=1,NB
      I=ITAC(II+1,NN1)
        IF(N2.EQ.ITAB(3,I).AND.NM.EQ.ITAB(2,I))THEN
        NUA(NA)=I
        GO TO 5
        ENDIF
 31   CONTINUE

      IF(NB.LT.(NBMAX-1))THEN
C l'arrete n existe pas dans la table ITAC qui n'est pas pleine
C On peut donc considerer que l'arrete est nouvelle
      NBAT=NBAT+1
      ITAC(1,NN1)=NB+1
      ITAC(NB+2,NN1)=NBAT
C     write(6,*)' NBAT=',nbat,' nn1=',nn1,' NB=',ITAC(1,NN1),nbat
      NUA(NA)=NBAT
      ITAB(1,NBAT)=N1
      ITAB(2,NBAT)=NM
      ITAB(3,NBAT)=N2
      GO TO 5
      ENDIF

      write(6,*)' Taille ITAC 1ere dime insuffisante NB=',NB
      NBMAX=NBMAX+2
      SEGADJ SITAC
C On fait une recherche parmis toutes les arretes existantes

      DO 3 I=1,NBAT
        IF(N2.EQ.ITAB(3,I).AND.NM.EQ.ITAB(2,I))THEN
        NUA(NA)=I
        GO TO 5
        ENDIF
 3    CONTINUE

      NBAT=NBAT+1
      ITAC(1,NN1)=NB+1
      ITAC(NB+2,NN1)=NBAT
      NUA(NA)=NBAT
      ITAB(1,NBAT)=N1
      ITAB(2,NBAT)=NM
      ITAB(3,NBAT)=N2

 5    CONTINUE

      DO 6 NF=1,NBF
C nb de pts par faces
      NPF=KNPF(NF,IK)
      IF(NPF.NE.8)GO TO 6
C nb d'arretes par faces
      NAF=KNAF(NF,IK)

      DO 61 NFA=1,NAF
      NUMA(NFA)=NUA(IFR(NFA,NF+IDF))
 61   CONTINUE

      CALL ORDOTA(NUMA,NAF)
      IF(MLENT1.LECT(NUMA(1)).EQ.0)THEN
C l'arrete n'a pas encore été touché
      NBFT4=NBFT4+1

      IF(NBFT4.GT.NBFTM)THEN
      write(6,*)' Taille TRAV NBFT4 insuffisante NBFT4=',NBFT4
      NBFTM=NBFTM+NBELEM
      SEGADJ TRAV,TRAV1
      ENDIF

      NUMF4(9,NBFT4)=NPF
      DO 621 I=1,NPF
      NUMF4(I,NBFT4)=IPT1.NUM(INF(I,NF+IDC3),K)
 621  CONTINUE
      MM=MM+1

      IF(MM.GT.MMMAX)THEN
      write(6,*)' Taille ITAF 2eme dime insuffisante MM=',MM
      MMMAX=MMMAX+NBELEM
      SEGADJ SITAF
      ENDIF

      MLENT1.LECT(NUMA(1))=MM
      ITAF(1,MM)=1
      ITAF(1+1,MM)=NBFT4
C     write(6,*)' NBFT4=',nbft4,' mm=',mm,' NB=',ITAF(1,mm)
      NUF(NF,NK)=NBFT4
      CALL RSETI(IFAC4(1,NBFT4),NUMA,NAF)
      GO TO 6
      ENDIF

C On cherche si la face existe déja dans la table ITAF
      MM1=MLENT1.LECT(NUMA(1))
      NB=ITAF(1,MM1)
      DO 631 II=1,NB
      I=ITAF(II+1,MM1)
        IF( NUMA(2).EQ.IFAC4(2,I).AND.NUMA(3).EQ.IFAC4(3,I)
     & .AND.NUMA(1).EQ.IFAC4(1,I))THEN
        NUF(NF,NK)=I
        GO TO 6
        ENDIF
 631  CONTINUE

      IF(NB.LT.(NBFAX-1))THEN
C la face n'existe pas dans la table ITAF qui n'est pas pleine
C On peut donc considerer que la face est nouvelle
      NBFT4=NBFT4+1

      IF(NBFT4.GT.NBFTM)THEN
      write(6,*)' Taille TRAV NBFT4 insuffisante NBFT4=',NBFT4
      NBFTM=NBFTM+NBELEM
      write(6,*)' NBFTM=',NBFTM,NBELEM
      SEGADJ TRAV,TRAV1
      ENDIF

      NUMF4(9,NBFT4)=NPF
      DO 622 I=1,NPF
      NUMF4(I,NBFT4)=IPT1.NUM(INF(I,NF+IDC3),K)
 622  CONTINUE
      ITAF(1,MM1)=NB+1
      ITAF(NB+2,MM1)=NBFT4
C     write(6,*)' NBFT4=',NBFT4,' mm1=',mm1,' NB=',ITAF(1,mm1)
      NUF(NF,NK)=NBFT4
      CALL RSETI(IFAC4(1,NBFT4),NUMA,NAF)
      GO TO 6
      ENDIF

      write(6,*)' Taille ITAF 1ere dime insuffisante NB=',NB
      NBFAX=NBFAX+2
      SEGADJ SITAF
C On fait une recherche parmis toutes les faces existantes


      IF(NAF.EQ.4)THEN
      DO 7 I=1,NBFT4
        IF( NUMA(2).EQ.IFAC4(2,I).AND.NUMA(3).EQ.IFAC4(3,I)
     & .AND.NUMA(1).EQ.IFAC4(1,I))THEN
        NUF(NF,NK)=I
        GO TO 6
        ENDIF
 7    CONTINUE
      ELSEIF(NAF.EQ.3)THEN
      DO 71 I=1,NBFT4
        IF( NUMA(2).EQ.IFAC4(2,I).AND.NUMA(3).EQ.IFAC4(3,I))THEN
        NUF(NF,NK)=I
        GO TO 6
        ENDIF
 71   CONTINUE
      ENDIF


 64   CONTINUE
      NBFT4=NBFT4+1

      IF(NBFT4.GT.NBFTM)THEN
      write(6,*)' Taille TRAV NBFT4 insuffisante NBFT4=',NBFT4
      NBFTM=NBFTM+NBELEM
      SEGADJ TRAV,TRAV1
      ENDIF

      ITAF(1,MM1)=NB+1
      ITAF(NB+2,MM1)=NBFT4
      NUMF4(9,NBFT4)=NPF
      DO 623 I=1,NPF
      NUMF4(I,NBFT4)=IPT1.NUM(INF(I,NF+IDC3),K)
 623  CONTINUE
      NUF(NF,NK)=NBFT4
      CALL RSETI(IFAC4(1,NBFT4),NUMA,NAF)

 6    CONTINUE
 2    CONTINUE
 1    CONTINUE

C**********************************************************

      SEGSUP SITAB,SITAC,SITAF,MLENTI,MLENT1,TRAV

C**********************************************************

      SEGACT MELEME

      NK=0
      DO 80 L=1,NBSOU1
      IK=KM(3,L)
      NBF=KTA(IK,3)
      IF(IK.EQ.1)GO TO 80

      IPT1=MELEME
      IF(NBSOU1.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1

      NBELEM=KM(1,L)
      NBNN  =KTA(IK,7)
      NP    =KTA(IK,8)

      NBSOUS=0
      NBREF=0
      SEGINI IPT2
      KM(2,L)=IPT2
      IPT2.ITYPEL=KTA(IK,2)

      IF(IK.GE.4)THEN
C CU27 & PR21
      IDC3=KTA(IK,6)

      DO 83 K=1,NBELEM
      NK=NK+1
      DO 8 I=1,NP
      IPT2.NUM(I,K)=IPT1.NUM(I,K)
 8    CONTINUE
      DO 81 I=1,NBF
      IPT2.NUM(I+NP,K)=NBPT+NBELT+NUF(I,NK)
 81   CONTINUE
      IPT2.NUM(NBNN,K)=NBPT+NK
 83   CONTINUE

      ELSEIF(IK.EQ.2.OR.IK.EQ.3)THEN
      IDC=7
      IF(IK.EQ.3)IDC=9
      NP=NBNN-1
      DO 84 K=1,NBELEM
      NK=NK+1
      DO 88 I=1,NP
      IPT2.NUM(I,K)=IPT1.NUM(I,K)
 88   CONTINUE
      IPT2.NUM(IDC,K)=NBPT+NK
 84   CONTINUE

      ENDIF

      SEGDES IPT1,IPT2
 80   CONTINUE

C     write(6,*)' NUF '
C     DO 783 K=1,NBELT
C     write(6,1011)K,(NUF(i,k),i=1,6)
C783  continue

C     do 784 k=1,nbft4
C     write(6,1011)K,(NUMF4(i,k),i=1,8)
C784  continue

C     write(6,*)' NBPT=',nbpt
C     write(6,*)' NBAT=',nbat
C     write(6,*)' NBELT=',NBELT,' NBFT4=',nbft4

      NBPTS=NBPT+NBELT+NBFT4
      IF(NBPTS.GT.NBPT)SEGADJ MCOORD

C************** On calcule les coordonnées des points ***********



      IF(NBFT4.NE.0)THEN

      DO 23 J=1,NBFT4

      NPF=NUMF4(9,J)
      DO 21 I=1,NPF
      N1=NUMF4(I,J)
      DO 21 M=1,IDIM
      XA(M,I)=XCOOR((N1-1)*(IDIM+1)    +M)
 21   CONTINUE

      CALL FFFACE(XA,NPF)

      N9=NBPT+NBELT+J
      DO 22 M=1,IDIM
      XCOOR((N9-1)*(IDIM+1)    +M)=XA(M,NPF+1)
 22   CONTINUE
 23   CONTINUE
      ENDIF


C     write(6,*)' On calcule les coordonnées du pt centre ',NBELT
      IF(NBELT.NE.0)THEN
      SEGACT MELEME
      NK=0
      DO 90 L=1,NBSOU1

      IK=KM(3,L)
      IF(IK.EQ.1)GO TO 90
      IPT1=MELEME
      IF(NBSOU1.NE.1)IPT1=LISOUS(L)
      SEGACT IPT1

      NBELEM=KM(1,L)
      NBNN  =KTA(IK,7)
      NP    =KTA(IK,8)

      IPT2=KM(2,L)

      DO 24 K=1,NBELEM
      NK=NK+1
      CALL INITD(XA(1,21),IDIM,0.D0)
      DO 25 I=1,NP
      N1=IPT1.NUM(I,K)
      DO 25 M=1,IDIM
      XA(M,I)=XCOOR((N1-1)*(IDIM+1)    +M)
      XA(M,21)=XA(M,21)+XA(M,I)
 25   CONTINUE

      N18=NBPT+NK
C     write(6,*)' *** NK=',nk,NBELT,L
      DO 26 M=1,IDIM
      XCOOR((N18-1)*(IDIM+1)    +M)=XA(M,21)/FLOAT(NP)
 26   CONTINUE
 24   CONTINUE

 90   CONTINUE

      ENDIF


C     write(6,*)' NBPT=',nbpt
C     write(6,*)' NBAT=',nbat
C     do 116 l=1,nbat
C     write(6,1011)L,(itab(i,l),i=1,3)
C116  continue

C     write(6,*)' LECT '
C     write(6,1001)(lect(ii),ii=1,nbpt)

C     write(6,*)' ITAC NN=',NN
C     do 118 l=1,NN
      nb=itac(1,l)
C     write(6,1011)L,itac(1,l),(itac(i+1,l),i=1,nb)
C118  continue

C     write(6,*)' NBFT4=',nbft4

C     do 117 l=1,nbft4
C     write(6,1011)L,(ifac4(i,l),i=1,4)
C117  continue

C     write(6,*)' MLENT1.LECT '
C     write(6,1001)(mlent1.lect(ii),ii=1,nbpt)

C     write(6,*)' ITAF MM=',MM
C     do 119 l=1,MM
C     nb=itaf(1,l)
C     write(6,1011)L,itaf(1,l),(itaf(i+1,l),i=1,nb)
C119  continue

      SEGSUP TRAV1
      IF(NBSOU1.EQ.1)THEN
      IPT3=KM(2,1)
      ELSE
      NBSOUS=NBSOU1
      NBELEM=0
      NBNN=0
      NBREF=0
      SEGINI IPT3
      DO 785 L=1,NBSOU1
      IPT3.LISOUS(L)=KM(2,L)
 785  CONTINUE
      ENDIF

C?    SEGACT IPT3
C?    nbref=IPT3.LISREF(/1)
C?    write(6,*)' nbref,meleme=',nbref,meleme
C?    IPT3.LISREF(1)=MELEME
      SEGDES IPT3

      CALL ECROBJ('MAILLAGE',IPT3)

      SEGSUP CARA,TRAV2

      RETURN
 1011 FORMAT('L=',I3,4X,15(1X,I5))
 1001 FORMAT(20(1X,I5))
 1002 FORMAT(10(1X,1PE11.4))
      END