C KDOM4A    SOURCE    PV        13/04/16    21:15:17     7765           
      SUBROUTINE KDOM4A(MTAB,MELEMQ)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KDOM4A
C                      
C DESCRIPTION       :  Subroutine called by KDOM2A
C                      Axial-symmetric case, TRI7 and QUA8
C                      We compute 
C                      MTAB . 'MAILLAGE'
C                      MTAB . 'CENTRE'
C                      MTAB . 'XCEN2D'
C                      MTAB . 'YCEN2D'
C                      MTAB . 'XXVOLUM'
C                      MTAB . 'XXSURF2D'
C                      and we change the position for the central points
C                      of MELEMQ
C     
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C      
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/LTMF                        
C
C************************************************************************
C
C INPUT/OUTPUT      : MTAB   : domaine table
C                     MELEMQ : QUAF mesh
C************************************************************************
C
C Created the 24/02/04
C
C
C**** Variables de COOPTIO
C
C      INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
C     &     ,IOTER,   IOLEC,  IOIMP,     IOCAR,  IOACQ
C     &     ,IOPER,   IOSGB,  IOGRA,     IOSAU,  IORES
C     &     ,IECHO,   IIMPI,  IOSPI
C     &     ,IDIM,    IFICLE, IPREFI     
CC     &     ,MCOORD
C     &     ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
C     &     ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
C     &     ,NORINC,NORVAL,NORIND,NORVAD
C     &     ,NUCROU, IPSAUV, IREFOR, ISAFOR
CC
C     
      IMPLICIT INTEGER(I-N)
-INC CCOPTIO
-INC SMCOORD
-INC SMELEME
-INC SMLMOTS
-INC SMCHPOI
-INC SMLENTI
      INTEGER MTAB, MELEMQ
     &     ,MCHPSU,MCHPNO,MCHPMR,MCHX2D,MCHY2D
     &     , NBS, NBREF, ITSOUS, NBNN, NBELEM, NBSOUS,NBE0
     &     , JGM, JGN, IGEOM, ISOUS, IELEM
     &     , NN1, NN2, NN3, NNC, NN4
     &     , NFAC, NSOM, NTP, JG, LAST, LASTS
     &     , LISFAC(4), LISSOM(4), NNS, NNOEU
     &     , INOE
     &     , LIFAC(3,4), MCHDIA
      REAL*8 X1,Y1,X2,Y2,X3,Y3,XCEN,YCEN, SURF, VOLU
     &     ,X4,Y4,SURF0,VOLU0,XCEN0,YCEN0,XC2D,YC2D,XC20,YC20
      CHARACTER*8 TYPI
      POINTEUR MELMAI.MELEME, MELCEN.MELEME, MELTFA.MELEME
     &     , MELSOM.MELEME, MELFAC.MELEME, MELFAL.MELEME, MELFAP.MELEME
      POINTEUR MCHVOL.MCHPOI, MPOVOL.MPOVAL
     &     , MCHS2D.MCHPOI, MPOSUR.MPOVAL, MPOX2D.MPOVAL, MPOY2D.MPOVAL
      POINTEUR MLRES.MLENTI, MLRESS.MLENTI, MLEFAC.MLENTI, MLETOF.MLENTI
C
C**** 'MAILLAGE'
C
      MELEME=MELEMQ
      SEGACT MELEME
      SEGINI, MELMAI=MELEME
      NBS=MELEME.LISOUS(/1)
      IF(NBS.EQ.0)NBS=1
      NBREF=0
      IF(NBS .EQ. 1)THEN
         ITSOUS=MELEME.ITYPEL
         IF(ITSOUS .EQ. 7)THEN
C           TRI7 -> TRI3
            MELMAI.ITYPEL=4
            NBNN=3
         ELSEIF(ITSOUS .EQ. 11)THEN
C           QUA9 -> QUA4
            MELMAI.ITYPEL=8
            NBNN=4
         ENDIF
         NBELEM=MELEME.NUM(/2)
         NBSOUS=0
         SEGADJ MELMAI
         NBE0=NBELEM
      ELSE
         NBE0=0
         DO ISOUS=1,NBS,1
            IPT1=MELEME.LISOUS(ISOUS)
            SEGACT IPT1
            ITSOUS=IPT1.ITYPEL
            NBELEM=IPT1.NUM(/2)
            IF(ITSOUS .EQ. 7)THEN
C              TRI7 -> TRI3
               NBNN=3
               MELMAI.ITYPEL=4
            ELSEIF(ITSOUS .EQ. 11)THEN
C              QUA9 -> QUA4
               MELMAI.ITYPEL=8
               NBNN=4
            ENDIF
            NBSOUS=0
            SEGINI IPT2
            MELMAI.LISOUS(ISOUS)=IPT2
            IPT2.ITYPEL=MELMAI.ITYPEL
            MELMAI.ITYPEL=0
            NBE0=NBE0+NBELEM
            IPT1=MELEME.LISOUS(ISOUS)
         ENDDO
      ENDIF
      CALL ECMO(MTAB,'MAILLAGE','MAILLAGE',MELMAI)
C
C**** 'CENTRE'
C
      NBELEM=NBE0
      NBNN=1
      NBSOUS=0
      NBREF=0
      SEGINI MELCEN
      MELCEN.ITYPEL=1
      CALL ECMO(MTAB,'CENTRE','MAILLAGE',MELCEN)
C
C**** 'XXVOLUM', 'XXSURF2D', 'XCEN2D', 'YCEN2D'
C
      TYPI='CENTRE  '
      JGN=4
      JGM=1
      SEGINI MLMOTS
      MLMOTS.MOTS(1)='SCAL'
      CALL KRCHP1(TYPI,MELCEN,MCHVOL,MLMOTS)
      IF(IERR.NE.0) GOTO 9999
      CALL ECMO(MTAB,'XXVOLUM','CHPOINT',MCHVOL)
      IF(IERR.NE.0) GOTO 9999
      CALL LICHT(MCHVOL,MPOVOL,TYPI,IGEOM)
      IF(IERR.NE.0) GOTO 9999
C     SEGACT MPOVOL
C
      CALL KRCHP1(TYPI,MELCEN,MCHS2D,MLMOTS)
      IF(IERR.NE.0) GOTO 9999
      CALL ECMO(MTAB,'XXSUR2D','CHPOINT',MCHS2D)
      IF(IERR.NE.0) GOTO 9999
      CALL LICHT(MCHS2D,MPOSUR,TYPI,IGEOM)
      IF(IERR.NE.0) GOTO 9999
C     SEGACT MPOSUR
C
      CALL KRCHP1(TYPI,MELCEN,MCHX2D,MLMOTS)
      IF(IERR.NE.0) GOTO 9999
      CALL ECMO(MTAB,'XCEN2D','CHPOINT',MCHX2D)
      IF(IERR.NE.0) GOTO 9999
      CALL LICHT(MCHX2D,MPOX2D,TYPI,IGEOM)
      IF(IERR.NE.0) GOTO 9999
C     SEGACT MPOX2D
C
      CALL KRCHP1(TYPI,MELCEN,MCHY2D,MLMOTS)
      IF(IERR.NE.0) GOTO 9999
      CALL ECMO(MTAB,'YCEN2D','CHPOINT',MCHY2D)
      IF(IERR.NE.0) GOTO 9999
      CALL LICHT(MCHY2D,MPOY2D,TYPI,IGEOM)
      IF(IERR.NE.0) GOTO 9999
C     SEGACT MPOY2D
      SEGSUP MLMOTS
C     
C     In KRIPAD
C     SEGDES MELCEN
      SEGACT MELCEN*MOD
C
C**** Filling
C
      NBE0=0
      DO ISOUS=1,NBS,1
         IF(NBS.EQ.1)THEN
            IPT1=MELEME
            IPT2=MELMAI
         ELSE
            IPT1=MELEME.LISOUS(ISOUS)
            IPT2=MELMAI.LISOUS(ISOUS)
         ENDIF
         NBELEM=IPT1.NUM(/2)
         ITSOUS=IPT1.ITYPEL
         IF(ITSOUS .EQ. 7)THEN
C           TRI
            DO IELEM=1,NBELEM,1
               NBE0=NBE0+1
               NN1=IPT1.NUM(1,IELEM)
               NN2=IPT1.NUM(3,IELEM)
               NN3=IPT1.NUM(5,IELEM)
               NNC=IPT1.NUM(7,IELEM)
C
               IPT2.NUM(1,IELEM)=NN1
               IPT2.NUM(2,IELEM)=NN2
               IPT2.NUM(3,IELEM)=NN3
               IPT2.ICOLOR(IELEM)=IPT1.ICOLOR(IELEM)
               MELCEN.NUM(1,NBE0)=NNC
               MELCEN.ICOLOR(NBE0)=IPT1.ICOLOR(IELEM)
C
C************* We compute the position of the center
C              the 'XXVOLUM'
C              the 'XXSUR2D'
C
               X1=XCOOR((NN1-1)*(IDIM+1)+1)
               Y1=XCOOR((NN1-1)*(IDIM+1)+2)
               X2=XCOOR((NN2-1)*(IDIM+1)+1)
               Y2=XCOOR((NN2-1)*(IDIM+1)+2)
               X3=XCOOR((NN3-1)*(IDIM+1)+1)
               Y3=XCOOR((NN3-1)*(IDIM+1)+2)
               CALL KDOM4B(X1,Y1,X2,Y2,X3,Y3,VOLU,SURF
     &              ,XCEN,YCEN,XC2D,YC2D)
C
               MPOVOL.VPOCHA(NBE0,1)=VOLU
               MPOSUR.VPOCHA(NBE0,1)=SURF
               MPOX2D.VPOCHA(NBE0,1)=XC2D
               MPOY2D.VPOCHA(NBE0,1)=YC2D
               XCOOR((NNC-1)*(IDIM+1)+1)=XCEN
               XCOOR((NNC-1)*(IDIM+1)+2)=YCEN
            ENDDO
         ELSEIF(ITSOUS .EQ. 11)THEN
C          QUA
            DO IELEM=1,NBELEM,1
               NBE0=NBE0+1
               NN1=IPT1.NUM(1,IELEM)
               NN2=IPT1.NUM(3,IELEM)
               NN3=IPT1.NUM(5,IELEM)
               NN4=IPT1.NUM(7,IELEM)
               NNC=IPT1.NUM(9,IELEM)
C
               IPT2.NUM(1,IELEM)=NN1
               IPT2.NUM(2,IELEM)=NN2
               IPT2.NUM(3,IELEM)=NN3
               IPT2.NUM(4,IELEM)=NN4
               IPT2.ICOLOR(IELEM)=IPT1.ICOLOR(IELEM)
               MELCEN.NUM(1,NBE0)=NNC
               MELCEN.ICOLOR(NBE0)=IPT1.ICOLOR(IELEM)
C
C************* We compute the position of the center
C              the 'XXVOLUM'
C              the 'XXSUR2D'
C
               X1=XCOOR((NN1-1)*(IDIM+1)+1)
               Y1=XCOOR((NN1-1)*(IDIM+1)+2)
               X2=XCOOR((NN2-1)*(IDIM+1)+1)
               Y2=XCOOR((NN2-1)*(IDIM+1)+2)
               X3=XCOOR((NN3-1)*(IDIM+1)+1)
               Y3=XCOOR((NN3-1)*(IDIM+1)+2)
               X4=XCOOR((NN4-1)*(IDIM+1)+1)
               Y4=XCOOR((NN4-1)*(IDIM+1)+2)
               CALL KDOM4B(X1,Y1,X2,Y2,X4,Y4,VOLU0,SURF0,XCEN0,YCEN0
     $              ,XC20,YC20)
               CALL KDOM4B(X2,Y2,X3,Y3,X4,Y4,VOLU,SURF,XCEN,YCEN
     $              ,XC2D,YC2D)
C
               MPOVOL.VPOCHA(NBE0,1)=VOLU+VOLU0
               MPOSUR.VPOCHA(NBE0,1)=SURF+SURF0
               MPOX2D.VPOCHA(NBE0,1)=((XC20*SURF0)+(XC2D*SURF))
     $              /(SURF+SURF0)
               MPOY2D.VPOCHA(NBE0,1)=((YC20*SURF0)+(YC2D*SURF))
     $              /(SURF+SURF0)
               XCOOR((NNC-1)*(IDIM+1)+1)=((XCEN0*VOLU0)+(XCEN*VOLU))
     $              /(VOLU+VOLU0)
               XCOOR((NNC-1)*(IDIM+1)+2)=((YCEN0*VOLU0)+(YCEN*VOLU))
     $              /(VOLU+VOLU0)
            ENDDO
         ENDIF
         SEGDES IPT2
      ENDDO
C
      IF(NBS.NE.1)THEN
         SEGDES MELMAI
      ENDIF
C
      SEGDES MPOSUR
      SEGDES MPOVOL
      SEGDES MELCEN
C
C     MELEME et ses "fils" sont toujours actifs
C
C**** We create ELTFA, FACE and SOMMET
C     N.B. The position of the noeud belonging to the FACE
C          is not correct
C
      SEGINI, MELTFA=MELEME
      NBREF=0
      IF(NBS .EQ. 1)THEN
         ITSOUS=MELEME.ITYPEL
         IF(ITSOUS .EQ. 7)THEN
C           TRI3
            NBNN=3
            MELTFA.ITYPEL=4
C           ELTFA TRI3
         ELSEIF(ITSOUS .EQ. 11)THEN
C           QUA4
            NBNN=4
            MELTFA.ITYPEL=8
C           ELTFA QUA4
         ENDIF
         NBELEM=MELEME.NUM(/2)
         NBSOUS=0
         SEGADJ MELTFA
      ELSE
         DO ISOUS=1,NBS,1
            IPT1=MELEME.LISOUS(ISOUS)
            NBELEM=IPT1.NUM(/2)
            ITSOUS=IPT1.ITYPEL
            IF(ITSOUS .EQ. 7)THEN
C              TRI3
               NBNN=3
               MELTFA.ITYPEL=4
            ELSEIF(ITSOUS .EQ. 11)THEN
C              QUA4
               NBNN=4
               MELTFA.ITYPEL=8
C              ELTFA QUA4
            ENDIF
            NBSOUS=0
            SEGINI IPT2
            MELTFA.LISOUS(ISOUS)=IPT2
C              
            IPT2.ITYPEL=MELTFA.ITYPEL
            MELTFA.ITYPEL=0
         ENDDO
      ENDIF
C
C**** We fill ELTFA
C     We also count:
C     NFAC  = number of non-triangular faces
C     NSOM  = number of SOMMET
C
      NTP=MCOORD.XCOOR(/1)/(IDIM+1)
      JG=NTP
      NFAC=0
      NSOM=0

      LAST=-1
      SEGINI MLRES
      LASTS=-1
      SEGINI MLRESS
C     LAST+MLRES = chaining list to find the faces
C     LASTS+MLRESS = chaining list to find the sommet

      DO ISOUS=1,NBS,1
         IF(NBS.EQ.1) THEN
            IPT1=MELEME
            IPT2=MELTFA
         ELSE
            IPT1=MELEME.LISOUS(ISOUS)
            IPT2=MELTFA.LISOUS(ISOUS)
         ENDIF
C         
         NBELEM=IPT1.NUM(/2)
         ITSOUS=IPT1.ITYPEL
         IF(ITSOUS .EQ. 7)THEN
C          TRI (2D)
            LISFAC(1)=2
            LISFAC(2)=4
            LISFAC(3)=6
            LISSOM(1)=1
            LISSOM(2)=3
            LISSOM(3)=5
            NNS=3
            NNOEU=3
         ELSEIF(ITSOUS .EQ. 11)THEN
C           QUA (2D)
            LISFAC(1)=2
            LISFAC(2)=4
            LISFAC(3)=6
            LISFAC(4)=8
            LISSOM(1)=1
            LISSOM(2)=3
            LISSOM(3)=5
            LISSOM(4)=7
            NNS=4
            NNOEU=4
         ENDIF
C
         DO IELEM=1,NBELEM,1
            DO INOE=1,NNOEU,1
               NN1=IPT1.NUM(LISFAC(INOE),IELEM)
               IPT2.NUM(INOE,IELEM)=NN1
               IF(MLRES.LECT(NN1) .EQ. 0)THEN
                  NFAC=NFAC+1
                  MLRES.LECT(NN1)=LAST
                  LAST=NN1
               ENDIF
            ENDDO
            DO INOE=1,NNS,1
               NN1=IPT1.NUM(LISSOM(INOE),IELEM)
               IF(MLRESS.LECT(NN1) .EQ. 0)THEN
                  NSOM=NSOM+1
                  MLRESS.LECT(NN1)=LASTS
                  LASTS=NN1
               ENDIF
            ENDDO
         ENDDO
C
         SEGDES IPT2
C     
      ENDDO
      IF(NBS. NE. 1)  SEGDES MELTFA
      CALL ECMO(MTAB,'ELTFA','MAILLAGE',MELTFA)
      IF(IERR .NE. 0) GOTO 9999
C
C******** Creation of SOMMET
C
      NBELEM=NSOM
      NBNN=1
      NBSOUS=0
      NBREF=0
      SEGINI MELSOM
      MELSOM.ITYPEL=1
      DO IELEM=1,NSOM,1
         MELSOM.NUM(1,IELEM)=LASTS
         LASTS=MLRESS.LECT(LASTS)
      ENDDO
      IF(LASTS .NE. -1)THEN
         WRITE(IOIMP,*) 'Subroutine kdom10.eso'
         CALL ERREUR(5)
      ENDIF
      call crech1(melsom,1)
      SEGDES MELSOM
      CALL ECMO(MTAB,'SOMMET','MAILLAGE',MELSOM)
      SEGSUP MLRESS
C     
C**** Creation of FACE
C
      NBELEM=NFAC
      NBNN=1
      NBSOUS=0
      NBREF=0
      SEGINI MELFAC
      MELFAC.ITYPEL=1
      DO IELEM=1,NFAC,1
         MELFAC.NUM(1,IELEM)=LAST
         LAST=MLRES.LECT(LAST)
      ENDDO
      IF(LAST .NE. -1)THEN
         WRITE(IOIMP,*) 'Subroutine kdom10.eso'
         CALL ERREUR(5)
      ENDIF
      SEGDES MELFAC
      SEGSUP MLRES
      CALL ECMO(MTAB,'FACE','MAILLAGE',MELFAC)
C
C******* Creation of FACEL and FACEP
C
      CALL KRIPAD(MELFAC,MLEFAC)
C     SEGINI MLEFAC
      JG=NFAC
      SEGINI MLETOF
C                          
C     MLETOF.LECT(I1) = how many times has the i-th face of MELFAP 
C                       already been touched?
C
      NBELEM=NFAC
      NBNN=3
      NBSOUS=0
      NBREF=0
      SEGINI MELFAL
      MELFAL.ITYPEL=3
C
C     FACEP is a SEG3
C
      NBELEM=NFAC
      NBNN=3
      NBSOUS=0
      NBREF=0
      SEGINI MELFAP
      MELFAP.ITYPEL=3
C
      DO ISOUS=1,NBS,1
C
C********** Loop on the elementary mesh of the QUAF
C
         IF(NBS.EQ.1) THEN
            IPT1=MELEME
         ELSE
            IPT1=MELEME.LISOUS(ISOUS)
         ENDIF
C        
         ITSOUS=IPT1.ITYPEL
         IF(ITSOUS .EQ. 7)THEN
C           TRI (2D)
            LIFAC(1,1)=2
            LIFAC(2,1)=1
            LIFAC(3,1)=3
            LIFAC(1,2)=4
            LIFAC(2,2)=3
            LIFAC(3,2)=5
            LIFAC(1,3)=6
            LIFAC(2,3)=5
            LIFAC(3,3)=1
C           Here we put the center point in LISSOM
            LISSOM(1)=7
            NNOEU=3
C               
         ELSEIF(ITSOUS .EQ. 11)THEN
C          QUA (2D)
            LIFAC(1,1)=2
            LIFAC(2,1)=1
            LIFAC(3,1)=3
            LIFAC(1,2)=4
            LIFAC(2,2)=3
            LIFAC(3,2)=5
            LIFAC(1,3)=6
            LIFAC(2,3)=5
            LIFAC(3,3)=7
            LIFAC(1,4)=8
            LIFAC(2,4)=7
            LIFAC(3,4)=1
            LISSOM(1)=9
            NNOEU=4
         ENDIF
C
         NBELEM=IPT1.NUM(/2)
         DO IELEM=1,NBELEM,1
C           NNOEU = number of quagrangular elements
            DO INOE=1,NNOEU,1
C              NN1 is the global number of the face
C              NN2 is the local number of the face in the MELEME
C                 'FACE'
C
               NN1=IPT1.NUM(LIFAC(1,INOE),IELEM)
               NN2=MLEFAC.LECT(NN1)
               IF(MLETOF.LECT(NN2).EQ.0)THEN
C
C                    MLETOF.LECT(NN2) = how many times the face NN2 has
C                    been touched?
C
                  MLETOF.LECT(NN2)=1
                  MELFAL.NUM(2,NN2)=NN1
                  MELFAL.NUM(1,NN2)=IPT1.NUM(LISSOM(1),IELEM)
                  MELFAL.NUM(3,NN2)=IPT1.NUM(LISSOM(1),IELEM)
                  MELFAP.NUM(1,NN2)=IPT1.NUM(LIFAC(2,INOE),IELEM)
                  MELFAP.NUM(2,NN2)=IPT1.NUM(LIFAC(3,INOE),IELEM)
                  MELFAP.NUM(3,NN2)=NN1
               ELSEIF(MLETOF.LECT(NN2).EQ.1)THEN
                  MLETOF.LECT(NN2)=2
                  MELFAL.NUM(3,NN2)=IPT1.NUM(LISSOM(1),IELEM)
               ELSE
                  WRITE(IOIMP,*) 'Subroutine kdom10.eso'
                  CALL ERREUR(5)
               ENDIF
            ENDDO
         ENDDO
      ENDDO
      SEGDES MELFAL
      SEGDES MELFAP
      SEGSUP MLETOF
      SEGSUP MLEFAC
      CALL ECMO(MTAB,'FACEL','MAILLAGE',MELFAL)
      IF(IERR .NE. 0) GOTO 9999
      CALL ECMO(MTAB,'FACEP','MAILLAGE',MELFAP)
      IF(IERR .NE. 0) GOTO 9999
      IF(NBS.NE.1)THEN
         SEGDES MELEME
      ENDIF
C     
C**** We have to create 
C     'XXSURFAC'
C     'XXNORMAF' 
C     'MATROT'
C     and to put the face centre in the right position!
C
      CALL KDOM4C(MELFAC,MELFAL,MELFAP,MCHPSU,MCHPNO,MCHPMR)
      IF(IERR.NE.0)GOTO 9999
C
      CALL ECMO(MTAB,'XXSURFAC','CHPOINT',MCHPSU)
      IF(IERR.NE.0) GOTO 9999
      CALL ECMO(MTAB,'XXNORMAF','CHPOINT',MCHPNO)
      IF(IERR.NE.0) GOTO 9999
      CALL ECMO(MTAB,'MATROT','CHPOINT',MCHPMR)
      IF(IERR.NE.0) GOTO 9999
C
C**** Finally, we compute XXDIEMIN
C
      CALL KDOM12(MELTFA,MELCEN,MELFAC,MCHPNO,MCHDIA)
      IF(IERR.NE.0) GOTO 9999
      CALL ECMO(MTAB,'XXDIEMIN','CHPOINT',MCHDIA)
      IF(IERR.NE.0) GOTO 9999
C
 9999 RETURN
      END