C DIVS      SOURCE    MAGN      10/05/18    21:15:21     6675           
      SUBROUTINE DIVS(IKAS,MTABX)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C************************************************************************       
C    Operateur KMAC
C
C    OBJET   : Cree un objet de type MATRIK
C
C
C  IKAS=1 KMAB calcul de B     (DIV U)
C  IKAS=2 KMBT calcul de Bt uniquement (GRAD P)
C  IKAS=3 KBBT calcul de B assemblage pour B et Bt
C
C***********************************************************************        
-INC CCOPTIO
-INC CCGEOME
-INC CCREEL
-INC SMCOORD
-INC SIZFFB
-INC SMCHPOI
      POINTEUR IZCH2.MCHPOI,IZCCH2.MPOVAL
      POINTEUR IZDV.MCHPOI,IZDDV.MPOVAL,IZTU1.MPOVAL,TETAN.MPOVAL
      POINTEUR IZTG1.MCHPOI,IZTGG1.MPOVAL,IZBETA.MPOVAL
      POINTEUR IZTGG.MCHPOI,IZANU.MPOVAL
-INC SMMATRIK

-INC SMLENTI
      POINTEUR IZIPAD.MLENTI,MLENTI1.MLENTI,MLENTI2.MLENTI,IPADU.MLENTI
-INC SMLMOTS
      POINTEUR LINCO.MLMOTS
-INC SMELEME
      POINTEUR MELEMZ.MELEME,MELEMB.MELEME,MELSTB.MELEME
      POINTEUR MELEM1.MELEME,MELES1.MELEME,MCTREI.MELEME
      POINTEUR IGEOM.MELEME,MELEMA.MELEME
      POINTEUR IZLEMC.MELEME,MELEMS.MELEME,MELEMC.MELEME
      POINTEUR MELEMI.MELEME,MELEMP.MELEME
      POINTEUR MACRO1.MELEME

      CHARACTER*8 TYPE,TYPC,NOMZ,NOMP,NOMD,NOM0
      CHARACTER*8 NOMPP,NOMDD
      CHARACTER*4 NOM
      INTEGER IPAD,IPAD2,IK
      REAL*8 KAUX,TETA1
      DIMENSION IXV(5)
C                                                                               
C*************************************************************************
CKMIC
C     write(6,*)' Operateur KMIC MTABX=',MTABX
C/MELEMS

C- Récupération de la table EQEX (pointeur MTAB1)
C
      CALL LEKTAB(MTABX,'EQEX',MTAB1)
      IF(MTAB1.EQ.0)THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  EQEX  '
         MOTERR( 9:16) = '  EQEX  '
         MOTERR(17:24) = '  KIZX  '
         CALL ERREUR(786)
         RETURN
      ENDIF
      CALL ACME(MTAB1,'NAVISTOK',NASTOK)
      IF(NASTOK.EQ.0)THEN
      CALL ZKMIC(IKAS,MTABX,MTAB1)
      RETURN
      ENDIF
C
C- Récupération de la table INCO (pointeur KINC)
C
      CALL LEKTAB(MTAB1,'INCO',KINC)
      IF(KINC.EQ.0)THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  INCO  '
         MOTERR( 9:16) = '  INCO  '
         MOTERR(17:24) = '  EQEX  '
         CALL ERREUR(786)
         RETURN
      ENDIF

C*************************************************************************
C OPTIONS
C KIMPL = 0 -> EXPL  1 -> IMPL  2 -> CN
C KFORM = 0 -> EFSI  1 -> EF    2 -> VF  3 -> EFMC
C KPRE=3 pression P0   KPRE=4 pression P1  KPRE=2 cas macro 1ère génération
C KPRE=5 pression MSOMMET

      IAXI=0
      IK=0
      IF(IFOMOD.EQ.0)IAXI=2
C
C- Récupération de la table des options KOPT (pointeur KOPTI)
C
      CALL LEKTAB(MTABX,'KOPT',KOPTI)
      IF (KOPTI.EQ.0) THEN
C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
         MOTERR( 1: 8) = '  KOPT  '
         MOTERR( 9:16) = '  KOPT  '
         MOTERR(17:24) = '  KIZX  '
         CALL ERREUR(786)
         RETURN
      ENDIF

      CALL ACME(KOPTI,'KIMPL',KIMPL)
      CALL ACME(KOPTI,'KPOIN',KPRE)
      CALL ACMF(KOPTI,'AIMPL',TETA1)
      CALL ACME(KOPTI,'KFORM',KFORM)
C     gbm rajoute option décentrement 18/10/99
      CALL ACME(KOPTI,'IDCEN',IDCEN)
C     gbm rajoute IKOMP = 2 pour l'opérateur (div) et
C     on ne modifie pas IKOMP pour grad
C     pour appel ultérieur à SUPGEF. 
C     IKOMP = 2 veut dire conservatif mais on a div(flux)
C     au lieu de div(ro*U).
      if (IKAS .EQ. 1 .OR. IKAS .EQ. 3) THEN
         IKOMP = 2
      endif
* GBM modif 22/11/99
      CALL ACMF(KOPTI,'CMD',CMD)


      IF (IERR.NE.0) RETURN
C     write(6,*)' Apres les options '
C*************************************************************************
C
C- Récupération de la table DOMAINE associée au domaine local
C
      CALL ACMM(MTABX,'NOMZONE',NOMZ)
      TYPE=' '
      CALL ACMO(MTABX,'DOMZ',TYPE,MMODEL)
      IF(TYPE.NE.'MMODEL')THEN
C On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40
         MOTERR( 1: 8) = ' MMODEL '
         MOTERR( 8:16) = ' MMODEL '
         MOTERR(17:24) = ' MMODEL '
         MOTERR(25:32) = ' MMODEL '
         MOTERR(33:40) = ' MMODEL '
         CALL ERREUR(471)
         RETURN
      ENDIF

      CALL LEKMOD(MMODEL,MTABZ,INEFMD)

      CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
      CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
      CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
C??   CALL LEKTAB(MTABZ,'MACRO',MACRO)
      MACRO1=0
      IF(INEFMD.EQ.2)CALL LEKTAB(MTABZ,'MACRO1',MACRO1)
C??   IF(INEFMD.EQ.3)CALL LEKTAB(MTABZ,'QUADRATI',MQUAD)
      IF (IERR.NE.0) RETURN
      IF(INEFMD.EQ.4.AND.KPRE.NE.5)THEN
C% Données incompatibles
      CALL ERREUR(21)
      RETURN
      ENDIF

      MELEMI=MELEME
      IF(MACRO1.NE.0.AND.KPRE.NE.2)THEN
      MELEMI=MELEME
      ENDIF

      IF(KPRE.EQ.2.AND.INEFMD.EQ.3)KPRE=3
      IF(INEFMD.EQ.1)KPRE=2

* ON MODIFIE POUR SOMMET -> SOMMET  GBM 18/10/99
       MELEMP=MELEME
       CALL LEKTAB(MTABZ,'CENTRE',MELEMC)

C*************************************************************************
C VERIFICATIONS SUR LES INCONNUES

C     write(6,*)' Verification sur les inconnues '
      TYPE='LISTMOTS'
      CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
      IF (IERR.NE.0) RETURN
      SEGACT LINCO
      NBINC=LINCO.MOTS(/2)
      IF(NBINC.NE.2)THEN
      WRITE(6,*)'Operateur KMAC '
      WRITE(6,*)'Nombre d''inconnues incorrecte : ',NBINC,' On attend 2'
C        Indice %m1:8 : contient plus de %i1 %m9:16
         MOTERR( 1:8) = 'LISTINCO'
         INTERR(1) = 2
         MOTERR(9:16) = ' MOTS   '
         CALL ERREUR(799)
      RETURN
      ENDIF

C     On recupere PHI n et TETA n pour Cranck-Nicholson
      NOMP=LINCO.MOTS(1)
      TYPE=' '
      CALL ACMO(KINC,NOMP,TYPE,MCHPOI)
      IF(TYPE.NE.'CHPOINT ')THEN
         WRITE(6,*)' Opérateur KMAC :'
         WRITE(6,*)' L objet CHPOINT ',NOMP,
     &        ' n existe pas dans la table'
C        Indice %m1:8 : ne contient pas un objet de type %m9:16
         MOTERR( 1: 8) = 'INC '//NOMP
         MOTERR( 9:16) = 'CHPOINT '
         CALL ERREUR(800)
         RETURN
      ELSE
         CALL LICHT(MCHPOI,IZTU1,TYPC,IGEOM0)
      ENDIF
C*************************************************************************
C Le domaine de definition est donne par le SPG  de la premiere inconnue
C Les inconnues suivantes devront posseder ce meme pointeur
C On verifie que les points de la zone sont tous inclus dans ce SPG
C Inconnue Primale

C     write(6,*)' Verification inconnue primale '
      CALL KRIPAD(IGEOM0,IPADU)
      IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN
      MELEMK=MELEMS
      ELSE
      MELEMK=MELEMC
      ENDIF

      CALL VERPAD(IPADU,MELEMK,IRET)
      IF(IRET.NE.0)THEN
      WRITE(6,*)' Opérateur KMAC '
      WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine'
     &        , ' de définition de l''inconnue ',NOMP
      WRITE(6,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0
C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
         MOTERR(1: 8) = 'INC '//NOMP
         MOTERR(9:16) = 'CHPOINT '
         CALL ERREUR(788)
      IPAS=0
      RETURN
      ENDIF


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

      NOMD=LINCO.MOTS(2)
      TYPE=' '
      CALL ACMO(KINC,NOMD,TYPE,MCHPOI)
      IF(TYPE.NE.'CHPOINT ')THEN
         WRITE(6,*)' Opérateur KMAC :'
         WRITE(6,*)' L objet CHPOINT ',NOMD,
     &        ' n existe pas dans la table'
C        Indice %m1:8 : ne contient pas un objet de type %m9:16
         MOTERR( 1: 8) = 'INC '//NOMD
         MOTERR( 9:16) = 'CHPOINT '
         CALL ERREUR(800)
         RETURN
      ELSE
         CALL LICHT(MCHPOI,TETAN,TYPC,IGEOM0)
      ENDIF

      NC=TETAN.VPOCHA(/2)
C*************************************************************************
C Le domaine de definition est donne par le SPG  de la premiere inconnue
C Les inconnues suivantes devront posseder ce meme pointeur
C On verifie que les points de la zone sont tous inclus dans ce SPG
C Inconnue Duale

C     write(6,*)' IGEOM0=',igeom0
      CALL KRIPAD(IGEOM0,MLENTI)
      IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN
C ON MODIFIE POUR CAS MEME ESPACE DIV U  SOMMET -> SOMMET GBM 18/10
      MELEMK=MELEMS
      ELSE
      MELEMK=MELEMK
      ENDIF

C     write(6,*)' Verification inconnue duale ',MELEMK
      CALL VERPAD(MLENTI,MELEMK,IRET)
      IF(IRET.NE.0)THEN
      WRITE(6,*)' Opérateur KMAC '
      WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine'
     &         ,' de définition de l''inconnue ',NOMD
      WRITE(6,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0
C     Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
         MOTERR(1: 8) = 'INC '//NOMD
         MOTERR(9:16) = 'CHPOINT '
         CALL ERREUR(788)
      IPAS=0
      RETURN
      ENDIF

C     GBM rajoute NPTU 18/10
      NPTU = IZTU1.VPOCHA(/1)

C*************************************************************************
C Lecture du ou des coefficients
C Type du coefficient :
C IK1=0 CHPOINT IK1=1 scalaire  IK1=2 vecteur

C     write(6,*)' Verification sur les coefficients '
      CALL ACME(MTABX,'IARG',IARG)
C     write(6,*)' IARG=',iarg

      IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)THEN
      TYPE=' '
      CALL ACMO(MTABZ,'MELSTB',TYPE,MELSTB)
      SEGACT MELSTB
      IF(IDIM.EQ.2)NBELEM=MELSTB.NUM(/2)/4
      IF(IDIM.EQ.3)NBELEM=MELSTB.NUM(/2)/8
      NBNN=MELSTB.NUM(/1)
      NBSOUS=0
      NBREF=0
      SEGINI MELEMA
      MELEMA.ITYPEL=MELSTB.ITYPEL

      NKPE=4
      IF(IDIM.EQ.3)NKPE=8
      do 4878 k=1,nbelem
      mi=(k-1)*NKPE+1
      do 4879 i=1,nbnn
      MELEMA.num(i,k)=melstb.num(i,mi)
 4879 continue
C     write(6,*)k,(MELEMA.num(i,k),i=1,nbnn)
 4878 continue

      TYPE=' '
      CALL ACMO(MTABZ,'MCHPOC',TYPE,MCHPOC)
      TYPE=' '
      CALL ACMO(MTABZ,'CENTRE',TYPE,MCTREI)
      ENDIF


C 1er COEF

      IXV(1)=MELEMC
      IXV(2)=1
      IXV(3)=0
      IXV(4)=MELEMS
      IXV(5)=-MELEMS
      IRET=5
      CALL LEKCOF('Opérateur KMAC :',
     & MTABX,KINC,1,IXV,IZTG1,IZTGG1,NPT1,NC1,IK1,IRET)
      IF(IRET.EQ.0)RETURN
      IF(IK1.GE.4)CALL KRIPAD(MELEMS,MLENTI)


      BETA0=1.D-1
      IF(IARG.EQ.2)THEN
C 2ème COEF
      IXV(1)=0
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur KMAC :',
     & MTABX,KINC,2,IXV,IZTG2,IZBETA,NPT2,NC2,IK2,IRET)
      IF(IRET.EQ.0)RETURN
      BETA0=IZBETA.VPOCHA(1,1)
      ENDIF

C     GBM rajoute 3eme coef mot clef 'SOMM'

C     GBM rajoute 4eme coef viscosité pour SUPGEF

      ANUK=1.D-15
      IF(IARG.EQ.4)THEN
C 4ème COEF
      IXV(1)=0
      IXV(2)=1
      IXV(3)=0
      CALL LEKCOF('Opérateur KMAC :',
     & MTABX,KINC,4,IXV,IZTGG,IZANU,NPT3,NC3,IK3,IRET)
      IF(IRET.EQ.0)RETURN
      ANUK=IZANU.VPOCHA(1,1)
      ENDIF

      NOMP=LINCO.MOTS(1)
      NOMD=LINCO.MOTS(2)(1:4)

      NRIGE=7
      NKID =9
      NKMT =7
      NMATRI=1
      IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)NMATRI=2
      IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.4.AND.IDIM.EQ.2)NMATRI=2
      SEGINI MATRIK

C Fin Stabilisation de toutes sortes
C      write(6,*)'C Fin Stabilisation de toutes sortes'

      NBME=IDIM
      CALL KRIPAD(MELEMI,MLENTI1)
      SEGACT MELEMI
      NBSOUS=MELEMI.LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
      SEGINI IMATRI

      IF(IKAS.EQ.2)THEN
      KSPGD=MELEMS
      KSPGP=MELEMC
      IRIGEL(2,1)=MELEME
      IRIGEL(1,1)=MELEME
      ELSE
      KSPGP=MELEMS
      KSPGD=MELEMS
      IRIGEL(1,1)=MELEME
      IRIGEL(2,1)=MELEME
      ENDIF
      SEGACT MELEMP

C     write(6,*)' ds kmac melemp=',IRIGEL(1,1)
C     write(6,*)' ds kmac melemd=',IRIGEL(2,1)

      IRIGEL(4,1)=IMATRI
* je modifie celui là pour retirer les multiplicateurs
      IF(IKAS.EQ.1)IRIGEL(7,1)= 5
      IF(IKAS.EQ.2)IRIGEL(7,1)=5
      IF(IKAS.EQ.3)IRIGEL(7,1)=5

      K0=0
      DO 11 L=1,NBSOUS
      IPT1=MELEMI
      IF(NBSOUS.NE.1)IPT1=MELEMI.LISOUS(L)
      SEGACT IPT1
      NBEL=IPT1.NUM(/2)

      IF(IKAS.EQ.2)THEN
      MP=IPT1.NUM(/1)
      NP=MELEMP.NUM(/1)
      ELSE
      NP=IPT1.NUM(/1)
      MP=MELEMP.NUM(/1)
      ENDIF

      DO 12 I=1,NBME
      SEGINI IZAFM
C     write(6,*)' ni izafm np=',np,' mp=',mp,' nbel=',nbel,izafm,l,i
      LIZAFM(L,I)=IZAFM
      IF(IKAS.EQ.2)THEN
      WRITE(NOM,FMT='(I1,A3)')I,NOMD(1:3)
      LISDUA(I)=NOM//'    '
      LISPRI(I)=NOMP
      ELSE
      WRITE(NOM,FMT='(I1,A3)')I,NOMP(1:3)
      LISPRI(I)=NOM//'    '
      LISDUA(I)=NOMD
      ENDIF
 12   CONTINUE
      IPM1=LIZAFM(L,1)
      IPM2=LIZAFM(L,2)
      IPM3=LIZAFM(L,2)
      IF(IDIM.EQ.3)IPM3=LIZAFM(L,3)

      IF (IKAS .EQ. 1 .OR. IKAS .EQ. 3) THEN
c     on code div roU -> T, 
         IF(IK1.LT.4)THEN
c        cas où le coef multiplicateur est scalaire, point
c        ou champ au centre.
            CALL KSPRUS
     & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,KO,
     &   ANUK,IDCEN,IKOMP,NPTU,MLENTI,IZTU1,TETAN,NP,NBEL,CMD)
         ELSE
C        MODIFé par GBM, cas coef au sommet.
            CALL KSPRJS
     & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1
     & ,IPADU,IK1)
         ENDIF
      ELSE
c     on code gradT -> U
c     on inverse la place de IZTU1 et de TETAN dans les arguments
c     ca l'utilisateur rentre la pression en premier, IZTU1 est donc
c     une pression (nom mal choisi). GBM 18/10/99
         IF(IK1.LT.4)THEN
c        cas où le coef multiplicateur est scalaire, point
c        ou champ au centre.
            CALL KSPRUS
     & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,KO,
     &   ANUK,IDCEN,IKOMP,NPTU,MLENTI,TETAN,IZTU1,NP,NBEL,CMD)
         ELSE
C        MODIFé par GBM, cas coef au sommet.
            CALL KSPRJS
     & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,
     &  IPADU,IK1)
         ENDIF
      ENDIF
C
      K0=K0+NBEL
      SEGDES IPM1,IPM2,IPM3
      SEGDES IPT1
 11   CONTINUE
      SEGDES MELEMI
      IF(IK1.EQ.4)SEGSUP MLENTI

      CALL ECROBJ('MATRIK',MATRIK)

            NAT=2
            NSOUPO=0
            SEGINI MCHPOI
            JATTRI(1)=2
            CALL ECROBJ('CHPOINT',MCHPOI)

C     write(6,*)' Fin operateur KMIC'
      SEGDES IMATRI,MATRIK
      RETURN
 1001 FORMAT(20(1X,I5))
      END