Télécharger raff.eso

Retour à la liste

raff
C RAFF      SOURCE    OF166741  26/04/22    21:15:02     12520          
      SUBROUTINE RAFF
C----------------------------------------------------------------------C
C                            Operateur RAFF                            C
C----------------------------------------------------------------------C
C
C Syntaxe 1 : GEO2 = RAFF GEO1 CHPO1 ;
C
C Syntaxe 2 : LRE2 = RAFF LRE1 ENT1 ;
C
C----------------------------------------------------------------------C
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
c-INC CCREEL
 
-INC SMCOORD
-INC SMELEME
-INC SMCHPOI
-INC SMMODEL
-INC SMCHAML
 
      INTEGER MCOL,NCOL
C--------------------------------
C B) VARIABLES RELATIVES A LA 2D !
C--------------------------------
      SEGMENT KARETE(NBNDS,NCOL)
      SEGMENT KARET2(NBNDS,NCOL)
      SEGMENT KMILIE(NBNDS,NCOL)
      SEGMENT KARPOS(NBNDS)
C
C--------------------------------
C C) VARIABLES RELATIVES A LA 3D !
C--------------------------------
      SEGMENT IRELA(NBNDS)
      SEGMENT JPLANS(JPLA1,JPLA2)
      SEGMENT JPLAN2(JPLA1,JPLA2)
      SEGMENT JPLAN3(JPLA1,JPLA2)
      SEGMENT JPLCOM(JPLA1)
      SEGMENT JNOEFA(JNBFA,5)
      SEGMENT IWORK2(JNBSOM)
      SEGMENT JFARAF(JNBFA,LLLL)
      SEGMENT IWORK1(IJKLMN)
C--------------------------
C A) VARIABLES 'GENERALES' !
C--------------------------
      SEGMENT ICPR(nbpts,2)
      SEGMENT LTYPE(NBLIS,2)
      SEGMENT LTYPE2(NBLIS2,4)
      SEGMENT LTYPE3(NBLIS3,2)
      SEGMENT LTYPE4(NBLIS3,2)
      SEGMENT IDENSI
        REAL*8 DENSI(NBNDS)
      ENDSEGMENT
      SEGMENT XDEN(nbpts)
      XPETI= 1E-30
C
C---------------------------------------------------
C D) VARIABLES RELATIVES AUX MAILLAGES DE RELATIONS !
C---------------------------------------------------
      POINTEUR IPT10.MELEME
      POINTEUR IPT11.MELEME
      POINTEUR IPT12.MELEME
      POINTEUR IPT13.MELEME
 
      POINTEUR IPT14.MELEME
      POINTEUR IPT15.MELEME
      POINTEUR IPT16.MELEME
      SEGMENT MAIREL(NBMAIL)
C=======================================================================
C                  Lecture des donnees utilisateurs
C=======================================================================
C
C     SYNTAXE 2 : cas du LISTREEL
C     ---------------------------
      CALL LIROBJ('LISTREEL',ILRE1,0,IRETOU)
      IF (IERR.NE.0) RETURN
C
      IF (IRETOU.NE.0) THEN
        ICAS = 0
        CALL LIRENT(KR1,0,IRETOU)
        IF (IERR.NE.0) RETURN
        IF (IRETOU.EQ.0) THEN
          CALL LIRREE(XR1,1,IRETOU)
          IF (IERR.NE.0) RETURN
          ICAS = 2
        ELSE
          ICAS = 1
        ENDIF
        CALL RAFLRE(ILRE1,KR1,XR1,ICAS,ILRE2)
        IF (IERR.NE.0) RETURN
        IF (ILRE2.NE.0) THEN
          CALL ECROBJ('LISTREEL',ILRE2)
        ELSE
          CALL ERREUR(5)
        ENDIF
        RETURN
      ENDIF
 
C     SYNTAXE 1 : cas du MAILLAGE
C     ---------------------------
C On lit le maillage a raffiner et le champoint de densite a affecter
c ipt1 : maillage initial a raffiner
c mchpo1 : chpoint de densite voulu (qui devrait etre defini sur ipt1)
      CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
      CALL ACTOBJ('MAILLAGE',IPT1,1)
      IF (IERR.NE.0) RETURN
      CALL LIROBJ('CHPOINT',MCHPO1,1,IRETOU)
      CALL ACTOBJ('CHPOINT',MCHPO1,1)
      IF (IERR.NE.0) RETURN
 
C=======================================================================
C                   Initialisations - Declarations
C=======================================================================
      NCOL=0
      MCOL=0
      NBOUCL=0
      NPTINI=nbpts
 
      IJKLMN=4
      SEGINI IWORK1
      IWORK1(1)=4
      IWORK1(2)=8
      IWORK1(3)=6
      IWORK1(4)=10
 
      SEGINI,XDEN
 
C=======================================================================
C                       Tests sur les donnees
C=======================================================================
C On teste pour savoir si le maillage IPT1 (maillage 1) a un type
C d'elements implemente.
c TRI3 TRI6 QUA4 QUA8 CUB8 CUB26 PRY6 PR15 TET4 TET10 SURE
C les pyramides 5 et 13 ont été débranchés
c creation du tableau
      NBSOU1 = IPT1.LISOUS(/1)
      NBLIS  = MAX(1,NBSOU1)
      NSURE=0
      SEGINI LTYPE
      IPT2=IPT1
      DO  1 INB=1,NBLIS
        IF (NBSOU1.NE.0) IPT2=IPT1.LISOUS(INB)
        NT=IPT2.ITYPEL
        IF (NT.EQ.48) NSURE=NSURE+1
        LTYPE(INB,1)=INB
        LTYPE(INB,2)=NT
        IF (NT.EQ.16) THEN
          write (*,*) 'Attention le raffinement PRI6 a debrancher !'
        ENDIF
        IF (NT.EQ.23) THEN
          write (*,*) 'Attention le raffinement TET4 a debrancher !'
        ENDIF
        IF ((NT.NE.4).AND.(NT.NE.6).AND.(NT.NE.8).AND.(NT.NE.10).AND.
     +   (NT.NE.14).AND.(NT.NE.15).AND.(NT.NE.16).AND.(NT.NE.17).AND.
     +   (NT.NE.23).AND.(NT.NE.24).AND.(NT.NE.48)) THEN
          WRITE(*,*) 'Erreur, type d element non pris en compte'
          RETURN
        ENDIF
c        write(*,*) 'type d elem du maillage de depart', LTYPE(INB,1),
c     +  LTYPE(INB,2)
 
C        IF ((NT.NE.4).AND.(NT.NE.6).AND.(NT.NE.8).AND.(NT.NE.10).AND.
C     +   (NT.NE.14).AND.(NT.NE.15).AND.(NT.NE.16).AND.(NT.NE.17).AND.
C     +   (NT.NE.23).AND.(NT.NE.24).AND.(NT.NE.48).AND.(NT.NE.25).AND.
C     +   (NT.NE.26)) RETURN
        IF ((NBSOU1.EQ.0).AND.(NT.EQ.48)) RETURN
   1  CONTINUE
 
C=======================================================================
C                  Separation du maillage initial
C=======================================================================
C Si le maillage initial comporte des lisous de type 48, on le separe
C en deux maillages dont un ne sera compose que des lisous de type 48.
C Maillage d'elements de type 48 : IPT6
C Maillage sans elements de type 48 : IPT5
C IPT1 reste le maillage initial complet
      LCONF2=0
      IF (NSURE.EQ.0) GOTO 7
      LCONF2=1
   7  CONTINUE
C-------------------------------
C A) Creation du maillage IPT6 !
C-------------------------------
C  A.1/** Cas ou il n'y a qu'un lisous de type 48
      IF (NSURE.EQ.1) THEN
        DO 2 K=1,LTYPE(/1)
          IF (LTYPE(K,2).EQ.48) IPT6=IPT1.LISOUS(LTYPE(K,1))
   2    CONTINUE
      ENDIF
C
C  A.2/** Cas ou il y a plusieurs lisous de type 48
      IF (NSURE.GT.1) THEN
        NBREF=0
        NBSOUS=NSURE
        NBNN=0
        NBELEM=0
        SEGINI IPT6
        KLISOU=0
        DO 3 K=1,LTYPE(/1)
          IF (LTYPE(K,2).EQ.48) THEN
            KLISOU=KLISOU+1
            IPT6.LISOUS(KLISOU)=IPT1.LISOUS(LTYPE(K,1))
          ENDIF
   3    CONTINUE
      ENDIF
      IF (NSURE.EQ.0) THEN
        NBREF=0
        NBSOUS=NSURE
        NBNN=0
        NBELEM=0
        SEGINI IPT6
      ENDIF
C
C------------------------------
C B) Creation de IPT5 !
C------------------------------
C
      NBSOUS=LTYPE(/1)-NSURE
      NBREF=0
      NBNN=0
      NBELEM=0
C B.0/** Cas ou il n'y a qu'un lisous en tout
      IF ((NBSOUS.EQ.1).AND.(NSURE.EQ.0)) THEN
        IPT5=IPT1
      ENDIF
C B.1/** Cas ou il n'y a qu'un lisous de type different de 48
      IF ((NBSOUS.EQ.1).AND.(NSURE.NE.0)) THEN
        DO 4 K=1,LTYPE(/1)
          IF (LTYPE(K,2).NE.48) IPT5=IPT1.LISOUS(LTYPE(K,1))
   4    CONTINUE
      ENDIF
C
C B.2/** Cas ou il y a plusieurs lisous de types differents de 48
      IF (NBSOUS.GT.1) THEN
        SEGINI IPT5
        KLISOU=0
        DO 5 K=1,LTYPE(/1)
          IF (LTYPE(K,2).NE.48) THEN
            KLISOU=KLISOU+1
            IPT5.LISOUS(KLISOU)=IPT1.LISOUS(LTYPE(K,1))
          ENDIF
   5    CONTINUE
      ENDIF
 
c compteur de relaton globale
      LCONF3=LCONF2
C=======================================================================
C         Debut de la boucle sur les itérations de Raffinement
C=======================================================================
   9  CONTINUE
 
c       write(*,*) 'RAFF: ITERATION boucle 9',NBOUCL
c       write(*,*) 'ipt5 : ', ipt5
c       write(*,*) 'ipt6 : ', ipt6
c       write(*,*) 'lconf3 : ', lconf3
C
C=======================================================================
C         Renumerotation locale des noeuds du maillage
C=======================================================================
C ICPR(I,1)=J : Le Ieme noeud global est le Jeme local
C ICPR(I,2)=J : Le Ieme noeud global appartient a J elements
 
c*      SEGACT IPT5
 
      NBOUCL=NBOUCL+1
      SEGINI ICPR
      NBNDS=0
      NCOL=0
      NBS5 = IPT5.LISOUS(/1)
      IPT2=IPT5
      DO 11 INB=1,MAX(1,NBS5)
        IF (NBS5.NE.0) THEN
          IPT2=IPT5.LISOUS(INB)
c*          SEGACT IPT2
        ENDIF
        IF (IPT2.ITYPEL.EQ.48) GOTO 11
        DO 10 J=1,IPT2.NUM(/2)
        DO 110 I=1,IPT2.NUM(/1)
          ITMP=IPT2.NUM(I,J)
          IF (ICPR(ITMP,1).NE.0) THEN
            ICPR(ITMP,2)=ICPR(ITMP,2)+1
            NCOL=MAX(NCOL,ICPR(ITMP,2))
          ELSE
            NBNDS=NBNDS+1
            ICPR(ITMP,1)=NBNDS
            ICPR(ITMP,2)=1
            NCOL=MAX(NCOL,ICPR(ITMP,2))
          ENDIF
 110    CONTINUE
  10    CONTINUE
  11  CONTINUE
      IF (IDIM.EQ.2) NCOL=NCOL*2
      IF (IDIM.EQ.3) NCOL=NCOL*4
 
C=======================================================================
C           Prise en compte des noeuds de relations
C=======================================================================
C IRELA(I,1) = 1 si le noeud I (en numerotation locale) est un noeud
C                   de relation ("hanging node")
C IRELA(I,1) = 0 sinon
      SEGINI IRELA
      IF (LCONF3.EQ.0) GOTO 1200
      SEGACT IPT6
      DO 1218 LAA=1,MAX(1,IPT6.LISOUS(/1))
        IF (IPT6.LISOUS(/1).EQ.0) THEN
          IPT2=IPT6
        ELSE
          IPT2=IPT6.LISOUS(LAA)
        ENDIF
        SEGACT IPT2
        IF (IPT2.ITYPEL.NE.48) GOTO 1218
        DO 1220 J=1,IPT2.NUM(/2)
          IRELA(ICPR(IPT2.NUM(1,J),1))=1
 1220   CONTINUE
 1218 CONTINUE
 1200 CONTINUE
 
C
C=======================================================================
C           Prise en compte du champoint de densite
C=======================================================================
C On affecte une densite a chacun des noeuds du maillage. On calcule
C cette densite a partir des densites existantes dans le champoint, a
C l'aide de la fonction RAFDEN.
C Remarque : on ne prend pas en compte les densites qui pourraient
C exister dans le tableau XCOOR.
 
C------------------------------------------
C A) Recuperation des valeurs du champoint !
C------------------------------------------
C XDEN(I,1) = densite affectee au noeud I (en numerotation globale)
      IF (NBOUCL.EQ.1) THEN
      SEGACT MCHPO1
      IF (MCHPO1.IPCHP(/1).NE.1) RETURN
      MSOUPO=MCHPO1.IPCHP(1)
      SEGACT MSOUPO
      IPT3=IGEOC
      MPOVA1=IPOVAL
      SEGACT IPT3
      SEGACT MPOVA1
      IF (MPOVA1.VPOCHA(/2).NE.1) RETURN
      NBPTCH=IPT3.NUM(/2)
c       SEGINI ICHPO
      DO 13 I=1,NBPTCH
        J=IPT3.NUM(1,I)
c         XR(I)=XCOOR((J-1)*(IDIM+1)+1)
c         YR(I)=XCOOR((J-1)*(IDIM+1)+2)
c         IF (IDIM.EQ.3) ZR(I)=XCOOR((J-1)*(IDIM+1)+3)
c         DEN(I)=MPOVA1.VPOCHA(I,1)
c           if(J.eq.115.or.J.eq.321.or.
c      .       J.eq.6261.or.J.eq.5815) then
C             write(*,*) 'RAFF: #',J,'DEN(',I,'/',NBPTCH,')=',DEN(I)
c           endif
        XDEN(J)=MPOVA1.VPOCHA(I,1)
  13  CONTINUE
C     Rem ;seulement pour les noeud "peres",
C     pour les "fils" raff2d remplit XDEN
      ENDIF
C
C--------------------------------------------------------------
C B) Affectation d'une densite a chacun des noeuds du maillage !
C--------------------------------------------------------------
C IDENSI(I,1) = densite affectee au noeud I (en numerotation locale)
      SEGINI IDENSI
      DO 15 INB=1,MAX(1,IPT5.LISOUS(/1))
        IF (IPT5.LISOUS(/1).EQ.0) THEN
          IPT2=IPT5
        ELSE
          IPT2=IPT5.LISOUS(INB)
          SEGACT IPT2
        ENDIF
c        write(*,*) 'ipt2 : ', ipt2, 'type', IPT2.ITYPEL
c        write(*,*)  IPT2.NUM(/2), IPT2.NUM(/1)
c        IF (IPT2.ITYPEL.EQ.48) GOTO 15
        DO 14 J=1,IPT2.NUM(/2)
        DO 114 I=1,IPT2.NUM(/1)
          NGLOB=IPT2.NUM(I,J)
          NLOC=ICPR(NGLOB,1)
c           XPT=XCOOR((NGLOB-1)*(IDIM+1)+1)
c           YPT=XCOOR((NGLOB-1)*(IDIM+1)+2)
c           ZPT=0
c           IF (IDIM.EQ.3) ZPT=XCOOR((NGLOB-1)*(IDIM+1)+3)
c           CALL RAFDEN(ICHPO,XPT,YPT,ZPT,IDIM,DENPT)
c           DENSI(NLOC)=DENPT
          DENSI(NLOC)=XDEN(NGLOB)
c          WRITE (*,*) 'densite au pt', NGLOB, XDEN(NGLOB), XPETI
          IF (XDEN(NGLOB).LT.XPETI) THEN
              WRITE (*,*) 'Champ de densité en entrée nul au point',
     +   NGLOB, 'impossible de raffiner'
              RETURN
          ENDIF
 114    CONTINUE
  14    CONTINUE
  15  CONTINUE
C
C=======================================================================
C                Remplissage du tableau KARETE
C=======================================================================
C On remplit le tableau KARETE avec les donnees provenant de l'analyse
C du maillage.
C KARETE(i,j)=k : le ieme noeud local forme une arete
C   avec le keme noeud local (k>i)
C KARPOS(i)=k : compteur de position (k=nb d'aretes touchant i)
C KARET2(i,j)=k : l'arete i-j appartient a k elements
      JPLA1=NBNDS
      JPLA2=NCOL
      SEGINI KARETE,KARPOS,KARET2
      NBARET=0
      NCOL=1
      DO 18 INB=1,MAX(1,IPT5.LISOUS(/1))
        IF (IPT5.LISOUS(/1).EQ.0) THEN
          IPT2=IPT5
        ELSE
          IPT2=IPT5.LISOUS(INB)
          SEGACT IPT2
        ENDIF
c        IF (IPT2.ITYPEL.EQ.48) GOTO 18
        ILPL=LPL(IPT2.ITYPEL)
        ILPT=LPT(IPT2.ITYPEL)
        DO 17 J=1,IPT2.NUM(/2)
        DO 117 K=1,ILPL*2-1,2
          NPTA=IPT2.NUM(KSEGM(ILPT+K-1),J)
          NPTB=IPT2.NUM(KSEGM(ILPT+K),J)
          NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
          NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
          NEXIST=0
          DO 16 I=1,MAX(1,KARPOS(NMIN))
            IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I
  16      CONTINUE
c création d'une nouvelle arrète
          IF (NEXIST.EQ.0) THEN
            IF ((KARPOS(NMIN)+1).GT.NCOL) THEN
              NCOL=KARPOS(NMIN)+1
              SEGADJ KARETE, KARET2
            ENDIF
            KARPOS(NMIN)=KARPOS(NMIN)+1
            KARETE(NMIN,KARPOS(NMIN))=NMAX
            KARET2(NMIN,KARPOS(NMIN))=1
            NBARET=NBARET+1
c            NCOL=MAX(NCOL,KARPOS(NMIN))
          ENDIF
          IF (NEXIST.NE.0) THEN
            KARET2(NMIN,NEXIST)=KARET2(NMIN,NEXIST)+1
          ENDIF
 117    CONTINUE
  17    CONTINUE
  18  CONTINUE
c      SEGADJ KARETE,KARET2
      SEGDES KARETE,KARPOS,KARET2
C
C=======================================================================
C          Remplissage des tableaux JPLANS et JNOEFA
C=======================================================================
C On remplit les tableaux de description des faces dans le cas 3D.
C JPLANS(i,j)=k : la jeme face dont le plus petit sommet en
C                 numerotation locale est le sommet i est la face k
C JNOEFA(i,1-3)=j-k-l : Les 3 plus petits sommets en numerotation
C                       locale de la face i sont les noeuds j,k et l
C JNOEFA(i,4-5)=j-k : La face identifiee par les 3 sommets precedents
C                     appartient aux lisous de type j et k
C JPLCOM(i)=k : compteur de position (k=nb de faces touchant i)
      IF (IDIM.NE.3) GOTO 40
C----------------------------------------------------
C A) Comptage du nombre d'elements total du maillage !
C----------------------------------------------------
      NELTOT=0
      DO 50 INB=1,MAX(1,IPT5.LISOUS(/1))
        IF (IPT5.LISOUS(/1).EQ.0) THEN
          IPT2=IPT5
        ELSE
          IPT2=IPT5.LISOUS(INB)
          SEGACT IPT2
        ENDIF
        NELTOT=NELTOT+IPT2.NUM(/2)
  50  CONTINUE
C
C-------------------------------------------------
C B) Initialisation des tableaux JPLANS et JNOEFA !
C-------------------------------------------------
      JNBFA=NELTOT*6
      JMAXPL=0
      JNUMFA=0
      SEGINI JPLANS,JPLCOM,JNOEFA
C
C-----------------------------------------------
C C)  Remplissage des tableaux JPLANS et JNOEFA !
C-----------------------------------------------
      DO 51 INB=1,MAX(1,IPT5.LISOUS(/1))
        IF (IPT5.LISOUS(/1).EQ.0) THEN
          IPT2=IPT5
        ELSE
          IPT2=IPT5.LISOUS(INB)
c*          SEGACT IPT2
        ENDIF
        JTYP2=IPT2.ITYPEL
C        IF (JTYP2.EQ.48) GOTO 51
        DO 52 KELEM=1,IPT2.NUM(/2)
          JLTEL1=LTEL(1,JTYP2)
          JLTEL2=LTEL(2,JTYP2)
          DO 53 IBB=1,JLTEL1
            JLDEL1=LDEL(1,JLTEL2-1+IBB)
            JTYFAC=IWORK1(JLDEL1)
            JLDEL2=LDEL(2,JLTEL2-1+IBB)
            JNBSOM=NBSOM(JTYFAC)
            JSPOS=NSPOS(JTYFAC)
            SEGINI IWORK2
            DO 54 IAA=1,JNBSOM
              NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),KELEM)
              IWORK2(IAA)=NGLOBA
  54        CONTINUE
C
C----------------------------------------------------
C D)  Tri des sommets de la face par ordre croissant !
C----------------------------------------------------
C On va ranger les 3 plus petits sommets de la face (en numerotation
C globale) dans les variables NPTA, NPTB, NPTC. On les classe par ordre
C croissant.
C
C 1/ Classement des 3 sommets par ordre croissant (num globale)
            NPTA=nbpts+1
            NPTB=NPTA+1
            NPTC=NPTB+1
            DO 55 ICC=1,JNBSOM
              IF (IWORK2(ICC).LT.NPTA) THEN
                NPTC=NPTB
                NPTB=NPTA
                NPTA=IWORK2(ICC)
              ELSEIF (IWORK2(ICC).LT.NPTB) THEN
                NPTC=NPTB
                NPTB=IWORK2(ICC)
              ELSEIF (IWORK2(ICC).LT.NPTC) THEN
                NPTC=IWORK2(ICC)
              ENDIF
  55        CONTINUE
C
C 2/ Classement des 3 sommets par ordre croissant (num locale)
            NPTA2=ICPR(NPTA,1)
            NPTB2=ICPR(NPTB,1)
            NPTC2=ICPR(NPTC,1)
            IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
              NPTA=NPTA2
              NPTB=MIN(NPTB2,NPTC2)
              NPTC=MAX(NPTB2,NPTC2)
            ENDIF
            IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
              NPTA=NPTB2
              NPTB=MIN(NPTA2,NPTC2)
              NPTC=MAX(NPTA2,NPTC2)
            ENDIF
            IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
              NPTA=NPTC2
              NPTB=MIN(NPTA2,NPTB2)
              NPTC=MAX(NPTA2,NPTB2)
            ENDIF
C
C----------------------------------------
C E)  Renseignement des nouvelles faces  !
C----------------------------------------
            SEGACT LTYPE
c            IHERED=0
c            DO 1022 MOI=1,IPT2.NUM(/1)
c              IHERED=IHERED+IRELA(ICPR(IPT2.NUM(MOI,KELEM),1))
c 1022       CONTINUE
 
C Si NPTA Ne correspond a aucune face : On cree une face
            IF (JPLCOM(NPTA).EQ.0) THEN
              JNUMFA=JNUMFA+1
              JPLCOM(NPTA)=JPLCOM(NPTA)+1
              JPLANS(NPTA,JPLCOM(NPTA))=JNUMFA
              JNOEFA(JNUMFA,1)=NPTA
              JNOEFA(JNUMFA,2)=NPTB
              JNOEFA(JNUMFA,3)=NPTC
              JNOEFA(JNUMFA,4)=LTYPE(INB,2)
c              IF (IHERED.NE.0) JNOEFA(JNUMFA,5)=LTYPE(INB,2)
C On signale les sous face de faces supportant un sure
C JNOEFA(JNUMFA,5)= 48
              IF (LCONF3.NE.0) THEN
                SEGACT IPT6
                DO 501 LAA1=1,MAX(1,IPT6.LISOUS(/1))
                  IF (IPT6.LISOUS(/1).EQ.0) THEN
                    IPT4=IPT6
                  ELSE
                    IPT4=IPT6.LISOUS(LAA1)
                  ENDIF
                  SEGACT IPT4
                  IF (IPT4.ITYPEL.NE.48) GOTO 501
                  IF (IPT4.ICOLOR(1).NE.3) GOTO 501
                  DO 502,MAA1=1,IPT4.NUM(/2)
                    NEXI1=0
                    NEXI2=0
                    DO 503 ICC1=1,JNBSOM
                      IF (IWORK2(ICC1).EQ.IPT4.NUM(1, MAA1)) THEN
                        NEXI1=NEXI1+1
                      ENDIF
  503               CONTINUE
                    IF (NEXI1.GT.0) THEN
                      DO 504 ICC2=1,JNBSOM
                      DO 1504 INN2=2,IPT4.NUM(/1)
                        IF (IWORK2(ICC2).EQ.IPT4.NUM(INN2,MAA1)) THEN
                          NEXI2=NEXI2+1
                        ENDIF
 1504                 CONTINUE
  504                 CONTINUE
                      IF (NEXI2.GT.0) JNOEFA(JNUMFA,5)= 48
                    ENDIF
  502             CONTINUE
  501           CONTINUE
              ENDIF
 
 
C Si NPTA  correspond a aumoins une face
            ELSE
              NEXIST=0
              DO 56 IDD=1,JPLCOM(NPTA)
                ITEMP1=JPLANS(NPTA,IDD)
                NAA=JNOEFA(ITEMP1,1)
                NBB=JNOEFA(ITEMP1,2)
                NCC=JNOEFA(ITEMP1,3)
C Si la face que actuelle existe déja : on renseigne JNOEFA(NEXIST,5)
                IF (NAA.EQ.NPTA.AND.NBB.EQ.NPTB.AND.NCC.EQ.NPTC) THEN
                  NEXIST=ITEMP1
                  JNOEFA(NEXIST,5)=LTYPE(INB,2)
                ENDIF
 
 
  56          CONTINUE
C Si la face que actuelle n'existe pas :  On cree une face
              IF (NEXIST.EQ.0) THEN
                JNUMFA=JNUMFA+1
                JPLCOM(NPTA)=JPLCOM(NPTA)+1
                JPLANS(NPTA,JPLCOM(NPTA))=JNUMFA
                JNOEFA(JNUMFA,1)=NPTA
                JNOEFA(JNUMFA,2)=NPTB
                JNOEFA(JNUMFA,3)=NPTC
                JNOEFA(JNUMFA,4)=LTYPE(INB,2)
c                IF (IHERED.NE.0) JNOEFA(JNUMFA,5)=LTYPE(INB,2)
 
C On signale les sous face de faces supportant un sure
C JNOEFA(JNUMFA,5)= 48
                IF (LCONF3.NE.0) THEN
                  SEGACT IPT6
                  DO 505 LAA1=1,MAX(1,IPT6.LISOUS(/1))
                    IF (IPT6.LISOUS(/1).EQ.0) THEN
                      IPT4=IPT6
                    ELSE
                      IPT4=IPT6.LISOUS(LAA1)
                    ENDIF
                    SEGACT IPT4
                    IF (IPT4.ITYPEL.NE.48) GOTO 505
                    IF (IPT4.ICOLOR(1).NE.3) GOTO 505
                    DO 506,MAA1=1,IPT4.NUM(/2)
                      NEXI1=0
                      NEXI2=0
                      DO 507 ICC1=1,JNBSOM
                        IF (IWORK2(ICC1).EQ.IPT4.NUM(1, MAA1)) THEN
                          NEXI1=NEXI1+1
                        ENDIF
  507                 CONTINUE
                      IF (NEXI1.GT.0) THEN
                        DO 508 ICC2=1,JNBSOM
                        DO 1508 INN2=2,IPT4.NUM(/1)
                          IF (IWORK2(ICC2).EQ.IPT4.NUM(INN2,MAA1)) THEN
                            NEXI2=NEXI2+1
                          ENDIF
 1508                   CONTINUE
  508                   CONTINUE
                        IF (NEXI2.GT.0) JNOEFA(JNUMFA,5)= 48
                      ENDIF
  506               CONTINUE
  505             CONTINUE
                ENDIF
              ENDIF
            ENDIF
 
            JMAXPL=MAX(JMAXPL,JPLCOM(NPTA))
  53      CONTINUE
  52    CONTINUE
  51  CONTINUE
      JPLA2=JMAXPL
      JNBFA=JNUMFA
      SEGADJ JPLANS,JPLCOM,JNOEFA
      SEGDES LTYPE
C
C=======================================================================
C                  Pre-remplissage du tableau KMILIE
C=======================================================================
C On remplit le tableau lorsqu'il y a des maillages de relations de
C type 48 (IPT6, LCONF3 non nul)
C KARETE(i,j)=k ET KMILIE(i,j)=n : le noeud global n est situe
C au milieu de l'arete [i,k]
  40  CONTINUE
      SEGINI KMILIE
      SEGACT KARETE*MOD,KARPOS*MOD,KARET2*MOD
      IF (LCONF3.EQ.0) THEN
        GOTO 28
      ENDIF
  22  CONTINUE
      NBRELA=0
      SEGACT IPT6
      DO 218 LAA=1,MAX(1,IPT6.LISOUS(/1))
        IF (IPT6.LISOUS(/1).EQ.0) THEN
          IPT2=IPT6
        ELSE
          IPT2=IPT6.LISOUS(LAA)
        ENDIF
        SEGACT IPT2
        IF (IPT2.ITYPEL.NE.48) GOTO 218
        IF ((IPT2.ICOLOR(1).NE.1).AND.(IPT2.ICOLOR(1).NE.2)) GOTO 218
        DO 220 J=1,IPT2.NUM(/2)
 
          NPTA=IPT2.NUM(2,J)
          NPTB=IPT2.NUM(3,J)
          NPTC=IPT2.NUM(1,J)
          NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
          NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
          NEXIST=0
          DO 219 I=1,MAX(1,KARPOS(NMIN))
            IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I
 219      CONTINUE
 
          IF (NEXIST.NE.0) then
              KMILIE(NMIN,NEXIST)=NPTC
c              WRITE (*,*) 'REPORT ANCIENNE RELA AU NOEUD:', NPTC
c              write (*,*) 'karret2(nmin nexist)', KARET2(NMIN,NEXIST)
              KARET2(NMIN,NEXIST)=2
          ELSE
              IF (IDIM.EQ.2) goto 27
              IF ((KARPOS(NMIN)+1).GT.NCOL) THEN
                NCOL=KARPOS(NMIN)+1
                SEGADJ KARETE, KARET2, KMILIE
              ENDIF
              JPOSI = KARPOS(NMIN)
              KARPOS(NMIN)=JPOSI+1
              KARETE(NMIN,JPOSI+1)=NMAX
              KARET2(NMIN,JPOSI+1)=2
              KMILIE(NMIN,JPOSI+1)=NPTC
c              WRITE (*,*) 'création de l arrete:', NMIN, NMAX
          ENDIF
C
C-----------------------------------------------
C Prise en compte des relations "hereditaires" !
C-----------------------------------------------
C Cette partie sert a traiter le cas ou un element serait raffine
C plusieurs fois alors qu'a sa frontiere les autres elements ne le
C seraient pas. Il y aurait donc plusieurs points sur une arete ayant
C des relations de conformite et plus seulement le noeud milieu.
C
c dans l'exenple ci dessous ou les noeuds C et D sont des hanging nodes
c L'arete [C B] n'a pas encore ete creee c'est ce qui est fait ici.
c
c |                   |
c |                   |
c |A        C    D    |B
c x---------X----X----X
c |         |    |    |
c |         |    |    |
   27     CONTINUE
          IF (IDIM.EQ.2) THEN
C 1/ Test sur la premiere arete : NMIN-NPTC
            IF (NEXIST.NE.0) THEN
              NEXIS2=0
              NMIN2=MIN(NMIN,ICPR(NPTC,1))
              NMAX2=MAX(NMIN,ICPR(NPTC,1))
              DO 24 I=1,MAX(1,KARPOS(NMIN2))
                IF (KARETE(NMIN2,I).EQ.NMAX2) NEXIS2=I
   24         CONTINUE
              IF (NEXIS2.EQ.0) THEN
                NBRELA=NBRELA+1
                JNBCOL=KARETE(/2)
                JPOSI=KARPOS(NMIN2)
                IF (JNBCOL.GT.JPOSI) THEN
                  KARPOS(NMIN2)=JPOSI+1
                  KARETE(NMIN2,JPOSI+1)=NMAX2
                  KARET2(NMIN2,JPOSI+1)=2
                ELSEIF (JNBCOL.EQ.JPOSI) THEN
                  NCOL=JNBCOL+1
                  SEGADJ KARETE,KMILIE,KARET2
                  KARETE(NMIN2,JPOSI+1)=NMAX2
                  KARPOS(NMIN2)=JPOSI+1
                  KARET2(NMIN2,JPOSI+1)=2
                ENDIF
c              WRITE (*,*) 'creation de l arrete:', NMIN2, NMAX2
              ELSE
                KARET2(NMIN2,NEXIS2)=2
              ENDIF
C
C 2/ Test sur la deuxieme arete : NMAX-NPTC
              NEXIS2=0
              NMIN2=MIN(NMAX,ICPR(NPTC,1))
              NMAX2=MAX(NMAX,ICPR(NPTC,1))
              DO 25 I=1,MAX(1,KARPOS(NMIN2))
                IF (KARETE(NMIN2,I).EQ.NMAX2) NEXIS2=I
  25          CONTINUE
              IF (NEXIS2.EQ.0) THEN
                NBRELA=NBRELA+1
                JNBCOL=KARETE(/2)
                JPOSI=KARPOS(NMIN2)
                IF (JNBCOL.GT.JPOSI) THEN
                  KARPOS(NMIN2)=JPOSI+1
                  KARETE(NMIN2,JPOSI+1)=NMAX2
                  KARET2(NMIN2,JPOSI+1)=2
                ELSEIF (JNBCOL.EQ.JPOSI) THEN
                  NCOL=JNBCOL+1
                  SEGADJ KARETE,KMILIE,KARET2
                  KARETE(NMIN2,JPOSI+1)=NMAX2
                  KARPOS(NMIN2)=JPOSI+1
                  KARET2(NMIN2,JPOSI+1)=2
                ENDIF
c              WRITE(*,*) 'creation de l arrete:', NMIN2, NMAX2
              ELSE
                KARET2(NMIN2,NEXIS2)=2
              ENDIF
            ENDIF
          ENDIF
 220    CONTINUE
 218  CONTINUE
      IF (NBRELA.NE.0) GOTO 22
  28  CONTINUE
 
C=======================================================================
C                  Pre-remplissage du tableau JFARAF
C=======================================================================
C On remplit le tableau JFARAF et JPLAN3 lorsqu'il y a des maillages de
C relations de type 48 (IPT6, LCONF3 non nul).
c
C Si JPLANS(i,j)=k et JFARAF(k , (l-1)*2+1)= n  : le noeud n (en
C numérotation globale) est le hanging node au centre de la jieme face
C dont le plus petit sommet en numerotation locale est le sommet i.
C
C Si de plus la face en question est un TRI6 et JFARAF(k , l*2)=m :
C le noeud m (en numérotation globale) est le 4ieme noeud de l'element
C 48 coul 4
C
C Si de plus la face en question est un QUA8-1 et JFARAF(k , l*2)=m :
C le noeud m (en numérotation globale) est le 5eme noeud de l'element
c 48 coul 4
C
C si JPLAN3(i,j)=k il y a k relation de compatibilité de faces associé a
C la jieme face dont le plus petit sommet en numerotation locale est le
C sommet i.
 
      IF (IDIM.NE.3) GOTO 711
      JNBFA=JNOEFA(/1)
      LLLL=10
      SEGINI JFARAF
      SEGINI JPLAN2,JPLAN3
      IF (LCONF3.EQ.0) GOTO 711
      MPOS25=0
      JNBSOM=0
 1287 CONTINUE
      NBRELA=0
      DO 670 MIE=1,MAX(1,IPT6.LISOUS(/1))
        IF (IPT6.LISOUS(/1).EQ.0) THEN
          IPT2=IPT6
        ELSE
          IPT2=IPT6.LISOUS(MIE)
        ENDIF
        SEGACT IPT2
        IF (IPT2.ITYPEL.NE.48) GOTO 670
        IF ((IPT2.ICOLOR(1).EQ.1).OR.(IPT2.ICOLOR(1).EQ.2)) GOTO 670
        LIGNUM=IPT2.NUM(/1)
        IF (IPT2.ICOLOR(1).EQ.4) MPOS25=4
        IF (IPT2.ICOLOR(1).EQ.5) MPOS25=5
        IF (IPT2.ICOLOR(1).EQ.3) JNBSOM=4
        IF (IPT2.ICOLOR(1).EQ.4) JNBSOM=3
        IF (IPT2.ICOLOR(1).EQ.5) JNBSOM=4
        IF (IPT2.ICOLOR(1).EQ.6) JNBSOM=4
        SEGINI IWORK2
C
C-------------------------
C A) Recherche de la face !
C-------------------------
C 1/ Classement des 3 sommets par ordre croissant (num globale)
        DO 671,MAA=1,IPT2.NUM(/2)
          DO 673 MBB=1,JNBSOM
            IWORK2(MBB)=IPT2.NUM(LIGNUM-JNBSOM+MBB,MAA)
 673      CONTINUE
          NPTA=nbpts+1
          NPTB=NPTA+1
          NPTC=NPTB+1
          DO 672 ICC=1,JNBSOM
            IF (IWORK2(ICC).LT.NPTA) THEN
              NPTC=NPTB
              NPTB=NPTA
              NPTA=IWORK2(ICC)
            ELSEIF (IWORK2(ICC).LT.NPTB) THEN
              NPTC=NPTB
              NPTB=IWORK2(ICC)
            ELSEIF (IWORK2(ICC).LT.NPTC) THEN
              NPTC=IWORK2(ICC)
            ENDIF
 672      CONTINUE
C
C 2/ Classement des 3 sommets par ordre croissant (num locale)
          NPTA1=NPTA
          NPTB1=NPTB
          NPTC1=NPTC
          NPTA2=ICPR(NPTA,1)
          NPTB2=ICPR(NPTB,1)
          NPTC2=ICPR(NPTC,1)
          IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
            NPTA=NPTA2
            NPTB=MIN(NPTB2,NPTC2)
            NPTC=MAX(NPTB2,NPTC2)
          ENDIF
          IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
            NPTA=NPTB2
            NPTB=MIN(NPTA2,NPTC2)
            NPTC=MAX(NPTA2,NPTC2)
          ENDIF
          IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
            NPTA=NPTC2
            NPTB=MIN(NPTA2,NPTB2)
            NPTC=MAX(NPTA2,NPTB2)
          ENDIF
C
C 3/ Recherche du numero de la face
          NEXIS2=0
          DO 674 IEE=1,JPLCOM(NPTA)
            MTMP=JPLANS(NPTA,IEE)
            JJ1=JNOEFA(MTMP,1)
            JJ2=JNOEFA(MTMP,2)
            JJ3=JNOEFA(MTMP,3)
            IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
              NEXIS2=IEE
            ENDIF
 674      CONTINUE
c 3.1/ si la face en question n'existe pas, dans le cas général on
c      s'arrète
          IF ((NEXIS2.EQ.0).AND.(IPT2.ICOLOR(1).NE.3)) THEN
c            WRITE (*,*) 'IPT2.ICOLOR(1)', IPT2.ICOLOR(1)
            WRITE (*,*) 'Raffinement impossible, le champ de densite
     + presente des variations trop brusques'
            WRITE (*,*) 'sure concerne', MAA
            WRITE (*,*) 'premier nd du sure concerne', IPT2.NUM(1,MAA)
            WRITE (*,*) 'numero du 1er nd de la face concernee', NPTA
            KARAF2=0
c            GOTO 679
          ENDIF
c 3.2/ si la face en question n'existe pas dans le cas d'une relation
c de face QUA4 il suffie de la créer
 
          IF ((NEXIS2.EQ.0).AND.(IPT2.ICOLOR(1).EQ.3)) THEN
            JPLCOM(NPTA)= JPLCOM(NPTA)+1
            JNUMFA=JNUMFA+1
            JNUMFS=JNUMFA
            JNUMF1=JPLCOM(NPTA)
            JMAXPL=MAX(JMAXPL,JPLCOM(NPTA))
c 3.2.1/ ajustement des segments
            JPLA2=JMAXPL
            JNBFA=JNUMFA
            SEGADJ JPLANS
            SEGADJ JNOEFA
            SEGADJ JPLAN3
            SEGADJ JFARAF
C 3.2.1 Renseignement des iformation de la nouvelle face dans JPLANS et
c JNOEFA, initialisation de  JPLAN3
            JPLANS(NPTA,JNUMF1)=JNUMFA
            JPLAN3(NPTA,JNUMF1)=0
            JNOEFA(JNUMFA,1)= NPTA
            JNOEFA(JNUMFA,2)= NPTB
            JNOEFA(JNUMFA,3)= NPTC
            JNOEFA(JNUMFA,4)= 48
            JNOEFA(JNUMFA,5)= 48
          ENDIF
 
          IF (NEXIS2.NE.0) THEN
            JNUMFS=JPLANS(NPTA,NEXIS2)
          ENDIF
C
C 4/ Recherche de la position de la face dans JPLANS
          NEXIS2=0
          DO 675 MAC=1,JPLCOM(NPTA)
            IF (JPLANS(NPTA,MAC).EQ.JNUMFS) NEXIS2=MAC
 675      CONTINUE
C
C-------------------------------------
C B) Renseignement de JFARAF et JPLAN3!
C-------------------------------------
 
          JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1
          LTEMPO=(JPLAN3(NPTA,NEXIS2)-1)*2+1
          JFARAF(JNUMFS,LTEMPO)=IPT2.NUM(1,MAA)
          IF (MPOS25.NE.0) JFARAF(JNUMFS,LTEMPO+1)=IPT2.NUM(MPOS25,MAA)
          IF (JNOEFA(JNUMFS,5).EQ.0) JNOEFA(JNUMFS,5)=48
 
 
 671    CONTINUE
 670  CONTINUE
C      IF (NBRELA.NE.0) GOTO 1287
 711  CONTINUE
 
C=======================================================================
C           Calcul du critere et raffinement si necessaire
C=======================================================================
C Pour calculer le critere, on a besoin d'un modele et d'un chamelem
C constant egal a 1.
C
C------------------------------------------
C A) Creation d'un modele simple : MMODE1 !
C------------------------------------------
      MN3 =12
      NFOR=1
      NMAT=2
      NPARMO=0
      NOBMOD=0
      N1=1
      SEGINI MMODE1
      SEGINI IMODE1
      MMODE1.KMODEL(1)=IMODE1
      IMODE1.FORMOD(1)='MECANIQUE'
      IMODE1.MATMOD(1)='ELASTIQUE'
      IMODE1.MATMOD(2)='ISOTROPE'
      IMODE1.CONMOD   ='                        '
 
      KARAF2=0
C=======================================================================
C           Boucle sur les eventuels sous-domaines de IPT5
C=======================================================================
 
c-----------------------------------------------------------------------
c      initialisation de LTYPE3 qui va contenir les tableau LTYP2
c      correspondant à chaque sous domaine de IPT5. Ces tableau LTYP2
c      contiendront eux même les maillages raffinés issue des
c      différents sous domaine de IPT5.
c      LTYPE3(I,1)=LTYP2 pour le sous domaine I de IPT5
c      LTYPE3(I,2)=IRR le permier indice d'un maillage 48 dans LTYP2
c-----------------------------------------------------------------------
      NBSOU5=IPT5.LISOUS(/1)
      NBLIS3=MAX(1,NBSOU5)
      SEGINI LTYPE3
      MCONF=0
      NSUR = 0
      DO 33 INB=1,NBLIS3
        IF (NBSOU5.EQ.0) THEN
          IPT2=IPT5
        ELSE
          IPT2=IPT5.LISOUS(INB)
        ENDIF
 
        IMODE1.IMAMOD=IPT2
        IMODE1.NEFMOD=IPT2.ITYPEL
        CALL PRQUOI(IMODE1)
c PRQUOI n'active que le dernier MINTE de INFMOD qui n'est pas forcement
c celui utilise par INTGCA. Appel a ACTOBJ pour activer tous les MINTE.
        CALL ACTOBJ('MMODEL  ',MMODE1,1)
 
C-------------------------------------------------------
C B) Creation d'un chamelem constant egal a 1 : MCHEL1 !
C-------------------------------------------------------
        CALL ZEROP(MMODE1,'NOEUD   ',MCHEL1)
c*        SEGACT MCHEL1
        MCHAM1=MCHEL1.ICHAML(1)
        SEGACT MCHAM1
        MELVA1=MCHAM1.IELVAL(1)
        SEGACT MELVA1*MOD
        MELVA1.VELCHE(1,1)=1.D0
C
C-------------------------------------------------------------
C C) Calcul de la surface ou du volume des elements : MCHEL2 !
C-------------------------------------------------------------
        ichm2 = 0
        CALL INTGCA(MMODE1,MCHEL1,ichm2,1,IRET,XRET,MCHEL2)
        IF (IRET.NE.1) RETURN
        IF (XRET.LT.XPETI) THEN
            WRITE (*,*) 'Volume(resp surf) de la zone a raffiner nul'
            RETURN
        ENDIF
c*        SEGACT IPT2
 
C---------------------------------------------------------------------
C D) Creation d'un chamelem indiquant les elements a raffiner MCHAM2 !
C---------------------------------------------------------------------
        SEGACT KMILIE*MOD
        SEGACT MCHEL2
        MCHAM2=MCHEL2.ICHAML(1)
        SEGACT MCHAM2
        MELVA2=MCHAM2.IELVAL(1)
        SEGACT MELVA2*MOD
        NACREE=0
        KARAF=0
 
        XDIM=1.D0/IDIM
 
        NBNN2=IPT2.NUM(/1)
        JTYP2=IPT2.ITYPEL
 
c       boucles sur les éléments de ipt2
        DO 32 J=1,IPT2.NUM(/2)
 
          XMOY=0.D0
          DO 29 I=1,NBNN2
            NGLOB=IPT2.NUM(I,J)
            NLOC=ICPR(NGLOB,1)
            DENPT=DENSI(NLOC)
            XMOY=XMOY+DENPT/NBNN2
  29      CONTINUE
 
          XINT=MELVA2.VELCHE(1,J)
C Un element est a raffiner si la moyenne des densites affectees a
C l'ensemble de ses noeuds est inférieure a la racine carree de sa
C surface (en 2D) ou a la racine cubique de son volume (en 3D)
C          WRITE (*,*) 'XMOY', XMOY, '(XINT**XDIM)',(XINT**XDIM)
          IF ((XMOY-(XINT**XDIM)).LT.0) THEN
c c         y a t'il des +epsilon...?
c           XDIFF = XMOY-(XINT**XDIM)
c           if(XDIFF.ge.0.and.XDIFF.lt.0.001*XMOY) then
c             write(*,*) 'RAFF : XMOY,XDIFF=',XMOY,XDIFF,
c      &     'pour l EF',J,(IPT2.NUM(iou,J),iou=1,IPT2.NUM(/1))
c           endif
c           IF (XDIFF.LT.0) THEN
            MELVA2.VELCHE(1,J)=1
c KARAF nombre d'éléments à rafiner dans le sous domaine
            KARAF=KARAF+1
            ILPT=LPT(JTYP2)
            ILPL=LPL(JTYP2)
 
c boucle sur les segments de l'élément à raffiner
            DO 31 K=1,ILPL*2-1,2
              NPTA=IPT2.NUM(KSEGM(ILPT+K-1),J)
              NPTB=IPT2.NUM(KSEGM(ILPT+K),J)
              NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
              NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
              NEXIST=0
              DO 30 I=1,MAX(1,KARPOS(NMIN))
                IF (KARETE(NMIN,I).EQ.NMAX) NEXIST=I
  30          CONTINUE
c On met dans kmilie les arêtes où il faudra mettre un nouveau noeud
c On leur associe la valeur -1
              IF (KMILIE(NMIN,NEXIST).EQ.0) THEN
                NACREE=NACREE+1
                KMILIE(NMIN,NEXIST)=-1
              ENDIF
  31        CONTINUE
C
C-------------------------------------------
C E) Comptage du nombre de faces a raffiner !
C-------------------------------------------
C 1/ Initialisations
            IF (IDIM.NE.3) GOTO 32
            JLTEL1=LTEL(1,JTYP2)
            JLTEL2=LTEL(2,JTYP2)
            DO 73 IBB=1,JLTEL1
              JLDEL1=LDEL(1,JLTEL2-1+IBB)
              JTYFAC=IWORK1(JLDEL1)
              JLDEL2=LDEL(2,JLTEL2-1+IBB)
C             JNBPTS=NBNNE(JTYFAC)
              JNBSOM=NBSOM(JTYFAC)
              JSPOS=NSPOS(JTYFAC)
              SEGINI IWORK2
              DO 74 IAA=1,JNBSOM
                NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),J)
                IWORK2(IAA)=NGLOBA
  74          CONTINUE
C
C 2/ Classement des 3 sommets par ordre croissant (num globale)
              NPTA=nbpts+1
              NPTB=NPTA+1
              NPTC=NPTB+1
              DO 75 ICC=1,JNBSOM
                IF (IWORK2(ICC).LT.NPTA) THEN
                  NPTC=NPTB
                  NPTB=NPTA
                  NPTA=IWORK2(ICC)
                ELSEIF (IWORK2(ICC).LT.NPTB) THEN
                  NPTC=NPTB
                  NPTB=IWORK2(ICC)
                ELSEIF (IWORK2(ICC).LT.NPTC) THEN
                  NPTC=IWORK2(ICC)
                ENDIF
  75          CONTINUE
C
C 3/ Classement des 3 sommets par ordre croissant (num locale)
              NPTA2=ICPR(NPTA,1)
              NPTB2=ICPR(NPTB,1)
              NPTC2=ICPR(NPTC,1)
              IF ((NPTA2.LT.NPTB2).AND.(NPTA2.LT.NPTC2)) THEN
                NPTA=NPTA2
                NPTB=MIN(NPTB2,NPTC2)
                NPTC=MAX(NPTB2,NPTC2)
              ENDIF
              IF ((NPTB2.LT.NPTA2).AND.(NPTB2.LT.NPTC2)) THEN
                NPTA=NPTB2
                NPTB=MIN(NPTA2,NPTC2)
                NPTC=MAX(NPTA2,NPTC2)
              ENDIF
              IF ((NPTC2.LT.NPTA2).AND.(NPTC2.LT.NPTB2)) THEN
                NPTA=NPTC2
                NPTB=MIN(NPTA2,NPTB2)
                NPTC=MAX(NPTA2,NPTB2)
              ENDIF
C
C 4/ Calcul du nombre de points a creer pour raffiner la face
              NEXIS2=0
              DO 76 IEE=1,JPLCOM(NPTA)
                MTMP=JPLANS(NPTA,IEE)
                JJ1=JNOEFA(MTMP,1)
                JJ2=JNOEFA(MTMP,2)
                JJ3=JNOEFA(MTMP,3)
                IF(JJ1.EQ.NPTA.AND.JJ2.EQ.NPTB.AND.JJ3.EQ.NPTC) THEN
                  NEXIS2=IEE
                ENDIF
  76          CONTINUE
              IF (NEXIS2.NE.0) THEN
                KFARAF=JPLAN2(NPTA,NEXIS2)
                IF (KFARAF.EQ.0) THEN
                  NACREE=NACREE+NBINTE(JTYFAC)
                  JPLAN2(NPTA,NEXIS2)=JPLAN2(NPTA,NEXIS2)+1
                ENDIF
              ENDIF
  73        CONTINUE
          ELSE
 
            MELVA2.VELCHE(1,J)=0
          ENDIF
  32    CONTINUE
 
c KARAF2 nombre d'éléments à rafiner sur tout les sous domaines
        KARAF2=KARAF2+KARAF
c si aucun élément à raffiner on stoque directement IPT2 dans LTYPE3
c on stoque également les anciennes relations et on passe au
c sous-domaine suivant.
 
c ATTENTION SI NBOUCL=1 il est possible que des relations surabondantes
 
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c   forcage  : nb max itération de raff
c       IF (NBOUCL.EQ.3) THEN
c           KARAF2=0
c           KARAF=0
c       ENDIF
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        IF (KARAF.EQ.0) THEN
            IF (NBOUCL.EQ.1) THEN
                IF (LCONF3.GT.0) THEN
                    SEGACT IPT6
                    NBLIS2=1+MAX(1,IPT6.LISOUS(/1))
                    SEGINI LTYPE2
                    LTYPE2(1,1)=1
                    LTYPE2(1,2)=IPT2.ITYPEL
                    LTYPE2(1,3)=IPT2
                    LTYPE2(1,4)=IPT2.ICOLOR(1)
                    LTYPE3(INB,2)=0
                    DO 676 MIE=1,MAX(1,IPT6.LISOUS(/1))
                        IF (IPT6.LISOUS(/1).EQ.0) THEN
                            IPT4=IPT6
                        ELSE
                            IPT4=IPT6.LISOUS(MIE)
                        ENDIF
                        NSUR=1
                       SEGACT IPT4
                       LTYPE2(1+MIE,1)=1+MIE
                       LTYPE2(1+MIE,2)=IPT4.ITYPEL
                       LTYPE2(1+MIE,3)=IPT4
                       LTYPE2(1+MIE,4)=IPT4.ICOLOR(1)
                       LTYPE3(INB,2)=2
  676               CONTINUE
                ELSE
                    NBLIS2=1
                    SEGINI LTYPE2
                    LTYPE2(1,1)=1
                    LTYPE2(1,2)=IPT2.ITYPEL
                    LTYPE2(1,3)=IPT2
                    LTYPE2(1,4)=IPT2.ICOLOR(1)
                    LTYPE3(INB,2)=0
                ENDIF
                LTYPE3(INB,1)=LTYPE2
            ELSE
                SEGACT LTYPE4
                LTYPE3(INB,1)=LTYPE4(INB,1)
                LTYPE3(INB,2)=LTYPE4(INB,2)
            ENDIF
            GOTO 33
        ENDIF
C
C----------------------------
C F) Raffinement du maillage !
C----------------------------
        NML=2
C
C 1/ Raffinement d'un maillage 2D
C On recupere dans MELEME le maillage raffine sans ses relations de
C conformite et dans IPT8 les relations de ce maillage
C (IPT8 encore incomplet en sortie de RAFF2D)
c entrees de raff2D:
c           - IPT2:   maillage élémentaire à raffiner
c                     (ici sous domaine de ipt5)
c           - ICPR:   tableau de passage noeuds local/global
c                     (ici à partir de ipt5 interne à la boucle 9)
c           - KARPOS: tableau du nb d'arretes par noeuds
c                     (ici à partir de ipt5, mis a jour avec ipt6)
c           - KARETE: tableau des d'arretes KARETE(i,j)=k
c                     (ici à partir de ipt5, mis a jour avec ipt6)
c           - KMILIE: tableau des hanging nodes associés au arrètes
c                     en entrée: à partir des relations de ipt6 et
c                     des nouveau noeuds à construires
c                     si relation issue de ipt6
c                        KMILIE(i,j)=n : le noeud global n est situe
c                        au milieu de l'arete [i,k]
c                     si nouvelle relation
c                        KMILIE(i,j)=-1 : il faut construire un noeud
c                        au milieu de l'arete [i,k]
c           - MELVA2: champ par elem valait 1 pour les élem à raffiner
c                     et 0 pour les autres
c                     (ici à partir de ipt2 et de densi)
c           - NACREE: nombre de noeuds à créer dans ipt2
c           - KARAF:  nombre d'élément à raffiner dans ipt2
c           - IDCP :  ??? variable non utilisé ???
c           - XDEN :  densité aux noeuds en notation globale
c                     (ici à partir de mchpo1)
c           - KARET2: tableau du nombre d'éléments touchant les arretes
c                     (ici à partir de ipt5, mis a jour avec ipt6)
csorties de raff2D :
c           - KMILIE: tableau des hanging nodes associés au arrètes
c                     en sortie: à partir des relations donnée en
c                     entrée et de celles crées
c                     si anciene relation à garder
c                       KMILIE(i,j)=n : le noeud global n est situe
c                       au milieu de l'arete [i,k]
c                     si anciene relation à suprimer
c                       KMILIE(i,j)=0 :
c                     si nouvelle relation
c                       KMILIE(i,j)=n : le noeud n à été construit au
c                       milieu de l'arete [i,k]
c           - MELEME: Maillage élémentaire raffiné à partir de ipt2
c                     sans relations
c           - IPT8 :  Maillage de relation incomplet
c                     si J nouvel élément 48:
c                       IPT8.NUM(1,j) = nouveau hanging node
c                       IPT8.NUM(2:4,j) = autres noeuds de la relation
c                       IPT8.NUM(5,j) = 1 si arête a 2 noeuds
c                                       2 si arête a 3 noeuds
c                     si J nouvel élément 48:
c                       IPT8.NUM(1,j) = nouveau hanging node
c                       IPT8.NUM(2:5,j) = 0
c           - XDEN :  densité aux noeuds en notation globale
c                     interpolé aux nouveaux noeuds.
        IF (IDIM.EQ.2) CALL RAFF2D(IPT2,ICPR,KARPOS,KARETE,KMILIE,
     +   MELVA2,NACREE,KARAF,MELEME,IDCP,IPT8,KARET2,XDEN)
C
C 2/ Raffinement d'un maillage 3D
        IF (IDIM.EQ.3) THEN
          NML=0
          IF ((IPT2.ITYPEL.EQ.25).OR.(IPT2.ITYPEL.EQ.26)) NML=1
        ENDIF
C   a) Maillage 3D contenant des pyramides
        IPT9=0
        IF (NML.EQ.1) CALL RAFPYR(IPT2,ICPR,KARPOS,KARETE,KMILIE,
     +    MELVA2,NACREE,KARAF,MELEME,JPLANS,JPLAN3,JPLCOM,JNOEFA,
     +    IPT8,JFARAF,KARET2,IPT9,XDEN)
C   b) Maillage 3D sans pyramides
        IF (NML.EQ.0) CALL RAFF3D(IPT2,IPT6,ICPR,KARPOS,KARETE,KMILIE,
     +    MELVA2,NACREE,KARAF,MELEME,JPLANS,JPLAN3,JPLCOM,JNOEFA,
     +    IPT8,JFARAF,KARET2,XDEN)
C
C=======================================================================
C          Creation du maillage de relations de conformite
C=======================================================================
C--------------------
C A) Initialisations !
C--------------------
      IPT10=0
      IPT11=0
      IPT12=0
      IPT13=0
      IPT14=0
      IPT15=0
      IPT16=0
      MCONF=0
c      IF (KARAF2.EQ.0) GOTO 33
      IF (IPT8.NE.0) THEN
C
C---------------------------------------
C B) Traitement des anciennes relations !
C---------------------------------------
C A ce stade, IPT8 contient les relations crees lors du dernier
C raffinement (noeud de relation + noeuds formant la relation) et les
C anciennes relations a reporter (noeud de relation seulement)
 
C 1/ Recherche de l'existence d'anciennes relations a reporter
        NB10=0
        MLONG=IPT8.NUM(/1)
        DO 775 MM=1,IPT8.NUM(/2)
          IF (IPT8.NUM(MLONG,MM).EQ.0) NB10=NB10+1
 775    CONTINUE
 
        IF (NB10.NE.0) THEN
C
C 2/ Creation de IPT10
c          MCOL=0
c          NBNN=9
c          NBELEM=NB10
c          NBREF=0
c          NBSOUS=0
c          SEGINI IPT10
c          IPT10.ITYPEL=48
C
C 3/ Recherche des anciennes relations a reporter
C Cette partie permet de compléter IPT8 avec les noeuds formant les
C relations, pour les anciennes relations à reporter
          DO 628 MM11=1,IPT8.NUM(/2)
            IF (IPT8.NUM(MLONG,MM11).NE.0) GOTO 628
            INRELA=IPT8.NUM(1,MM11)
c On active pous le cas ou IPT6=IPT2.
            SEGACT IPT6
            DO 625 KOU=1,MAX(1,IPT6.LISOUS(/1))
              IF (IPT6.LISOUS(/1).EQ.0) THEN
                IPT2=IPT6
              ELSE
                IPT2=IPT6.LISOUS(KOU)
              ENDIF
              SEGACT IPT2
              IF (IPT2.ITYPEL.NE.48) GOTO 625
              DO 626 JKG=1,IPT2.NUM(/2)
                INOEGL=IPT2.NUM(1,JKG)
                IF (INOEGL.NE.INRELA) GOTO 626
                MCOL=MCOL+1
                DO 627 MTMP=1,IPT2.NUM(/1)
                  IPT8.NUM(MTMP,MM11)=IPT2.NUM(MTMP,JKG)
                  IPT8.NUM(MLONG,MM11)= IPT2.ICOLOR(JKG)
 627            CONTINUE
                GOTO 628
 626          CONTINUE
 625        CONTINUE
 628      CONTINUE
        ENDIF
c à ce stade le maillage de relation ipt8 est complet
C
C---------------------------------------
C C) Traitement des nouvelles relations !
C---------------------------------------
C Ici, IPT8 contient toutes les relations du nouveau maillage
C On les trie par type de relations
 
C 1/ Recherche de l'existence de nouvelles relations
      NB11=0
      NB12=0
      NB13=0
      NB14=0
      NB15=0
      NB16=0
      MLONG=IPT8.NUM(/1)
      DO 776 MM=1,IPT8.NUM(/2)
        IF (IPT8.NUM(MLONG,MM).EQ.1) NB11=NB11+1
        IF (IPT8.NUM(MLONG,MM).EQ.2) NB12=NB12+1
        IF (IPT8.NUM(MLONG,MM).EQ.3) NB13=NB13+1
        IF (IPT8.NUM(MLONG,MM).EQ.4) NB14=NB14+1
        IF (IPT8.NUM(MLONG,MM).EQ.5) NB15=NB15+1
        IF (IPT8.NUM(MLONG,MM).EQ.6) NB16=NB16+1
 776  CONTINUE
C
C 2/ Cas des relations 'aretes a 2 noeuds'
        IF (NB11.NE.0) THEN
          MCONF=MCONF+1
          MCOL=0
          NBNN=3
          NBELEM=NB11
          NBREF=0
          NBSOUS=0
          SEGINI IPT11
          IPT11.ITYPEL=48
          DO 778 MM11=1,IPT8.NUM(/2)
            IF (IPT8.NUM(MLONG,MM11).NE.1) GOTO 778
            MCOL=MCOL+1
            DO 777 MTMP=1,3
              IPT11.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
              IPT11.ICOLOR(MCOL)=1
 777        CONTINUE
 778      CONTINUE
        ENDIF
C
C 3/ Cas des relations 'aretes a 3 noeuds'
        IF (NB12.NE.0) THEN
          MCONF=MCONF+1
          MCOL=0
          NBNN=4
          NBELEM=NB12
          NBREF=0
          NBSOUS=0
          SEGINI IPT12
          IPT12.ITYPEL=48
          DO 668 MM11=1,IPT8.NUM(/2)
            IF (IPT8.NUM(MLONG,MM11).NE.2) GOTO 668
            MCOL=MCOL+1
            DO 667 MTMP=1,4
              IPT12.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
              IPT12.ICOLOR(MCOL)=2
 667        CONTINUE
 668      CONTINUE
        ENDIF
C
C 4/ Cas des relations 'face QUA4'
        IF (NB13.NE.0) THEN
          MCONF=MCONF+1
          MCOL=0
          NBNN=5
          NBELEM=NB13
          NBREF=0
          NBSOUS=0
          SEGINI IPT13
          IPT13.ITYPEL=48
          DO 658 MM11=1,IPT8.NUM(/2)
            IF (IPT8.NUM(MLONG,MM11).NE.3) GOTO 658
            MCOL=MCOL+1
            DO 657 MTMP=1,5
              IPT13.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
              IPT13.ICOLOR(MCOL)=3
 657        CONTINUE
 658      CONTINUE
        ENDIF
C
C 5/ Cas des relations 'face TRI6'
        IF (NB14.NE.0) THEN
          MCONF=MCONF+1
          MCOL=0
          NBNN=7
          NBELEM=NB14
          NBREF=0
          NBSOUS=0
          SEGINI IPT14
          IPT14.ITYPEL=48
          DO 648 MM11=1,IPT8.NUM(/2)
            IF (IPT8.NUM(MLONG,MM11).NE.4) GOTO 648
            MCOL=MCOL+1
            DO 647 MTMP=1,7
              IPT14.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
              IPT14.ICOLOR(MCOL)=4
 647        CONTINUE
 648      CONTINUE
        ENDIF
C
C 6/ Cas des relations 'face QUA8-1'
        IF (NB15.NE.0) THEN
          MCONF=MCONF+1
          MCOL=0
          NBNN=9
          NBELEM=NB15
          NBREF=0
          NBSOUS=0
          SEGINI IPT15
          IPT15.ITYPEL=48
          DO 638 MM11=1,IPT8.NUM(/2)
            IF (IPT8.NUM(MLONG,MM11).NE.5) GOTO 638
            MCOL=MCOL+1
            DO 637 MTMP=1,9
              IPT15.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
              IPT15.ICOLOR(MCOL)=5
 637        CONTINUE
 638      CONTINUE
        ENDIF
C
C 7/ Cas des relations 'face QUA8-2'
        IF (NB16.NE.0) THEN
          MCONF=MCONF+1
          MCOL=0
          NBNN=9
          NBELEM=NB16
          NBREF=0
          NBSOUS=0
          SEGINI IPT16
          IPT16.ITYPEL=48
          DO 678 MM11=1,IPT8.NUM(/2)
            IF (IPT8.NUM(MLONG,MM11).NE.6) GOTO 678
            MCOL=MCOL+1
            DO 677 MTMP=1,9
              IPT16.NUM(MTMP,MCOL)=IPT8.NUM(MTMP,MM11)
              IPT16.ICOLOR(MCOL)=6
 677        CONTINUE
 678      CONTINUE
        ENDIF
C
      ENDIF
C
C--------------------------------------
C D) Creation du maillage de relations !
C--------------------------------------
C 1/ Initialisation de MAIREL
C MAIREL contient tous les maillages de relations, certains étant vides
C car il n'y a aucune relation du type considéré dans le nveau maillage
      NBMAIL=6
      SEGINI MAIREL
      MAIREL(1)=IPT11
      MAIREL(2)=IPT12
      MAIREL(3)=IPT13
      MAIREL(4)=IPT14
      MAIREL(5)=IPT15
      MAIREL(6)=IPT16
C
C 2/ Initialisation de LTYPE2
C LTYPE2 va contenir tous les sous-objets du nouveau maillage,
C relations et non relations
      NBLIS2=MCONF+MAX(1,MELEME.LISOUS(/1))
      SEGINI LTYPE2
C
C 3/ Renseignement des lisous non 48 dans LTYPE2
C S'il y a N sous-objets de non-relations dans le nouveau maillage,
C les N-ièmes premières colonnes de LTYPE2 contiennent ces sous-objets
      DO 685 IZZ=1,MAX(1,MELEME.LISOUS(/1))
        IF (MELEME.LISOUS(/1).EQ.0) THEN
          IPT2=MELEME
        ELSE
          IPT2=MELEME.LISOUS(IZZ)
          SEGACT IPT2
        ENDIF
        LTYPE2(IZZ,1)=IZZ
        LTYPE2(IZZ,2)=IPT2.ITYPEL
        LTYPE2(IZZ,3)=IPT2
        LTYPE2(IZZ,4)=0
 685  CONTINUE
C
C 4/ Renseignement des lisous 48 dans LTYPE2
C Les colonnes suivantes contiennent les sous-objets de relations,
C chaque sous-objet correspondant a un type de relations de conformité
      IZZ=NBLIS2-MCONF
      DO 686 IXX=1,6
        IF (MAIREL(IXX).EQ.0) GOTO 686
        IPT2=MAIREL(IXX)
        IZZ=IZZ+1
        LTYPE2(IZZ,1)=IZZ
        LTYPE2(IZZ,2)=IPT2.ITYPEL
        LTYPE2(IZZ,3)=IPT2
        LTYPE2(IZZ,4)=IPT2.ICOLOR(1)
 686  CONTINUE
C--------------------------------------
C E)  Renseignement de LTYPE2 interne a
c    la boucle 33 dans LTYPE3 externe a
c     la boucle 33
C--------------------------------------C
C
 
      LTYPE3(INB,1)=LTYPE2
      IF (MCONF.EQ.0) THEN
         LTYPE3(INB,2)=0
      ELSE
         LTYPE3(INB,2)=NBLIS2-MCONF+1
      ENDIF
c fin de boucle sur les sous domaines d'ipt5
  33  CONTINUE
c      SEGSUP IPT2
 
 
C          Mise a jour des maillages en entrée sortie de la boucle 9
C          IPT5, IPT6, IPT7
C=======================================================================
 679  CONTINUE
C--------------------------------------C
C A)  initialisation
C--------------------------------------C
 
c 1/ initialisaton de IPT7 maillage stockant tout les sous maillage
c à la fin d'une itération de la boucle 9
      NBS7=0
      DO 331 I3=1,LTYPE3(/1)
        LTYPE2=LTYPE3(I3, 1)
        SEGACT LTYPE2
        NBS7=NBS7+LTYPE2(/1)+1
 331  CONTINUE
      IF (NBS7.EQ.1) NBS7=0
      NBNN=0
      NBELEM=0
      NBREF=0
      NBSOUS=NBS7
      SEGINI IPT7
 
c 2/ initialisaton de IPT6 maillage stockant tout les sous maillage
c de relation à la fin d'une itération de la boucle 9
      NBS6=0
 
      DO 332 I3=1,LTYPE3(/1)
        LTYPE2=LTYPE3(I3, 1)
        SEGACT LTYPE2
        IF (LTYPE3(I3,2).NE.0) NBS6=NBS6+(LTYPE2(/1)-LTYPE3(I3,2)+1)
 332  CONTINUE
      IF (NBS6.EQ.1) NBS6=0
      NBNN=0
      NBELEM=0
      NBREF=0
      NBSOUS=NBS6
      SEGINI IPT6
c 3/ initialisaton de IPT5 maillage stockant tout les sous maillage
c non relation à la fin d'une itération de la boucle 9
      NBS5=0
      DO 333 I3=1,LTYPE3(/1)
        LTYPE2=LTYPE3(I3, 1)
        SEGACT LTYPE2
        IF (LTYPE3(I3,2).EQ.0) THEN
          NBS5=NBS5+LTYPE2(/1)
        ELSE
          NBS5=NBS5+(LTYPE3(I3,2)-1)
        ENDIF
 333  CONTINUE
      IF (NBS5.EQ.1) NBS5=0
      NBNN=0
      NBELEM=0
      NBREF=0
      NBSOUS=NBS5
      SEGINI IPT5
C--------------------------------------C
C B)  assemblage des nouveaux maillages a
c     partir des sous maillages stockes
c     dans LTYP3
C--------------------------------------C
      IND5=0
      IND6=0
      IND7=0
c 1/assemblage des nouveaux maillages non relation a partir des sous
c  maillages stockes dans LTYP3
 
      DO 35 I3=1,LTYPE3(/1)
c extraction d'un LTYPE2 de LTYPE3
        LTYPE2=LTYPE3(I3, 1)
        SEGACT LTYPE2
c        NBNN=0
c        NBELEM=0
c        NBREF=0
c        NBSOUS=0
        DO 613 I2=1,MAX(1,(LTYPE3(I3,2)-1))
c extraction d'un MAILLAGE non relation de LTYPE2
          IPT2=LTYPE2(I2,3)
          SEGACT IPT2
          IND5=IND5+1
          IF (NBS5.EQ.0) THEN
            IPT5=IPT2
c            write (*,*) 'IPT5.' , IPT5
c            write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
          ELSE
            IPT5.LISOUS(IND5)=IPT2
c            write (*,*) 'IPT5.LISOUS(',IND5,')' , IPT5.LISOUS(IND5)
c            write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
          ENDIF
          IND7=IND7+1
          IF (NBS7.EQ.0) THEN
            IPT7=IPT2
          ELSE
            IPT7.LISOUS(IND7)=IPT2
c            write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7)
c            write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
          ENDIF
 613    CONTINUE
 35   CONTINUE
 
c 2/assemblage des nouveaux maillage de relation à partir des derniers
c  sous maillages stockés dans LTYPE3
 
      DO 36 I3=1,LTYPE3(/1)
c extraction d'un LTYPE2 de LTYPE3
c        write (*,*) 'I3', I3
c        write (*,*) 'LTYPE3(I3,2)', LTYPE3(I3,2)
          IF (LTYPE3(I3,2).EQ.0) GOTO 36
          SEGACT LTYPE2
          LTYPE2=LTYPE3(I3, 1)
          NBNN=0
          NBELEM=0
          NBREF=0
          NBSOUS=0
          DO 614 I2=LTYPE3(I3,2), LTYPE2(/1)
c extraction d'un MAILLAGE de relation de LTYPE2
c            write (*,*) 'I2', I2
            SEGINI IPT2
            IPT2=LTYPE2(I2,3)
            SEGACT IPT2
c            write (*,*) 'IPT2: Type', IPT2.ITYPEL,'nbelem', IPT2.NUM(/2)
c            write (*,*) 'coul IPT2' , IPT2.ICOLOR(1)
            IND6=IND6+1
            IF (NBS6.EQ.0) THEN
              IPT6=IPT2
c              write (*,*) 'IPT6' , IPT6
c              write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
           ELSE
              IPT6.LISOUS(IND6)=IPT2
 
c          write (*,*) 'IPT6.LISOUS(',IND6,')' , IPT6.LISOUS(IND6)
c          write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
        ENDIF
 614    CONTINUE
 36   CONTINUE
 
c 3/ Concatenation des sous maillages de relation de même couleur
 
c  GG 11 juin 2018 : attention il peut encore y avoir des problemes a la
c jonctions entre deux sous zones dans le cas ou on raffine plusieurs
c types d'element differents (passer cette etape de concatenation dans
c la boucle 33 ?
 
      IF (NBS6.EQ.0) then
         IND7=IND7+1
         IPT7.LISOUS(IND7) = IPT6
c          write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7)
c          write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
 
         goto 370
      ENDIF
      NBS60 = 0
      SEGACT IPT6
 
c boucle sur les couleurs
      NBELEM = 0
      NBNN= 0
      NBREF= 0
      NBSOUS=6
      SEGINI IPT4
      DO 37 ICOUL = 1, 6
c boucle sur les sous maillages de ipt6
          NCOl = 0
          DO 38 I6 = 1, NBS6
                IPT3 = IPT6.LISOUS(I6)
                SEGACT IPT3
                IF ((IPT3.icolor(1)).EQ.ICOUL) then
                   NCOL = NCOL+1
                   IF (NCOL.EQ.1) THEN
                      SEGINI IPT2
                      IPT2 = IPT6.LISOUS(I6)
 
                   ENDIF
                   IF (NCOL.GE.2) THEN
 
c                      IPT3 = IPT6.LISOUS(I6)
c                      SEGACT IPT3
c Plusieurs sous maillages de même couleur on concatene IPT2 et IPT3
                      NBELEM = IPT2.NUM(/2)
                      NBNN= IPT2.NUM(/1)
                      NBREF=0
                      NBSOUS=0
                      NBELE0 = NBELEM
                      DO 381 IE3 = 1, IPT3.NUM(/2)
                      DO 382  IE2 = 1, NBELE0
                    IF ((IPT3.NUM(1,IE3)).EQ.(IPT2.NUM(1,IE2))) GOTO 381
 382                  CONTINUE
                      NBELEM = NBELEM +1
                      SEGADJ IPT2
                      DO IN2 = 1, NBNN
                         IPT2.NUM(IN2 , NBELEM) = IPT3.NUM(IN2, IE3)
                         IPT2.icolor(NBELEM) = IPT3.icolor(IE3)
                      ENDDO
 381                  CONTINUE
                   ENDIF
                ENDIF
 38        CONTINUE
          IF (NCOL.GE.1) THEN
              IPT4.LISOUS(ICOUL) = IPT2
              NBS60= NBS60+1
          ENDIF
 
 37   CONTINUE
 
c mise a jour de ipt6 et ipt7
      IF (NBS60.EQ.1) then
         IPT6 = IPT2
         segact ipt6
 
c         write (*,*) 'IPT6' , IPT6
c         write (*,*) 'type', IPT6.ITYPEL, 'Nbel' , IPT6.num(/2)
         IND7=IND7+1
         IPT7.LISOUS(IND7) = IPT2
c          write (*,*) 'IPT7.LISOUS(',IND7,')' , IPT7.LISOUS(IND7)
c          write (*,*) 'type', IPT2.ITYPEL, 'Nbel' , IPT2.num(/2)
      ELSE
        SEGACT IPT4
        NBELEM = 0
        NBNN = 0
        NBREF = 0
        NBSOUS = NBS60
        SEGINI IPT6
        IS6=0
        DO ICOUL = 1, 6
          IF ((IPT4.LISOUS(ICOUL)).NE.0) THEN
            IS6 = IS6+1
            IPT6.LISOUS(IS6)= IPT4.LISOUS(ICOUL)
c            write (*,*) 'IPT6.LISOUS(',IND6,')' , IPT6.LISOUS(IND6)
c            write (*,*) 'type', IPT4.ITYPEL, 'Nbel' , IPT4.num(/2)
 
            IND7=IND7+1
            IPT7.LISOUS(IND7) = IPT4.LISOUS(ICOUL)
          ENDIF
        ENDDO
      ENDIF
 370   CONTINUE
      NBELEM = 0
      NBNN = 0
      NBREF = 0
      NBSOUS = IND7
 
      SEGADJ IPT7
 
      SEGINI LTYPE4
      LTYPE4 = LTYPE3
 
c  mise a jour du nombre global de relations
      LCONF3 = MCONF
 
C=======================================================================
C              Test sur la progression du raffinement
C=======================================================================
C
C     CB215821 Retrait du caractere *MOD le cas echeant
      SEGACT IPT5,MCHPO1,MSOUPO,MPOVA1,IPT3,KMILIE
c      SEGSUP ICPR,KARETE,KARPOS,KARET2,ICHPO,IDENSI
c      SEGSUP ICPR,KARETE,KARPOS,KARET2,IDENSI
      SEGSUP ICPR,KARETE,KARPOS,KARET2
      SEGSUP MCHAM1,MELVA1,MCHAM2,MELVA2,IMODE1
      SEGSUP MMODE1,MCHEL1,MCHEL2
 
      IF ((IDIM.EQ.3).AND.(NBOUCL.NE.1)) SEGSUP JPLANS,JPLAN2,JPLAN3,
     +JPLCOM,JNOEFA,IWORK2,JFARAF, IPT9
 
C S'il y a eu a raffiner des elements lors de l'iteration courante,
C on passe a l'iteration suivante
C
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c   forcage  : nb max itération de raff
c      IF (NBOUCL.EQ.3) THEN
c           KARAF2=0
c           KARAF=0
c      ENDIF
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      IF (KARAF2.NE.0) GOTO 9
 
C ### SORTIE TEMPORAIRE DU PROGRAMME ###
C Si on est a la premiere iteration de RAFF et qu'il n'y a pas
C d'elements a raffiner, on retourne en sortie le maillage de depart,
C donné en entree
      IF ((KARAF2.EQ.0).AND.(NBOUCL.EQ.1)) THEN
        WRITE(*,*) 'il n est pas necessaire de raffiner le maillage'
        MELEME=IPT1
        GOTO 998
      ENDIF
 
C=======================================================================
C                     Mises a jour diverses
C=======================================================================
C **  On affecte une densite a chaque noeud que l'on a cree
c      DO 37 I=NPTINI+1,XCOOR(/1)/(IDIM+1)
c        XCOOR(I*(IDIM+1))=DENSI(ICPR(I,1))
c  37  CONTINUE
 
C=======================================================================
C                      Fin du programme
C=======================================================================
 
      MELEME=IPT7
 
      SEGSUP XDEN
 998  CONTINUE
 
      CALL ACTOBJ('MAILLAGE',MELEME,1)
      CALL ECROBJ('MAILLAGE',MELEME)
 
c      RETURN
      END