Sources Esope | Cast3M

konja7
C KONJA7    SOURCE    OF166741  24/12/13    21:16:35     12097          
      SUBROUTINE KONJA7(ILINC,IRN,IUN,IPN,IGAMN,INORM,ICHPVO
     $     ,ICHPSU,IUINF,IUPRI,MELEMC,MELEFE,MELLIM,IMAT)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KONJA7
C
C DESCRIPTION       :  Voir KON12
C                      Calcul du jacobien du résidu pour la méthode de
C                      AUSM+low mach
C
C                      Cas 3D, gaz "calorically perfect"
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. BECCANTINI, SFME/LTMF
C
C************************************************************************
C
C
C APPELES (Outils
C          CASTEM)  :
C
C APPELES (Calcul)  :
C
C************************************************************************
C
C ENTREES
C
C     ILINC  : liste des inconnues (pointeur d'un objet de type LISTMOTS)
C
C  1) Pointeurs des CHPOINT
C
C     IRN     : CHPOINT CENTRE  contenant la masse volumique ;
C
C     IUN     : CHPOINT CENTRE  contenant la vitesse ;
C
C     IPN     : CHPOINT CENTRE  contenant la pression ;
C
C     IGAMN   : CHPOINT CENTRE  contenant le gamma ;
C
C     INORM   : CHPOINT FACE    contenant les normales aux faces ;
C
C     ICHPOVO : CHPOINT VOLUME  contenant le volume
C
C     ICHPOSU : CHPOINT FACE    contenant la surface des faces
C
C     IUINF   : CHPOINT, one component, "cut-off velocity"
C               0 if no Bas MACH
C
C     IUPRI   : CHPOINT, one component, "minimum reference velocity"
C               0 if no BAs Mach
C
C  2) Pointeurs de MELEME de la table DOMAINE
C
C     MELEMC  : MELEME 'CENTRE' du SPG des CENTRES
C
C     MELEFE  : MELEME 'FACEL' du connectivité Faces -> Elts
C
C     MELLIM  : MELEME SPG des conditions aux bords
C
C SORTIES
C
C     IMAT    : pointeur de la MATRIK du jacobien du residu
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :
C
C************************************************************************
C
C
C N.B.: On suppose qu'on a déjà controllé RO, P > 0
C                                         GAMMA \in (1,3)
C       Si non il faut le faire!!!
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)
      INTEGER ILINC, IRN,IUN,IPN,IGAMN,INORM,ICHPVO,ICHPSU
     &        , IMAT, IGEOMC, IGEOMF, IUINF, IUPRI
     &        , NFAC, NBSOUS, NBREF, NBELEM, NBNN, NRIGE, NMATRI, NKID
     &        , NKMT, NBME, NBEL,  MP, NP
     &        , IFAC, NGCF, NLCF, NGCG, NGCD, NLCG, NLCD
     &        , NLFL
      REAL*8 ROG, PG, UXG, UYG, UZG, GAMG, VOLG
     &       , ROD, PD, UXD, UYD, UZD, VOLD
     &       , SURF
     &       , FUNCEL
      REAL*8 WVEC_L(5), WVEC_R(5), NVECT(3), TVECT1(3),TVECT2(3)
      REAL*8 JTL(5,5), JTR(5,5), V_INF
      REAL*8 C11,C12,C13,C21,C22,C23,C31,C32,C33,DET
      REAL*8 ZC11,ZC12,ZC13,ZC21,ZC22,ZC23,ZC31,ZC32,ZC33
      CHARACTER*8 TYPE
C
C**** LES INCLUDES
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
-INC SMMATRIK
-INC SMELEME
-INC SMLMOTS
-INC SMLENTI
      POINTEUR MPRN.MPOVAL, MPUN.MPOVAL, MPPN.MPOVAL, MPGAMN.MPOVAL,
     &         MPNORM.MPOVAL, MPVOLU.MPOVAL, MPOVSU.MPOVAL
      POINTEUR MELEMC.MELEME, MELEMF.MELEME, MELEFE.MELEME,
     &         MELEDU.MELEME, MELLIM.MELEME
      POINTEUR MLENTC.MLENTI, MLENTF.MLENTI, MLELIM.MLENTI
      POINTEUR RR.IZAFM, RUX.IZAFM, RUY.IZAFM, RUZ.IZAFM, RRET.IZAFM,
     &     UXR.IZAFM, UXUX.IZAFM, UXUY.IZAFM, UXUZ.IZAFM, UXRET.IZAFM,
     &     UYR.IZAFM, UYUX.IZAFM, UYUY.IZAFM, UYUZ.IZAFM, UYRET.IZAFM,
     &     UZR.IZAFM, UZUX.IZAFM, UZUY.IZAFM, UZUZ.IZAFM, UZRET.IZAFM,
     &     RETR.IZAFM, RETUX.IZAFM, RETUY.IZAFM, RETUZ.IZAFM,
     &     RETRET.IZAFM
      POINTEUR MLMINC.MLMOTS, MPUPRI.MPOVAL, MPUINF.MPOVAL
C
C**** KRIPAD pour la correspondance global/local des conditions limits
C
      CALL KRIPAD(MELLIM,MLELIM)
C     SEGACT MELLIM
C
C**** KRIPAD pour la correspondance global/local des centres
C
      CALL KRIPAD(MELEMC,MLENTC)
C
C     SEGACT MLENTC
      SEGACT MELEMC
C
      SEGACT MELEFE
C
      CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
      CALL LICHT(INORM,MPNORM,TYPE,IGEOMF)
      CALL LICHT(ICHPVO,MPVOLU,TYPE,IGEOMC)
C
C**** LICHT active les MPOVALs en *MOD
C
C     i.e.
C
C     SEGACT MPOVSU*MOD
C     SEGACT MPOVNO*MOD
C     SEGACT MPVOLU*MOD
C
      MELEMF = IGEOMF
      CALL KRIPAD(MELEMF,MLENTF)
C
C     SEGACT MLENTF
      SEGACT MELEMF
C
      CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
      CALL LICHT(IPN,MPPN,TYPE,IGEOMC)
      CALL LICHT(IUN,MPUN,TYPE,IGEOMC)
      CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC)
C
C     SEGACT MPRN*MOD
C     SEGACT MPPN*MOD
C     SEGACT MPUN*MOD
C     SEGACT MPGAMN*MOD
C
      NFAC = MELEFE.NUM(/2)
C
C**** Maillage des inconnues primales
C
      NBSOUS = 0
      NBREF = 0
      NBELEM =  NFAC
      NBNN = 2
C
      SEGINI MELEDU
C     MELEPR = MELEDU
C
C**** MELEDU = 'SEG2'
C
      MELEDU.ITYPEL =  2
C
      NRIGE = 7
      NMATRI = 1
      NKID = 9
      NKMT = 7
C
      SEGINI MATRIK
      IMAT = MATRIK
      MATRIK.IRIGEL(1,1) = MELEDU
      MATRIK.IRIGEL(2,1) = MELEDU
C
C**** Matrice non symetrique
C
      MATRIK.IRIGEL(7,1) = 2
C
      NBME = 25
      NBSOUS = 1
      SEGINI IMATRI
      MLMINC = ILINC
      SEGACT MLMINC
      MATRIK.IRIGEL(4,1) = IMATRI
C
      IMATRI.LISPRI(1) = MLMINC.MOTS(1)
      IMATRI.LISPRI(2) = MLMINC.MOTS(2)
      IMATRI.LISPRI(3) = MLMINC.MOTS(3)
      IMATRI.LISPRI(4) = MLMINC.MOTS(4)
      IMATRI.LISPRI(5) = MLMINC.MOTS(5)
C
      IMATRI.LISPRI(6)  = MLMINC.MOTS(1)
      IMATRI.LISPRI(7)  = MLMINC.MOTS(2)
      IMATRI.LISPRI(8)  = MLMINC.MOTS(3)
      IMATRI.LISPRI(9)  = MLMINC.MOTS(4)
      IMATRI.LISPRI(10) = MLMINC.MOTS(5)
C
      IMATRI.LISPRI(11) = MLMINC.MOTS(1)
      IMATRI.LISPRI(12) = MLMINC.MOTS(2)
      IMATRI.LISPRI(13) = MLMINC.MOTS(3)
      IMATRI.LISPRI(14) = MLMINC.MOTS(4)
      IMATRI.LISPRI(15) = MLMINC.MOTS(5)
C
      IMATRI.LISPRI(16) = MLMINC.MOTS(1)
      IMATRI.LISPRI(17) = MLMINC.MOTS(2)
      IMATRI.LISPRI(18) = MLMINC.MOTS(3)
      IMATRI.LISPRI(19) = MLMINC.MOTS(4)
      IMATRI.LISPRI(20) = MLMINC.MOTS(5)
C
      IMATRI.LISPRI(21) = MLMINC.MOTS(1)
      IMATRI.LISPRI(22) = MLMINC.MOTS(2)
      IMATRI.LISPRI(23) = MLMINC.MOTS(3)
      IMATRI.LISPRI(24) = MLMINC.MOTS(4)
      IMATRI.LISPRI(25) = MLMINC.MOTS(5)
C
      IMATRI.LISDUA(1) = MLMINC.MOTS(1)
      IMATRI.LISDUA(2) = MLMINC.MOTS(1)
      IMATRI.LISDUA(3) = MLMINC.MOTS(1)
      IMATRI.LISDUA(4) = MLMINC.MOTS(1)
      IMATRI.LISDUA(5) = MLMINC.MOTS(1)
C
      IMATRI.LISDUA(6)  = MLMINC.MOTS(2)
      IMATRI.LISDUA(7)  = MLMINC.MOTS(2)
      IMATRI.LISDUA(8)  = MLMINC.MOTS(2)
      IMATRI.LISDUA(9)  = MLMINC.MOTS(2)
      IMATRI.LISDUA(10) = MLMINC.MOTS(2)
C
      IMATRI.LISDUA(11) = MLMINC.MOTS(3)
      IMATRI.LISDUA(12) = MLMINC.MOTS(3)
      IMATRI.LISDUA(13) = MLMINC.MOTS(3)
      IMATRI.LISDUA(14) = MLMINC.MOTS(3)
      IMATRI.LISDUA(15) = MLMINC.MOTS(3)
C
      IMATRI.LISDUA(16) = MLMINC.MOTS(4)
      IMATRI.LISDUA(17) = MLMINC.MOTS(4)
      IMATRI.LISDUA(18) = MLMINC.MOTS(4)
      IMATRI.LISDUA(19) = MLMINC.MOTS(4)
      IMATRI.LISDUA(20) = MLMINC.MOTS(4)
C
      IMATRI.LISDUA(21) = MLMINC.MOTS(5)
      IMATRI.LISDUA(22) = MLMINC.MOTS(5)
      IMATRI.LISDUA(23) = MLMINC.MOTS(5)
      IMATRI.LISDUA(24) = MLMINC.MOTS(5)
      IMATRI.LISDUA(25) = MLMINC.MOTS(5)
C
      NBEL = NBELEM
      NBSOUS = 1
      NP = 2
      MP = 2
      SEGINI   RR  , RUX , RUY , RUZ, RRET ,
     &         UXR , UXUX , UXUY , UXUZ, UXRET ,
     &         UYR , UYUX , UYUY , UYUZ, UYRET ,
     &         UZR , UZUX , UZUY , UZUZ, UZRET ,
     &         RETR , RETUX , RETUY , RETUZ, RETRET
C
C**** Duale = IMATRI.LISDUA(1) = 'RN'
C     Primale = IMATRI.LISPRI(1) = 'RN'
C     -> IMATRI.LIZAFM(1,1) = RR
C
C     Duale = IMATRI.LISDUA(2) = 'RN'
C     Primale = IMATRI.LISPRI(1) = 'RUXN'
C     -> IMATRI.LIZAFM(1,2) = RUX
C     ...
C
      IMATRI.LIZAFM(1,1) = RR
      IMATRI.LIZAFM(1,2) = RUX
      IMATRI.LIZAFM(1,3) = RUY
      IMATRI.LIZAFM(1,4) = RUZ
      IMATRI.LIZAFM(1,5) = RRET
C
      IMATRI.LIZAFM(1,6) = UXR
      IMATRI.LIZAFM(1,7) = UXUX
      IMATRI.LIZAFM(1,8) = UXUY
      IMATRI.LIZAFM(1,9) = UXUZ
      IMATRI.LIZAFM(1,10) = UXRET
C
      IMATRI.LIZAFM(1,11) = UYR
      IMATRI.LIZAFM(1,12) = UYUX
      IMATRI.LIZAFM(1,13) = UYUY
      IMATRI.LIZAFM(1,14) = UYUZ
      IMATRI.LIZAFM(1,15) = UYRET
C
      IMATRI.LIZAFM(1,16) = UZR
      IMATRI.LIZAFM(1,17) = UZUX
      IMATRI.LIZAFM(1,18) = UZUY
      IMATRI.LIZAFM(1,19) = UZUZ
      IMATRI.LIZAFM(1,20) = UZRET
C
      IMATRI.LIZAFM(1,21) = RETR
      IMATRI.LIZAFM(1,22) = RETUX
      IMATRI.LIZAFM(1,23) = RETUY
      IMATRI.LIZAFM(1,24) = RETUZ
      IMATRI.LIZAFM(1,25) = RETRET
C**************************************************************
C    Bas Mach
C**************************************************************
      CALL LICHT(IUPRI,MPUPRI,TYPE,IGEOMC)
      CALL LICHT(IUINF,MPUINF,TYPE,IGEOMC)
C**************************************************************
      DO IFAC = 1, NFAC, 1
         NGCF = MELEFE.NUM(2,IFAC)
         NLCF = MLENTF.LECT(NGCF)
         IF(NLCF .NE. IFAC)THEN
            WRITE(IOIMP,*) 'Il ne faut pas jouer avec la table domaine'
            CALL ERREUR(5)
            GOTO 9999
         ENDIF
         NLFL = MLELIM.LECT(NGCF)
         NGCG = MELEFE.NUM(1,IFAC)
         NGCD = MELEFE.NUM(3,IFAC)
         IF(NLFL .NE. 0)THEN
C
C********** The point belongs on BC -> No contribution to jacobian!
C
            MELEDU.NUM(1,IFAC) = NGCG
            MELEDU.NUM(2,IFAC) = NGCD
         ELSEIF(NGCG .NE. NGCD)THEN
C
C********** Les MELEMEs
C
            MELEDU.NUM(1,IFAC) = NGCG
            MELEDU.NUM(2,IFAC) = NGCD
C
C********** Les etats G et D
C
            NLCG = MLENTC.LECT(NGCG)
            NLCD = MLENTC.LECT(NGCD)
C
            ROG = MPRN.VPOCHA(NLCG,1)
            PG = MPPN.VPOCHA(NLCG,1)
            UXG = MPUN.VPOCHA(NLCG,1)
            UYG = MPUN.VPOCHA(NLCG,2)
            UZG = MPUN.VPOCHA(NLCG,3)
            GAMG = MPGAMN.VPOCHA(NLCG,1)
            VOLG = MPVOLU.VPOCHA(NLCG,1)
C-----------------------------------------
            WVEC_L(1)=ROG
            WVEC_L(2)=UXG
            WVEC_L(3)=UYG
            WVEC_L(4)=UZG
            WVEC_L(5)=PG
C-----------------------------------------
            ROD = MPRN.VPOCHA(NLCD,1)
            PD = MPPN.VPOCHA(NLCD,1)
            UXD = MPUN.VPOCHA(NLCD,1)
            UYD = MPUN.VPOCHA(NLCD,2)
            UZD = MPUN.VPOCHA(NLCD,3)
            VOLD = MPVOLU.VPOCHA(NLCD,1)
c-----------------------------------------------------
            WVEC_R(1)=ROD
            WVEC_R(2)=UXD
            WVEC_R(3)=UYD
            WVEC_R(4)=UZD
            WVEC_R(5)=PD
C-----------------------------------------------------
C********** La normale G->D
C           La tangente
C-----------------------------------------------------
c            SURF = MPOVSU.VPOCHA(NLCF,1)
c            CNX = MPNORM.VPOCHA(NLCF,7)
c            CNY = MPNORM.VPOCHA(NLCF,8)
c            CNZ = MPNORM.VPOCHA(NLCF,9)
cC
cC********** Cosinus directeurs de tangente 1
cC
c            CT1X = MPNORM.VPOCHA(NLCF,1)
c            CT1Y = MPNORM.VPOCHA(NLCF,2)
c            CT1Z = MPNORM.VPOCHA(NLCF,3)
cC
cC********** Cosinus directeurs de tangente 2
cC
c            CT2X = MPNORM.VPOCHA(NLCF,4)
c            CT2Y = MPNORM.VPOCHA(NLCF,5)
c            CT2Z = MPNORM.VPOCHA(NLCF,6)
C-------------------------------------------------------
            SURF = MPOVSU.VPOCHA(NLCF,1)
            NVECT(1) = MPNORM.VPOCHA(NLCF,7)
            NVECT(2) = MPNORM.VPOCHA(NLCF,8)
            NVECT(3) = MPNORM.VPOCHA(NLCF,9)
c-----------------------------------------------
            TVECT1(1) = MPNORM.VPOCHA(NLCF,1)
            TVECT1(2) = MPNORM.VPOCHA(NLCF,2)
            TVECT1(3) = MPNORM.VPOCHA(NLCF,3)
c----------------------------------------------
            TVECT2(1) = MPNORM.VPOCHA(NLCF,4)
            TVECT2(2) = MPNORM.VPOCHA(NLCF,5)
            TVECT2(3) = MPNORM.VPOCHA(NLCF,6)
C-------------------------------------------------------
C********** La contribution de Gauche
C--------------------------------------------------------
            V_INF=MAX(MPUINF.VPOCHA(NLCG,1),
     &              MPUINF.VPOCHA(NLCD,1),
     &              MPUPRI.VPOCHA(NLCG,1),
     &              MPUPRI.VPOCHA(NLCD,1))
            CALL JA3DBM(jtl,jtr,wvec_l,wvec_r,nvect,tvect1,tvect2,
     &                  gamg,v_inf)
C
C
C********** AB.AM(IFAC,IPRIM,IDUAL)
C           A = nom de l'inconnu duale (Ro,rUX,rUY,RET)
C           B = nom de l'inconnu primale (Ro,rUX,rUY,RET)
C           IPRIM = 1, 2  -> G, D
C           IDUAL = 1, 2  -> G, D
C           i.e.
C           A_IDUAL = AB.AM(IFAC,IPRIM,IDUAL) * B_IPRIM + ...
C
C
C********** Dual RN
            FUNCEL = SURF * JTL(1,1)
            RR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RR.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(1,2)
            RUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RUX.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(1,3)
            RUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RUY.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(1,4)
            RUZ.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RUZ.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(1,5)
            RRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RRET.AM(IFAC,1,2) = FUNCEL / VOLD
c----------------------------------------------------
C********** Dual RUXN
C----------------------------------------------------
            FUNCEL = SURF * JTL(2,1)
            UXR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXR.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(2,2)
            UXUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXUX.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(2,3)
            UXUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXUY.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(2,4)
            UXUZ.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXUZ.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(2,5)
            UXRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXRET.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
C********** Dual RUYN
C----------------------------------------------------
            FUNCEL = SURF * JTL(3,1)
            UYR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYR.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(3,2)
            UYUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYUX.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(3,3)
            UYUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYUY.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(3,4)
            UYUZ.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYUZ.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(3,5)
            UYRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYRET.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
C********** Dual RUZN
C----------------------------------------------------
            FUNCEL = SURF * JTL(4,1)
            UZR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UZR.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(4,2)
            UZUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UZUX.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(4,3)
            UZUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UZUY.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(4,4)
            UZUZ.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UZUZ.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
            FUNCEL = SURF * JTL(4,5)
            UZRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UZRET.AM(IFAC,1,2) = FUNCEL / VOLD
C----------------------------------------------------
C********** Dual RETN
C----------------------------------------------------
             FUNCEL = SURF * JTL(5,1)
            RETR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RETR.AM(IFAC,1,2) = FUNCEL / VOLD
C-----------------------------------------------------
            FUNCEL = SURF * JTL(5,2)
            RETUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RETUX.AM(IFAC,1,2) = FUNCEL / VOLD
C-----------------------------------------------------
            FUNCEL = SURF * JTL(5,3)
            RETUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RETUY.AM(IFAC,1,2) = FUNCEL / VOLD
C-----------------------------------------------------
            FUNCEL = SURF * JTL(5,4)
            RETUZ.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RETUZ.AM(IFAC,1,2) = FUNCEL / VOLD
C-----------------------------------------------------
            FUNCEL = SURF * JTL(5,5)
            RETRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RETRET.AM(IFAC,1,2) = FUNCEL / VOLD
C-----------------------------------------------------
C
C********** La contribution de D
C
c            CNX = -1.0D0 * CNX
c            CNY = -1.0D0 * CNY
c            CNZ = -1.0D0 * CNZ
c            CT1X = -1.0D0 * CT1X
c            CT1Y = -1.0D0 * CT1Y
c            CT1Z = -1.0D0 * CT1Z
c            CT2X = -1.0D0 * CT2X
c            CT2Y = -1.0D0 * CT2Y
c            CT2Z = -1.0D0 * CT2Z
C
c            CALL VLHJ4(ROD,UXD,UYD,UZD,PD,RETD,GAMD,CNX,CNY,CNZ,
c     &           CT1X,CT1Y,CT1Z,CT2X,CT2Y,CT2Z,
c     &           DFRO,DFRUN,DFRUT1,DFRUT2,DFRET)
C
C
C********** Dual RN
C------------------------------------------------------
            FUNCEL = SURF * JTR(1,1)
            RR.AM(IFAC,2,2) = FUNCEL / VOLD
            RR.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(1,2)
            RUX.AM(IFAC,2,2) = FUNCEL / VOLD
            RUX.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(1,3)
            RUY.AM(IFAC,2,2) = FUNCEL / VOLD
            RUY.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(1,4)
            RUZ.AM(IFAC,2,2) = FUNCEL / VOLD
            RUZ.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(1,5)
            RRET.AM(IFAC,2,2) = FUNCEL / VOLD
            RRET.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
C********** Dual RUXN
C--------------------------------------------
            FUNCEL = SURF * JTR(2,1)
            UXR.AM(IFAC,2,2) = FUNCEL / VOLD
            UXR.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(2,2)
            UXUX.AM(IFAC,2,2) = FUNCEL / VOLD
            UXUX.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(2,3)
            UXUY.AM(IFAC,2,2) = FUNCEL / VOLD
            UXUY.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(2,4)
            UXUZ.AM(IFAC,2,2) = FUNCEL / VOLD
            UXUZ.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(2,5)
            UXRET.AM(IFAC,2,2) = FUNCEL / VOLD
            UXRET.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
C********** Dual RUYN
C--------------------------------------------
            FUNCEL = SURF * JTR(3,1)
            UYR.AM(IFAC,2,2) = FUNCEL / VOLD
            UYR.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(3,2)
            UYUX.AM(IFAC,2,2) = FUNCEL / VOLD
            UYUX.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(3,3)
            UYUY.AM(IFAC,2,2) = FUNCEL / VOLD
            UYUY.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(3,4)
            UYUZ.AM(IFAC,2,2) = FUNCEL / VOLD
            UYUZ.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(3,5)
            UYRET.AM(IFAC,2,2) = FUNCEL / VOLD
            UYRET.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
C********** Dual RUZN
C--------------------------------------------
            FUNCEL = SURF * JTR(4,1)
            UZR.AM(IFAC,2,2) = FUNCEL / VOLD
            UZR.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(4,2)
            UZUX.AM(IFAC,2,2) = FUNCEL / VOLD
            UZUX.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(4,3)
            UZUY.AM(IFAC,2,2) = FUNCEL / VOLD
            UZUY.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(4,4)
            UZUZ.AM(IFAC,2,2) = FUNCEL / VOLD
            UZUZ.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
            FUNCEL = SURF * JTR(4,5)
            UZRET.AM(IFAC,2,2) = FUNCEL / VOLD
            UZRET.AM(IFAC,2,1) = -FUNCEL / VOLG
C--------------------------------------------
C********** Dual RETN
C--------------------------------------------
            FUNCEL = SURF * JTR(5,1)
            RETR.AM(IFAC,2,2) = FUNCEL / VOLD
            RETR.AM(IFAC,2,1) = -FUNCEL / VOLG
C---------------------------------------------
            FUNCEL = SURF * JTR(5,2)
            RETUX.AM(IFAC,2,2) = FUNCEL / VOLD
            RETUX.AM(IFAC,2,1) = -FUNCEL / VOLG
C---------------------------------------------
            FUNCEL = SURF * JTR(5,3)
            RETUY.AM(IFAC,2,2) = FUNCEL / VOLD
            RETUY.AM(IFAC,2,1) = -FUNCEL / VOLG
C---------------------------------------------
            FUNCEL = SURF * JTR(5,4)
            RETUZ.AM(IFAC,2,2) = FUNCEL / VOLD
            RETUZ.AM(IFAC,2,1) = -FUNCEL / VOLG
C---------------------------------------------
            FUNCEL = SURF * JTR(5,5)
            RETRET.AM(IFAC,2,2) = FUNCEL / VOLD
            RETRET.AM(IFAC,2,1) = -FUNCEL / VOLG
C
         ELSE
C----------------------------------------------
C********** Murs (NGCG = NGCD)
C
C
C********** Les MELEMEs
C----------------------------------------------
            MELEDU.NUM(1,IFAC) = NGCG
            MELEDU.NUM(2,IFAC) = NGCD
            NLCG = MLENTC.LECT(NGCG)
C---------------------------------------------
            ROG = MPRN.VPOCHA(NLCG,1)
            PG = MPPN.VPOCHA(NLCG,1)
            UXG = MPUN.VPOCHA(NLCG,1)
            UYG = MPUN.VPOCHA(NLCG,2)
            UZG = MPUN.VPOCHA(NLCG,3)
            GAMG = MPGAMN.VPOCHA(NLCG,1)
            VOLG = MPVOLU.VPOCHA(NLCG,1)
C-------------------------------------------
            WVEC_L(1)=ROG
            WVEC_L(2)=UXG
            WVEC_L(3)=UYG
            WVEC_L(4)=UZG
            WVEC_L(5)=PG
C--------------------------------------------
C********** La normale sortante
C-------------------------------------------
c            SURF = MPOVSU.VPOCHA(NLCF,1)
c            CNX = MPNORM.VPOCHA(NLCF,7)
c            CNY = MPNORM.VPOCHA(NLCF,8)
c            CNZ = MPNORM.VPOCHA(NLCF,9)
cC
c            CALL VLHJ6(ROG,UXG,UYG,UZG,PG,GAMG,CNX,CNY,CNZ,
c     &           DFRUN)
C---------------------------------------------------------
            SURF = MPOVSU.VPOCHA(NLCF,1)
            NVECT(1) = MPNORM.VPOCHA(NLCF,7)
            NVECT(2) = MPNORM.VPOCHA(NLCF,8)
            NVECT(3) = MPNORM.VPOCHA(NLCF,9)
c--------------------------------------------
            TVECT1(1) = MPNORM.VPOCHA(NLCF,1)
            TVECT1(2) = MPNORM.VPOCHA(NLCF,2)
            TVECT1(3) = MPNORM.VPOCHA(NLCF,3)
c----------------------------------------------
            TVECT2(1) = MPNORM.VPOCHA(NLCF,4)
            TVECT2(2) = MPNORM.VPOCHA(NLCF,5)
            TVECT2(3) = MPNORM.VPOCHA(NLCF,6)
C------- COEFFICIENTS ----------------------------
            C11=TVECT1(2)*TVECT2(3)-TVECT1(3)*TVECT2(2)
            C12=NVECT(2)*TVECT2(3)-TVECT2(2)*NVECT(3)
            C13=NVECT(2)*TVECT1(3)-TVECT1(2)*NVECT(3)
C---------------------------------
            C21=TVECT1(1)*TVECT2(3)-TVECT1(3)*TVECT2(1)
            C22=NVECT(1)*TVECT2(3)-TVECT2(1)*NVECT(3)
            C23=NVECT(1)*TVECT1(3)-TVECT1(1)*NVECT(3)
C---------------------------------
            C31=TVECT1(1)*TVECT2(2)-TVECT1(2)*TVECT2(1)
            C32=NVECT(1)*TVECT2(2)-TVECT2(1)*NVECT(2)
            C33=NVECT(1)*TVECT1(2)-TVECT1(1)*NVECT(2)
            DET=NVECT(1)*C11-NVECT(2)*C21+NVECT(3)*C31
C---------------------------------
            ZC11=-NVECT(1)*C11-TVECT1(1)*C12+TVECT2(1)*C13
            ZC12=-NVECT(2)*C11-TVECT1(2)*C12+TVECT2(2)*C13
            ZC13=-NVECT(3)*C11-TVECT1(3)*C12+TVECT2(3)*C13
C---------------------------------
            ZC21=NVECT(1)*C21+TVECT1(1)*C22-TVECT2(1)*C23
            ZC22=NVECT(2)*C21+TVECT1(2)*C22-TVECT2(2)*C23
            ZC23=NVECT(3)*C21+TVECT1(3)*C22-TVECT2(3)*C23
C---------------------------------
            ZC31=-NVECT(1)*C31-TVECT1(1)*C32+TVECT2(1)*C33
            ZC32=-NVECT(2)*C31-TVECT1(2)*C32+TVECT2(2)*C33
            ZC33=-NVECT(3)*C31-TVECT1(3)*C32+TVECT2(3)*C33
C-------------------------------------------------
            ROD = ROG
            PD = PG
            UXD = (ZC11*UXG+ZC12*UYG+ZC13*UZG)/DET
            UYD = (ZC21*UXG+ZC22*UYG+ZC23*UZG)/DET
            UZD = (ZC31*UXG+ZC32*UYG+ZC33*UZG)/DET
            VOLD = VOLG
C------------------------------------------------
            WVEC_R(1)=ROD
            WVEC_R(2)=UXD
            WVEC_R(3)=UYD
            WVEC_R(4)=UZD
            WVEC_R(5)=PD
C-------------------------------------------
C********** La normale sortante
C-------------------------------------------
             V_INF=MAX(MPUINF.VPOCHA(NLCG,1),
     &              MPUPRI.VPOCHA(NLCG,1))
            CALL jw3Dbm(jtl,wvec_l,wvec_r,nvect,tvect1,tvect2,
     &                  gamg,v_inf)
C---------------------------------------------------------
C********** Dual RN
C---------------------------------------------------------
            RR.AM(IFAC,1,1) = 0.0D0
            RR.AM(IFAC,1,2) = 0.0D0
C---------------------------------------------------------
            RUX.AM(IFAC,1,1) = 0.0D0
            RUX.AM(IFAC,1,2) = 0.0D0
C---------------------------------------------------------
            RUY.AM(IFAC,1,1) = 0.0D0
            RUY.AM(IFAC,1,2) = 0.0D0
C---------------------------------------------------------
            RUZ.AM(IFAC,1,1) = 0.0D0
            RUZ.AM(IFAC,1,2) = 0.0D0
C--------------------------------------------------------
            RRET.AM(IFAC,1,1) = 0.0D0
            RRET.AM(IFAC,1,2) = 0.0D0
C---------------------------------------------------------
C********** Dual RUXN
C---------------------------------------------------------
            FUNCEL = SURF * JTL(2,1)
            UXR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXR.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(2,2)
            UXUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXUX.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(2,3)
            UXUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXUY.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(2,4)
            UXUZ.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXUZ.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(2,5)
            UXRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXRET.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
C********** Dual RUYN
C----------------------------------------------------
            FUNCEL = SURF * JTL(3,1)
            UYR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYR.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(3,2)
            UYUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYUX.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(3,3)
            UYUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYUY.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(3,4)
            UYUZ.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYUZ.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(3,5)
            UYRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYRET.AM(IFAC,1,2) = 0.0D0
C-----------------------------------------------------
C********** Dual RUZN
C-----------------------------------------------------
            FUNCEL = SURF * JTL(4,1)
            UZR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UZR.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(4,2)
            UZUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UZUX.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(4,3)
            UZUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UZUY.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(4,4)
            UZUZ.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UZUZ.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            FUNCEL = SURF * JTL(4,5)
            UZRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UZRET.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
C********** Dual RETN
C----------------------------------------------------
            RETR.AM(IFAC,1,1) = 0.0D0
            RETR.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            RETUX.AM(IFAC,1,1) = 0.0D0
            RETUX.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            RETUY.AM(IFAC,1,1) = 0.0D0
            RETUY.AM(IFAC,1,2) = 0.0D0
C----------------------------------------------------
            RETUZ.AM(IFAC,1,1) = 0.0D0
            RETUZ.AM(IFAC,1,2) = 0.0D0
C---------------------------------------------------
            RETRET.AM(IFAC,1,1) = 0.0D0
            RETRET.AM(IFAC,1,2) = 0.0D0
C---------------------------------------------------
C********** Dual RN
C---------------------------------------------------
            RR.AM(IFAC,2,2) = 0.0D0
            RR.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            RUX.AM(IFAC,2,2) = 0.0D0
            RUX.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            RUY.AM(IFAC,2,2) = 0.0D0
            RUY.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            RUZ.AM(IFAC,2,2) = 0.0D0
            RUZ.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            RRET.AM(IFAC,2,2) = 0.0D0
            RRET.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
C********** Dual RUXN
C---------------------------------------------------
            UXR.AM(IFAC,2,2) = 0.0D0
            UXR.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UXUX.AM(IFAC,2,2) = 0.0D0
            UXUX.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UXUY.AM(IFAC,2,2) = 0.0D0
            UXUY.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UXUZ.AM(IFAC,2,2) = 0.0D0
            UXUZ.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UXRET.AM(IFAC,2,2) = 0.0D0
            UXRET.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
C********** Dual RUYN
C---------------------------------------------------
            UYR.AM(IFAC,2,2) = 0.0D0
            UYR.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UYUX.AM(IFAC,2,2) = 0.0D0
            UYUX.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UYUY.AM(IFAC,2,2) = 0.0D0
            UYUY.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UYUZ.AM(IFAC,2,2) = 0.0D0
            UYUZ.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UYRET.AM(IFAC,2,2) = 0.0D0
            UYRET.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
C********** Dual RUZN
C---------------------------------------------------
            UZR.AM(IFAC,2,2) = 0.0D0
            UZR.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UZUX.AM(IFAC,2,2) = 0.0D0
            UZUX.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UZUY.AM(IFAC,2,2) = 0.0D0
            UZUY.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UZUZ.AM(IFAC,2,2) = 0.0D0
            UZUZ.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            UZRET.AM(IFAC,2,2) = 0.0D0
            UZRET.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
C********** Dual RETN
C---------------------------------------------------
            RETR.AM(IFAC,2,2) = 0.0D0
            RETR.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            RETUX.AM(IFAC,2,2) = 0.0D0
            RETUX.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            RETUY.AM(IFAC,2,2) = 0.0D0
            RETUY.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            RETUZ.AM(IFAC,2,2) = 0.0D0
            RETUZ.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
            RETRET.AM(IFAC,2,2) = 0.0D0
            RETRET.AM(IFAC,2,1) = 0.0D0
C---------------------------------------------------
         ENDIF
      ENDDO
C
      SEGDES MELEMC
      SEGDES MELEFE
      SEGDES MELEMF
C
      SEGDES MPOVSU
      SEGDES MPVOLU
      SEGDES MPNORM
C
      SEGDES MPRN
      SEGDES MPPN
      SEGDES MPUN
      SEGDES MPGAMN
C
      SEGDES MELEDU
      SEGDES MATRIK
      SEGDES IMATRI
C
      SEGDES   RR  , RUX , RUY , RUZ, RRET ,
     &         UXR , UXUX , UXUY , UXUZ, UXRET ,
     &         UYR , UYUX , UYUY , UYUZ, UYRET ,
     &         UZR , UZUX , UZUY , UZUZ, UYRET ,
     &         RETR , RETUX , RETUY , RETUZ, RETRET
 
      SEGSUP MLENTC
      SEGSUP MLENTF
      SEGDES MLMINC
      SEGSUP MLELIM
C
      SEGDES MPUPRI
      SEGDES MPUINF
      IF(MELLIM .NE.0) SEGDES MELLIM
C
 9999 CONTINUE
      RETURN
      END