C RAFF3D    SOURCE    CB215821  26/10/26    21:16:02     12682          
      SUBROUTINE RAFF3D(IPT2,IPT3,ICPR,KARPOS,KARETE,KMILIE,MELVA2,
     .NACREE,KARAF,IPT4,JPLANS,JPLAN3,JPLCOM,JNOEFA,IPT7,JFARAF,KARET2,
     .XDEN)
c entrees de raff3D:
c           - IPT2:   maillage élémentaire à raffiner 
c           - IPT3:   maillage de SURE précédant 
c           - ICPR:   tableau de passage noeuds local/global
c           - KARPOS: tableau du nb d'arretes par noeuds
c           - KARETE: tableau des d'arretes KARETE(i,j)=k 
c           - KMILIE: tableau des hanging nodes associés au arêtes
c                     si ancienne relation 
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           - NACREE: nombre de noeuds à créer dans ipt2
c           - KARAF:  nombre d'élément à raffiner dans ipt2
c           - JPLANS : tableau des faces JPLAN(i,j)=k
c           - JPLAN3 : tableau du nombre de relation de compat par face
c           - JPLCOM : tableau du nb de faces par noeuds
c           - JNOEFA : tableau des noeuds et éléments associé au faces
c           - JFARAF : tableau des hanging nodes associés au faces pour
c                      les anciennes relations.
c           - XDEN :  densité aux noeuds en notation globale
c           - KARET2: tableau du nombre d'éléments touchant les arêtes
c sorties: 
c           - KMILIE: tableau des hanging nodes associés au arê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 à supprimer
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           - JPLAN3 : tableau du nombre de relation de compat par face
c                      mis a jour
c           - JFARAF : tableau des hanging nodes associés au faces pour
c                      mis a jour.
c           - MELEME: Maillage élémentaire raffiné à partir de ipt2
c                     sans relations
c           - IPT7 :  Maillage de relations créé incomplet 
c           - XDEN :  densité aux noeuds en notation globale
c                     interpolé aux nouveaux noeuds.
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER(I-N)

-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC CCGEOME
-INC SMELEME
-INC SMCHPOI
-INC SMMODEL
-INC SMCHAML
C
C======================================================================
C                           Declarations
C======================================================================
      SEGMENT ICPR(nbpts,2)
      SEGMENT XDEN(nbpts)
      SEGMENT KARETE(NBNDS,NCOL)
      SEGMENT KARET2(NBNDS,NCOL)
      SEGMENT KMILIE(NBNDS,NCOL)
      SEGMENT KARPOS(NBNDS)
      SEGMENT JPLANS(JPLA1,JPLA2)
      SEGMENT JPLAN3(JPLA1,JPLA2)
      SEGMENT JPLCOM(JPLA1)
      SEGMENT JNOEFA(JNBFA,5)
      SEGMENT JFARAF(JNBFA,LLLL)
      SEGMENT NUMNOE(INUMNO)
      SEGMENT IWORK2(JNBSOM)
      SEGMENT IWORK1(IJKLMN)
      SEGMENT IWORK4(JNBSOM)
      SEGMENT IWORK3(IJKLMN)

      REAL*8 LON68,LON510,LON47,LON29
C
C======================================================================
C                        Initialisations
C======================================================================
      SEGACT JPLANS,JPLCOM,JNOEFA*MOD,JPLAN3*MOD, IPT3
      SEGACT IPT2,ICPR,KARPOS,KARETE,KMILIE*MOD,MELVA2,JFARAF*MOD
      IJKLMN=4
      ipt7=0
      SEGINI IWORK1
      IWORK1(1)=4
      IWORK1(2)=8
      IWORK1(3)=6
      IWORK1(4)=10
C      JNBFA=JNOEFA(/1)
C      LLLL=2
C      IF (IPT2.ITYPEL.EQ.15) LLLL=10
C      IF ((IPT2.ITYPEL.EQ.17).OR.(IPT2.ITYPEL.EQ.24)) LLLL=6
C      LLLL=10
C      SEGINI JFARAF
      NCOMPL=LNELM(1,(IPT2.ITYPEL-1)*2+2)
      IF (NCOMPL.NE.0) GOTO 999
      NBNN=NBNNE(LNELM(2,(IPT2.ITYPEL-1)*2+1))
      INELM=LNELM(1,(IPT2.ITYPEL-1)*2+1)
C
C Creation du squelette du maillage resultat
      NBELEM=IPT2.NUM(/2)-KARAF+INELM*KARAF
      NBSOUS=0
      NBREF=0
      SEGINI IPT4
      IPT4.ITYPEL=LNELM(2,(IPT2.ITYPEL-1)*2+1)
      LPO2=LPOS2(IPT2.ITYPEL)
      segact mcoord*mod
      NBPT0=nbpts
      NBPT1=nbpts
      NBPTS=NBPT1+NACREE+KARAF*NBINTE(IPT2.ITYPEL)
      SEGADJ MCOORD,XDEN
      INUMNO=NBRAF(IPT2.ITYPEL)
      SEGINI NUMNOE
      LCOMP=1
C
C======================================================================
C                   Phase de raffinement 3D
C======================================================================
C=======================================
C A) Boucle sur les elements a raffiner !
C=======================================
      DO 6 IARAF=1,MELVA2.VELCHE(/2)
        IF (MELVA2.VELCHE(1,IARAF).NE.1) THEN
          DO 1 IJKL=1,IPT2.NUM(/1)
            IPT4.NUM(IJKL,NBELEM)=IPT2.NUM(IJKL,IARAF)
   1      CONTINUE
          NBELEM=NBELEM-1
          GOTO 6
        ENDIF
        JCOMPT=0
        JPOS5=LPOS5(IPT2.ITYPEL)
C
C==========================================================
C B) Boucle sur les noeuds a creer pour raffiner l'element !
C==========================================================
        DO 4 I=1,NBRAF(IPT2.ITYPEL)
          JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
          JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
          JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
          LTYPNO=JTYPNO(JPOS5-1+I)
c          write(*,*) 'LTYPNO', LTYPNO
C
C-------------------------------
C ** B.1 / On est sur une arete !
C-------------------------------
          IF (LTYPNO.EQ.0) THEN
            NPTA=IPT2.NUM(LISNOE(JLISN),IARAF)
            NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF)
            NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
            NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
            DO 2 K=1,MAX(1,KARPOS(NMIN))
              IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K
   2        CONTINUE
            IF (KMILIE(NMIN,NEXIST).GT.0) THEN
c si le noeud existe deja on rempli numnoe avec Kmilie
              NUMNOE(I)=KMILIE(NMIN,NEXIST)
              JCOMPT=JCOMPT+JLONG
              GOTO 4
            ELSE
c sinon on ajoute un noeud a la fin de numnoe et on le met dans  kmilie
              NBPT1=NBPT1+1
              NUMNOE(I)=NBPT1
              KMILIE(NMIN,NEXIST)=NBPT1
            ENDIF

c            IF ((NBPT1.eq.10330).or.(NBPT1.eq.10341)) THEN
c               WRITE(*,*) 'NUMNOE(I)', NUMNOE(I) 
c               WRITE(*,*) 'NPTA', NPTA, 'NPTB', NPTB
c               WRITE(*,*) 'NMIN', NMIN, 'NMAX', NMAX ,'NEXI', NEXIST
c               WRITE(*,*) 'KMILIE(NMIN,NEXIST)', KMILIE(NMIN,NEXIST)
c               WRITE(*,*) 'KARETE(NMIN,NEXIST)', KARETE(NMIN,NEXIST)
c            endif
C
C------------------------------
C ** B.2 / On est sur une face !
C------------------------------
C Il faut identifier cette face (numero, sommets...)
          ELSEIF ((LTYPNO.GT.0).AND.(LTYPNO.LT.7)) THEN
C => a) On initialise toutes les donnees relatives a cette face
            JTYPEL=IPT2.ITYPEL
            JLTEL2=LTEL(1,JTYPEL)-1+LTYPNO
            JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO
            JLDEL1=LDEL(1,JLTEL2)
            JTYFAC=IWORK1(JLDEL1)
            JLDEL2=LDEL(2,JLTEL2)
            JNBSOM=NBSOM(JTYFAC)
            JSPOS=NSPOS(JTYFAC)
C => b) On identifie les sommets de la face (n° global)
            SEGINI IWORK2
            DO 10 IAA=1,JNBSOM
              NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF)
              IWORK2(IAA)=NGLOBA
  10        CONTINUE
C => c) On classe ces sommets par ordre croissant (NPTA < NPTB < NPTC)
            NPTA=nbpts+1
            NPTB=NPTA+1
            NPTC=NPTB+1
            DO 11 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
  11        CONTINUE
C => d) On passe ces sommets en n° 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 => e) On cherche le numero de la face
            NEXIS2=0
            DO 12 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
  12        CONTINUE
            JNUMFA=JPLANS(NPTA,NEXIS2)
C
C---------------------------------
C ** B.3 / Raffinement de la face !
C---------------------------------
            KTEST1=JPLAN3(NPTA,NEXIS2)
            KTEST2=NBINTE(JTYFAC)
C => a) Si la face est de type QUA4 (un seul noeud a creer, au milieu)
            IF (JTYFAC.EQ.8) THEN
              IF (KTEST1.LT.KTEST2) THEN
c    a.1) si le noeud milieu n'existe pas on l'ajoute a la fin de numnoe
c          et on rempli JFARAF et JPLAN3
                NBPT1=NBPT1+1
                NUMNOE(I)=NBPT1
                JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1
                JFARAF(JNUMFA,1)=NBPT1
              ELSEIF (KTEST1.EQ.KTEST2) THEN
c    a.2) si le noeud milieu existe deja on rempli numnoe grace à JFARAF
c         et on met celui-ci a zéro
                NUMNOE(I)=JFARAF(JNUMFA,1)
                JFARAF(JNUMFA,1)=0
                JCOMPT=JCOMPT+JLONG
                GOTO 4
              ELSE
                WRITE(*,*) 'ERREUR'
              ENDIF
c              IF ((NBPT1.eq.10248).or.(NBPT1.eq.10276)) THEN
c              IF ((NPTA.eq.235).or.(NPTA.eq.236).or.(NPTA.eq.237)) THEN
c                WRITE(*,*) 'NUMNOE(I)', NUMNOE(I) 
c                WRITE(*,*) 'NPTA', NPTA1, 'NPTB', NPTB1, 'NPTC', NPTC1
c                WRITE(*,*) 'NPTA', NPTA, 'NPTB', NPTB, 'NPTC', NPTC
c                WRITE(*,*) 'JNUMFA', JNUMFA, 'NEXI', NEXIS2
c                WRITE(*,*) 'JPLAN3(NPTA,NEXIS2)', JPLAN3(NPTA,NEXIS2)
c                WRITE(*,*) 'JFARAF(JNUMFA,1)', JFARAF(JNUMFA,1)
c                WRITE(*,*) 'JNOEFA(JNUMFA,1)', JNOEFA(JNUMFA,1)
c                WRITE(*,*) 'JNOEFA(JNUMFA,2)', JNOEFA(JNUMFA,2)
c                WRITE(*,*) 'JNOEFA(JNUMFA,3)', JNOEFA(JNUMFA,3)
c                WRITE(*,*) 'JNOEFA(JNUMFA,4)', JNOEFA(JNUMFA,4)
c                WRITE(*,*) 'JNOEFA(JNUMFA,5)', JNOEFA(JNUMFA,5)
c              endif
            ENDIF
C => b) Si la face n'est pas de type QUA4 (donc de type TRI6 ou QUA8)
            IF (JTYFAC.NE.8) THEN
C    b.1)   Si il manque des noeuds dans la face
C           On ajoute un noeud à la fin de numnoe
              IF (KTEST1.LT.KTEST2) THEN
                NBPT1=NBPT1+1
                NUMNOE(I)=NBPT1
                JPLAN3(NPTA,NEXIS2)=JPLAN3(NPTA,NEXIS2)+1
                NEXIS3=0
                XCO2=0.25
                DO 13 KBB=1,JLONG
                  XCO1=XCOEFF(JPOS1-1+KBB)
                  IF (XCO1.EQ.XCO2) NEXIS3=KBB
  13            CONTINUE
                JFARAF(JNUMFA,2*(KTEST1)+1)=NBPT1
                IF (NEXIS3.NE.0) THEN
C    b.1.1)     si le point est sur une arête du triangle interieur 
C               (pour un TRI6)
C               si le point est sur une arête de la coie interieure 
C               (pour un QUA8)

C               on renseigne la partie paire de JFARAF qui doit déja 
C               exister ????
                  NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF)
                  JFARAF(JNUMFA,2*(KTEST1)+2)=NGLOB
                ELSE
C    b.1.2)     si le point est le milieu d'un QUA8

C               on met la partie paire de JFARAF à zéro
                  JFARAF(JNUMFA,2*(KTEST1)+2)=0
                ENDIF
              ELSEIF (KTEST1.EQ.KTEST2) THEN
C    b.2)     Si il ne manque pas de noeuds dans la face
                NEXIS3=0
                NEXIS4=0
                XCO2=0.25
                DO 14 KBB=1,JLONG
                  XCO1=XCOEFF(JPOS1-1+KBB)
                  IF (XCO1.EQ.XCO2) NEXIS3=KBB
  14            CONTINUE
                IF (NEXIS3.NE.0) THEN
C    b.2.1)     si le point est sur une arête du triangle interieur 
C               (pour un TRI6)
C               si le point est sur une arête de la croie interieure 
C               (pour un QUA8)
                  NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF)
                  DO 15 KAA=2,JFARAF(/2),2
                     IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA
  15              CONTINUE
c                 on rempli numnoe grace à JFARAF et on met la partie 
c                 associée au noeud en question de celui-ci a 0
                  NUMNOE(I)=JFARAF(JNUMFA,NEXIS4-1)
                  JFARAF(JNUMFA,NEXIS4-1)=0
                ELSE
C    b.2.2)     si le point est le milieu d'un QUA8
                  DO 16 KAA=2,JFARAF(/2),2
                     IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA
  16              CONTINUE
c                 on rempli numnoe grace à JFARAF et on met la partie 
c                 associée au noeud en question de celui-ci a 0
                  NUMNOE(I)=JFARAF(JNUMFA,NEXIS4-1)
                  JFARAF(JNUMFA,NEXIS4-1)=0
                ENDIF
                JCOMPT=JCOMPT+JLONG
                GOTO 4
              ENDIF
            ENDIF
C
C------------------------------------------------------
C ** B.4 / On est a l'interieur du volume de l'element !
C------------------------------------------------------
          ELSEIF (LTYPNO.EQ.7) THEN
            NBPT1=NBPT1+1
            NUMNOE(I)=NBPT1
          ENDIF
C allocation de plus de mémoire dans xcoord et xden si nécéssaire
          IF (NBPT1.EQ.NBPTS) THEN
            NBPTS=NBPTS+200
            SEGADJ MCOORD,XDEN
          ENDIF
C
C==============================
C C) Creation du nouveau point !
C==============================
C On continue ici que lorsque l'on doit creer un nouveau point
          XPT=0.
          YPT=0.
          ZPT=0.
          XDEN1=0.D0
          XDENMIN = 1000.
          DO 3 J=1,JLONG
            NGLOB=IPT2.NUM(LISNOE(JLISN-1+J),IARAF)
            XINI=XCOOR((NGLOB-1)*(IDIM+1)+1)
            YINI=XCOOR((NGLOB-1)*(IDIM+1)+2)
            ZINI=XCOOR((NGLOB-1)*(IDIM+1)+3)
            XPT=XPT+XINI*XCOEFF(JPOS1-1+J)
            YPT=YPT+YINI*XCOEFF(JPOS1-1+J)
            ZPT=ZPT+ZINI*XCOEFF(JPOS1-1+J)
            
            XDEN1=XDEN1+XDEN(NGLOB)*XCOEFF(JPOS1-1+J)
            IF (XDEN(NGLOB).LT.XDENMIN) XDENMIN = XDEN(NGLOB) 
c          IF (NBPT1.eq.10219) THEN
c            write (*,*) 'nglob', nglob ,'nbpt1', nbpt1
c            write (*,*) 'xpt ', xpt
c            write (*,*) 'ypt ', ypt
c            write (*,*) 'zpt ', zpt
c            write (*,*) 'xden1 ', xden1
c            write (*,*) 'xcoeff1 ', XCOEFF(JPOS1-1+J)
c          endif
   3      CONTINUE

c seuil de densité a xdenmin
          IF (XDEN1.LT.XDENMIN) XDEN1 = XDENMIN

          XCOOR((NBPT1-1)*(IDIM+1)+1)=XPT
          XCOOR((NBPT1-1)*(IDIM+1)+2)=YPT
          XCOOR((NBPT1-1)*(IDIM+1)+3)=ZPT
          XCOOR((NBPT1-1)*(IDIM+1)+4)=XDEN1
          XDEN(NBPT1)=XDEN1
          JCOMPT=JCOMPT+JLONG

C======================================================================
   4    CONTINUE
C======================================================================

        JPOS4=LPOS4(IPT2.ITYPEL)
C
C Cas des tetraedres
        LONDIA=0
        JTYPEL=IPT2.ITYPEL
        IF (JTYPEL.EQ.23) THEN
          LPT6=NUMNOE(6)
          LPT5=NUMNOE(5)
          LPT8=NUMNOE(8)
          LPT10=NUMNOE(10)
          XPT6=XCOOR((LPT6-1)*(IDIM+1)+1)
          XPT5=XCOOR((LPT5-1)*(IDIM+1)+1)
          XPT8=XCOOR((LPT8-1)*(IDIM+1)+1)
          XPT10=XCOOR((LPT10-1)*(IDIM+1)+1)
          YPT6=XCOOR((LPT6-1)*(IDIM+1)+2)
          YPT5=XCOOR((LPT5-1)*(IDIM+1)+2)
          YPT8=XCOOR((LPT8-1)*(IDIM+1)+2)
          YPT10=XCOOR((LPT10-1)*(IDIM+1)+2)
          ZPT6=XCOOR((LPT6-1)*(IDIM+1)+3)
          ZPT5=XCOOR((LPT5-1)*(IDIM+1)+3)
          ZPT8=XCOOR((LPT8-1)*(IDIM+1)+3)
          ZPT10=XCOOR((LPT10-1)*(IDIM+1)+3)
          LON68=SQRT((XPT8-XPT6)**2+(YPT8-YPT6)**2+(ZPT8-ZPT6)**2)
          LON510=SQRT((XPT10-XPT5)**2+(YPT10-YPT5)**2+(ZPT10-ZPT5)**2)
          IF (LON510.LT.LON68) LONDIA=4*4
        ENDIF
        IF (JTYPEL.EQ.24) THEN
          LPT2=NUMNOE(2)
          LPT4=NUMNOE(4)
          LPT7=NUMNOE(7)
          LPT9=NUMNOE(9)
          XPT2=XCOOR((LPT2-1)*(IDIM+1)+1)
          XPT4=XCOOR((LPT4-1)*(IDIM+1)+1)
          XPT7=XCOOR((LPT7-1)*(IDIM+1)+1)
          XPT9=XCOOR((LPT9-1)*(IDIM+1)+1)
          YPT2=XCOOR((LPT2-1)*(IDIM+1)+2)
          YPT4=XCOOR((LPT4-1)*(IDIM+1)+2)
          YPT7=XCOOR((LPT7-1)*(IDIM+1)+2)
          YPT9=XCOOR((LPT9-1)*(IDIM+1)+2)
          ZPT2=XCOOR((LPT2-1)*(IDIM+1)+3)
          ZPT4=XCOOR((LPT4-1)*(IDIM+1)+3)
          ZPT7=XCOOR((LPT7-1)*(IDIM+1)+3)
          ZPT9=XCOOR((LPT9-1)*(IDIM+1)+3)
          LON47=SQRT((XPT7-XPT4)**2+(YPT7-YPT4)**2+(ZPT7-ZPT4)**2)
          LON29=SQRT((XPT9-XPT2)**2+(YPT9-YPT2)**2+(ZPT9-ZPT2)**2)
          IF (LON29.LT.LON47) LONDIA=4*10
        ENDIF
C
C===================================
C D) Creation des nouveaux elements !
C===================================
C On remplit la portion de IPT4 relative aux elements crees a partir
C de la division de l'element IARAF (indice de boucle 1).
C Cette portion de IPT4 contient les colonnes dont la valeur s'etend
C de INELM*(LCOMP-1)+1 a INELM*LCOMP.
        DO 1000 J=1,INELM
        DO 5 I=1,IPT4.NUM(/1)
          NTEMP=LIELM(JPOS4-1+NBNN*(J-1)+I)
          IF (((JTYPEL.EQ.23).OR.(JTYPEL.EQ.24)).AND.(J.GT.4)) THEN
            NTEMP=LIELM(JPOS4-1+NBNN*(J-1)+I+LONDIA)
          ENDIF
          IF (NTEMP.GT.NBNN) THEN
            IPT4.NUM(I,INELM*(LCOMP-1)+J)=NUMNOE(NTEMP-NBNN)
          ELSE
            IPT4.NUM(I,INELM*(LCOMP-1)+J)=IPT2.NUM(NTEMP,IARAF)
          ENDIF
   5    CONTINUE
 1000   CONTINUE
        LCOMP=LCOMP+1
   6  CONTINUE
C
      NBPT =NBPT1
C      SEGADJ MCOORD,XDEN
C
C=======================================================================
C              Preparation du maillage de relations
C=======================================================================
      SEGSUP KARET2
      NBNDS=KARETE(/1)
      NCOL=KARETE(/2)
      SEGINI KARET2
C==================================
C A) Relations dues aux faces (3D) !
C==================================
C Tous les noeuds qui restent dans JFARAF sont a creer en tant que
C relations de conformite ou des bords du domaine

C-----------------------------------------------------------------------
C ** A.1 / Prise en compte des relation 
C          d'arretes pour les arretes de chaque faces 
C-----------------------------------------------------------------------

C-----------------------------------------------------------------------
C** A.1.1/ faces raffinée
C-----------------------------------------------------------------------
 
c      write (*,*) 'relations dues au faces'
      JTYPEL=IPT2.ITYPEL
      JLTEL1=LTEL(1,JTYPEL)
      JLTEL2=LTEL(2,JTYPEL)
      IF (JTYPEL.NE.14) GOTO 141
c Boucle sur les elements de l'ancien maillage
      DO 121 KELEM1=1,IPT2.NUM(/2)
c Boucle sur les faces de l'ancien maillage
        DO 122 IBB1=1,JLTEL1
          JLDEL1=LDEL(1,JLTEL2-1+IBB1)
          JTYFAC=IWORK1(JLDEL1)
          JLDEL2=LDEL(2,JLTEL2-1+IBB1)
          JNBSOM=NBSOM(JTYFAC)
          JSPOS=NSPOS(JTYFAC)
          SEGINI IWORK2
C
C 1/ trouver le numéro de la face
          DO 123 IAA=1,JNBSOM
              NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),KELEM1)
              IWORK2(IAA)=NGLOBA
 123      CONTINUE
C
C 1.1/ Classement des 3 sommets par ordre croissant de n° globale
          NPTA=NBPT +1
          NPTB=NPTA+1
          NPTC=NPTB+1
          DO 124 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
 124      CONTINUE
C
C1.2/ Classement des 3 sommets par ordre croissant de n° 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 1.3/ Recherche du numero de la face
          NEXIS2=0
          DO 125 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
 125      CONTINUE
          
          JNUMF1=JPLANS(NPTA,NEXIS2)

C
C 2/ Exclusion des faces non voulues 
   
C on exlu les faces où un noeud existait déja
          IF (JFARAF(JNUMF1,1).EQ.0) GOTO 122
C on exlu les faces de bords
          IF (JNOEFA(JNUMF1,5).EQ.0) GOTO 122

C 3/ Recherche des arretes de la faces pour y mettre des relations 

          IF (JTYPEL.NE.14) GOTO 137
          NBAR1=0
          DO 1001 ISS=1, (JNBSOM-1)
          DO 138 JSS=(ISS+1), JNBSOM

            NMIN=MIN(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
            NMAX=MAX(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))

            DO 139 IRR=1,MAX(1,KARPOS(NMIN))                   
              IF (KARETE(NMIN,IRR).EQ.NMAX) THEN
                NBAR1=NBAR1+1
                KARET2(NMIN,IRR)=KARET2(NMIN,IRR)+1
c                write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
c                IF (JFARAF(JNUMF1,1).EQ.15748) THEN
c                  write (*,*) ISS,JSS 
c                  write (*,*) '!!! NMIN', NMIN, 'NMAX', NMAX, '!!!'
c                  write (*,*) IWORK2(ISS), IWORK2(JSS)
c                  write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
c                  write (*,*) 'KARET2(NMIN,IRR)',KARET2(NMIN,IRR) 
c                ENDIF
              ENDIF
 139        CONTINUE
 138      CONTINUE 
 1001     CONTINUE
c          WRITE (*,*) 'Nbarrete', NBAR1
c          IF (NBAR1.NE.4) THEN
c             write (*,*) IWORK2(1), IWORK2(2)
c             write (*,*) IWORK2(3), IWORK2(4) 
c          ENDIF
 137      CONTINUE    
 122    CONTINUE        
 121  CONTINUE

 141  CONTINUE

C-----------------------------------------------------------------------
C** A.1.2/ relations de faces conservées dans le maillage précédent
C-----------------------------------------------------------------------

      SEGINI IPT5
      DO 670 MIE=1,MAX(1,IPT3.LISOUS(/1))
        IF (IPT3.LISOUS(/1).EQ.0) THEN
          IPT5=IPT3
        ELSE
          IPT5=IPT3.LISOUS(MIE)
        ENDIF
        SEGACT IPT5
c        write (*,*) 'ipt3 raff 3D' , ipt3, IPT3.ITYPEL, IPT3.ICOLOR(1)
c        write (*,*) 'ipt5 raff 3D' , ipt5, IPT5.ITYPEL, IPT5.ICOLOR(1)
        IF (IPT5.ITYPEL.NE.48) GOTO 670
        IF ((IPT5.ICOLOR(1).NE.3)) GOTO 670
        LIGNUM=IPT2.NUM(/1)
        JNBSOM=4
        SEGINI IWORK2

C
C 1/ trouver le numéro de la face
C 1.1/ Classement des 3 sommets par ordre croissant (num globale)
        DO 671,MAA=1,IPT5.NUM(/2)
          DO 673 MBB=1,JNBSOM
c            write (*,*) 'MAA' , MAA, 'MBB', MBB
            IWORK2(MBB)=IPT5.NUM(MBB+1,MAA)
            
 673      CONTINUE
          NPTA=nbpt +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 1.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 1.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
          JNUMF1=JPLANS(NPTA,NEXIS2)
c          IF (MAA.EQ.5) THEN
c                  write (*,*) 'MAA', MAA, 'nexis2', nexis2 
c                  write (*,*) IWORK2(1), IWORK2(2), IWORK2(3), IWORK2(4) 
c                  write (*,*) 'JFARAF(JNUMF1,1)', JFARAF(JNUMF1,1)
c                  write (*,*) 'JNOEFA(JNUMF1,5)',JNOEFA(JNUMF1,5) 
c          ENDIF
C
C 2/ Exclusion des faces non voulues
          IF (NEXIS2.EQ.0) THEN
c            WRITE (*,*)'FACE NON TROUVE SURE:', MAA 
            GOTO 671
          ENDIF
C on exlu les faces où un noeud existait déja
          IF (JFARAF(JNUMF1,1).EQ.0) GOTO 671

C on exlu les faces de bords
          IF (JNOEFA(JNUMF1,5).EQ.0) GOTO 671
          NBAR1=0
C
C 3/ Recherche des arretes de la faces pour y mettre des relations
         
          DO 1002 ISS=1, (JNBSOM-1)
          DO 135 JSS=(ISS+1), JNBSOM

            NMIN=MIN(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))
            NMAX=MAX(ICPR(IWORK2(ISS),1),ICPR(IWORK2(JSS),1))

            DO 136 IRR=1,MAX(1,KARPOS(NMIN))                   
              IF (KARETE(NMIN,IRR).EQ.NMAX) THEN
                NBAR1=NBAR1+1
                KARET2(NMIN,IRR)=KARET2(NMIN,IRR)+1
c                write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
c                IF (MAA.EQ.5) THEN
c                  write (*,*) ISS,JSS 
c                  write (*,*) '!!! NMIN', NMIN, 'NMAX', NMAX, '!!!'
c                  write (*,*) IWORK2(ISS), IWORK2(JSS)
c                  write (*,*) 'KMILIE(NMIN,IRR)', KMILIE(NMIN,IRR)
c                  write (*,*) 'KARET2(NMIN,IRR)',KARET2(NMIN,IRR) 
c                ENDIF
              ENDIF
 136        CONTINUE
 135      CONTINUE 
 1002     CONTINUE

 671    CONTINUE
   
 670  CONTINUE
      SEGDES IPT5
C---------------------------------
C ** A.2 / Initialisation de IPT5 !
C---------------------------------



C 1/ Comptage du nombre de noeuds soumis a des relations
      NBRELA=0
      DO 114 IHF=1,JFARAF(/1)
C on exlu les faces de bords
        IF (JNOEFA(IHF,5).EQ.0) GOTO 114
C on exlu les faces où un noeud existait déja
        IF (JFARAF(IHF,1).EQ.0) GOTO 114
        DO 113 JF=1,JFARAF(/2),2
          IF (JFARAF(IHF,JF).NE.0) NBRELA=NBRELA+1
 113    CONTINUE
 114  CONTINUE
C
C 2/ Creation de IPT5
      IF (NBRELA.EQ.0) THEN
        IPT5=0
        GOTO 210
      ENDIF
C      IF (IPT2.ITYPEL.EQ.14) NBNN=4
C      IF (IPT2.ITYPEL.EQ.15) NBNN=8
C      IF (IPT2.ITYPEL.EQ.17) NBNN=5
C      IF (IPT2.ITYPEL.EQ.24) NBNN=5
      NBNN=10
      NBELEM=NBRELA
      NBSOUS=0
      NBREF=0
      SEGINI IPT5
      IPT5.ITYPEL=48
C
C 3/ Renseignement des noeuds support des relations
      NBRELA=0
      DO 116 IPF=1,JFARAF(/1)
        IF (JNOEFA(IPF,5).EQ.0) GOTO 116
        IF (JFARAF(IPF,1).EQ.0) GOTO 116
        DO 115 JF=1,JFARAF(/2),2
          IF (JFARAF(IPF,JF).EQ.0) GOTO 115
          NBRELA=NBRELA+1
          IPT5.NUM(1,NBRELA)=JFARAF(IPF,JF)
 115    CONTINUE
 116  CONTINUE
c      write (*,*) 'NBRELA FACE ', NBRELA
C
C-----------------------------------------------------
C ** A.3 / Recherche des noeuds formant les relations !
C-----------------------------------------------------
C 1/ Boucle sur l'ensemble des noeuds a creer
      DO 200 IARAF=1,MELVA2.VELCHE(/2)
        IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 200
        JCOMPT=0
        JPOS5=LPOS5(IPT2.ITYPEL)
        DO 190 I=1,NBRAF(IPT2.ITYPEL)
          JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
          JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
          JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
          LTYPNO=JTYPNO(JPOS5-1+I)
          IF ((LTYPNO.EQ.0).OR.(LTYPNO.EQ.7)) GOTO 189
C
C 2/ Preparation pour trouver le noeud et la face en question
          JTYPEL=IPT2.ITYPEL
          JLTEL2=LTEL(2,JTYPEL)-1+LTYPNO
          JLDEL1=LDEL(1,JLTEL2)
          JTYFAC=IWORK1(JLDEL1)
          JLDEL2=LDEL(2,JLTEL2)
          JNBSOM=NBSOM(JTYFAC)
          JSPOS=NSPOS(JTYFAC)
          SEGINI IWORK2
C
C 3/ Classement des 3 sommets par ordre croissant de n° globale
          DO 100 IAA=1,JNBSOM
            NGLOBA=IPT2.NUM(LFAC(JLDEL2-1+IBSOM(JSPOS-1+IAA)),IARAF)
            IWORK2(IAA)=NGLOBA
 100      CONTINUE
          NPTA=nbpt +1
          NPTB=NPTA+1
          NPTC=NPTB+1
          DO 110 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
 110      CONTINUE
C
C 4/ Classement des 3 sommets par ordre croissant de n° 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 5/ Recherche du numero de la face
          NEXIS2=0
          DO 120 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
 120      CONTINUE
          JNUMFA=JPLANS(NPTA,NEXIS2)

C
C 6/ Recherche du numero global du point
          IF (JNOEFA(JNUMFA,5).EQ.0) GOTO 189
          IF (JTYFAC.EQ.8) INOEGL=JFARAF(JNUMFA,1)
          IF (JTYFAC.NE.8) THEN
            NEXIS3=0
            NEXIS4=0
            XCO2=0.25
            DO 140 KBB=1,JLONG
              XCO1=XCOEFF(JPOS1-1+KBB)
              IF (XCO1.EQ.XCO2) NEXIS3=KBB
 140        CONTINUE
            IF (NEXIS3.NE.0) THEN
              NGLOB=IPT2.NUM(LISNOE(JLISN-1+NEXIS3),IARAF)
              DO 150 KAA=2,JFARAF(/2),2
                IF (JFARAF(JNUMFA,KAA).EQ.NGLOB) NEXIS4=KAA
 150          CONTINUE
              INOEGL=JFARAF(JNUMFA,NEXIS4-1)
            ELSE
              DO 160 KAA=2,JFARAF(/2),2
                IF (JFARAF(JNUMFA,KAA).EQ.0) NEXIS4=KAA
 160          CONTINUE
              INOEGL=JFARAF(JNUMFA,NEXIS4-1)
            ENDIF
          ENDIF

C
C------------------------------
C ** A.4 / Remplissage de IPT5 !
C------------------------------
C 1/ Recherche de la position du point dans IPT5
          NEXIS5=0
          DO 170 IGG=1,IPT5.NUM(/2)
            IF (INOEGL.EQ.IPT5.NUM(1,IGG)) NEXIS5=IGG
 170      CONTINUE
          IF (NEXIS5.EQ.0) GOTO 189
          IF (INOEGL.EQ.139) THEN
c            write (*,*) INOEGL, IWORK2(1), IWORK2(2), IWORK2(3)  
          ENDIF

C
C 2/ Renseignement des points formant les relations
          DO 180 IHH=1,JLONG
            IPT5.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF)
 180      CONTINUE
          IF (JLONG.EQ.4) IPT5.NUM(10,NEXIS5)=3
          IF (JLONG.EQ.5) IPT5.NUM(10,NEXIS5)=4
          IF (JLONG.EQ.8) THEN
            IF (JPOS1.EQ.16) IPT5.NUM(10,NEXIS5)=5
            IF (JPOS1.EQ.24) IPT5.NUM(10,NEXIS5)=6
          ENDIF
 189      CONTINUE
        JCOMPT=JCOMPT+JLONG
 190    CONTINUE
 200  CONTINUE
 210  CONTINUE
C
C===================================
C B) Relations dues aux aretes (2D) !
C===================================
C On cree un maillage IPT6 contenant tous les noeuds soumis a des
C relations
C

c        ILPL=LPL(IPT4.ITYPEL)
c        ILPT=LPT(IPT4.ITYPEL)
c        DO 317 J=1,IPT4.NUM(/2)
c        DO 317 K=1,ILPL*2-1,2
c          NPTA=IPT4.NUM(KSEGM(ILPT+K-1),J)
c          NPTB=IPT4.NUM(KSEGM(ILPT+K),J)
c
c          if ((NPTB.eq.129).and.(NPTA.eq.9087)) then 
c            write(*,*) '!ici!',NBPT0 
c          endif
c          IF((NPTA.GT.NBPT0).OR.(NPTB.GT.NBPT0)) THEN
c             GOTO 317
c          ENDIF
c          NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
c          NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
c          if ((NMIN.eq.235).and.(NMAX.eq.236)) then
c              write(*,*) 'nmin', nmin, 'nmax', nmax
c              write(*,*) 'NPTA',NPTA, 'NPTB', NPTB
c          endif
c          NEXIST=0
c          DO 316 I=1,MAX(1,KARPOS(NMIN))
c            IF (KARETE(NMIN,I).EQ.NMAX) THEN
c              KARET2(NMIN,I)=KARET2(NMIN,I)+1
c
c              IF (KMILIE(NMIN,I).gt.0) then
c prise en compte des arretes multiples sur le nouveau maillage
c dans l'exemple ci dessous où les noeuds C et D sont des hanging nodes
c L'arette [C B] n'a pas encore été prise en compte dans karet2 c'est 
c ce qui est fait ici.
c
c |                   |
c |                   |
c |A        C    D    |B
c x---------X----X----X
c |         |    |    |
c |         |    |    |

c1/ Test sur la premiere arete : NMIN-KMILIE(NMIN,I)
c                NMILI = ICPR(KMILIE(NMIN,I),1)
c                NEXIS2=0
c                NMIN2=MIN(NMIN,NMILI)
c                NMAX2=MAX(NMIN,NMILI)
c                DO 318 I2=1,KMILIE(/2)
c                  IF (KMILIE(NMIN2,I2).GT.0) THEN
c                    KARET2(NMIN2,I2)=KARET2(NMIN2,I2)+1
c                  ENDIF
c 318            CONTINUE
C
C 2/ Test sur la deuxieme arete : NMAX-NPTC
c                NEXIS2=0
c                NMIN2=MIN(NMAX,NMILI)
c                NMAX2=MAX(NMAX,NMILI)
c                DO 319 I2=1,KMILIE(/2)
c                  IF (KMILIE(NMIN2,I2).GT.0) THEN
c                    KARET2(NMIN2,I2)=KARET2(NMIN2,I2)+1
c                  ENDIF
c 319            CONTINUE
c              ENDIF
c            ENDIF
c 316      CONTINUE
c 317    CONTINUE

C 1/ Comptage du nombre de noeuds soumis a des relations
      NBELEM=0
      DO 1003 J=1,KMILIE(/2)
      DO 27 I=1,KMILIE(/1)

c        IF (KMILIE(I,J).eq.10330) write(*,*) '!kmilie', I, J, '=10330!'
c        IF (KMILIE(I,J).eq.10330) write(*,*) 'KARET2(I,J)', KARET2(I,J)
c        IF (KMILIE(I,J).eq.10341) write(*,*) '!kmilie', I, J, '=10341!'
c        IF (KMILIE(I,J).eq.10341) write(*,*) 'KARET2(I,J)', KARET2(I,J)
        IF (KARET2(I,J).EQ.0) GOTO 27
        IF (KMILIE(I,J).GT.0) NBELEM=NBELEM+1
        
  27  CONTINUE
 1003 CONTINUE
c      write (*,*) 'NBELEM', NBELEM
C
C 2/ Creation de IPT6
      IPT6=0
      IF (NBELEM.EQ.0) GOTO 999
      NBNN=5
      NBREF=0
      NBSOUS=0
      SEGINI IPT6
      IPT6.ITYPEL=48
C
C 3/ Renseignement des noeuds support des relations
      DO 1004 J=1,KMILIE(/2)
      DO 28 I=1,KMILIE(/1)
        IF (KARET2(I,J).EQ.0) GOTO 28
        IF (KMILIE(I,J).GT.0) THEN
          NBREF=NBREF+1
          IPT6.NUM(1,NBREF)=KMILIE(I,J)
        ENDIF
   28  CONTINUE
 1004 CONTINUE
C       WRITE (*,*) 'nombre de noeuds supports de rela seg', NBREF
C
C 4/ Recherche des noeuds formant les relations
      DO 24 IARAF=1,MELVA2.VELCHE(/2)
        IF (MELVA2.VELCHE(1,IARAF).NE.1) GOTO 24
        JCOMPT=0
        JPOS5=LPOS5(IPT2.ITYPEL)
        DO 23 I=1,NBRAF(IPT2.ITYPEL)
          JPOS1=LPOS1(1,I+LPOS2(IPT2.ITYPEL)-1)
          JLONG=LPOS1(2,I+LPOS2(IPT2.ITYPEL)-1)
          JLISN=LPOS3(IPT2.ITYPEL)+JCOMPT
          LTYPNO=JTYPNO(JPOS5-1+I)
          IF (LTYPNO.NE.0) GOTO 22
          NPTA=IPT2.NUM(LISNOE(JLISN),IARAF)
          NPTB=IPT2.NUM(LISNOE(JLISN+1),IARAF)
          NMIN=MIN(ICPR(NPTA,1),ICPR(NPTB,1))
          NMAX=MAX(ICPR(NPTA,1),ICPR(NPTB,1))
          DO 29 K=1,MAX(1,KARPOS(NMIN))
            IF (KARETE(NMIN,K).EQ.NMAX) NEXIST=K
  29      CONTINUE
          IF (KARET2(NMIN,NEXIST).EQ.0) GOTO 22
          IF (KMILIE(NMIN,NEXIST).EQ.0) GOTO 22
          NEXIS5=0
          DO 20 MM=1,IPT6.NUM(/2)
            INOEGL=KMILIE(NMIN,NEXIST)
            INRELA=IPT6.NUM(1,MM)
            IF (INOEGL.EQ.INRELA) NEXIS5=MM
  20      CONTINUE
C
C 5/ Renseignement des noeuds formant les relations
          DO 21 IHH=1,JLONG
            IPT6.NUM(1+IHH,NEXIS5)=IPT2.NUM(LISNOE(JLISN-1+IHH),IARAF)
  21      CONTINUE
          IF (JLONG.EQ.2) IPT6.NUM(5,NEXIS5)=1
          IF (JLONG.EQ.3) IPT6.NUM(5,NEXIS5)=2
  22      CONTINUE
        JCOMPT=JCOMPT+JLONG
  23    CONTINUE
  24  CONTINUE

C 444 CONTINUE

C
C============================================
C C) Creation du maillage de relations final !
C============================================
      IF (IPT5.EQ.0) THEN
        IPT7=IPT6
        GOTO 999
      ENDIF
      NBELEM=IPT5.NUM(/2)+IPT6.NUM(/2)
C      NBNN=MAX(IPT5.NUM(/1),IPT6.NUM(/1))
      NBNN=10
      NBREF=0
      NBSOUS=0
      SEGINI IPT7
      IPT7.ITYPEL=48
      DO 1005 NEO=1,IPT5.NUM(/2)
      DO 42 MOR=1,10
        IPT7.NUM(MOR,NEO)=IPT5.NUM(MOR,NEO)
        IPT7.ICOLOR(NEO)=IPT5.ICOLOR(NEO)
  42  CONTINUE
 1005 CONTINUE
      NN5 = IPT5.NUM(/2)
      DO 1006 NEO=IPT5.NUM(/2)+1,IPT6.NUM(/2)+IPT5.NUM(/2)
      DO 43 MOR=1,IPT6.NUM(/1)
        IF (MOR.LT.IPT6.NUM(/1)) THEN
          IPT7.NUM(MOR,NEO)=IPT6.NUM(MOR,NEO-NN5)
        ELSE
          IPT7.NUM(10,NEO)=IPT6.NUM(MOR,NEO-NN5)
        ENDIF
        IPT7.ICOLOR(NEO)=IPT6.ICOLOR(NEO-NN5)
  43  CONTINUE
1006  CONTINUE
      SEGDES IPT5, IPT6

C=======================================================================
C                      Fin du programme
C=======================================================================
 999  CONTINUE
c      write (*,*) 'IPT2', IPT2
c      write (*,*) 'IPT3', IPT3
c      write (*,*) 'IPT4', IPT4
c      write (*,*) 'IPT5', IPT5
c      write (*,*) 'IPT6', IPT6
c      write (*,*) 'IPT7', IPT7
      RETURN
      END








 
 
 
 
 
 
 
 
 
