Procédures | Cast3M

* INITEFMH  PROCEDUR  GOUNAND   05/02/16    21:15:51     5029           
**********************************************************************
'DEBP'  INITEFMH MoDARCY*'MMODEL' Porosite*'CHPOINT' MateDiff*'CHPOINT'
                 ChPSour*'CHPOINT' DeltaT*'FLOTTANT' Cini*'CHPOINT'
                 TetaDiff*'FLOTTANT' TetaConv*'FLOTTANT'
                 TetaLin*'FLOTTANT'  Qface/'CHPOINT' QELEM/'CHPOINT'
                 DISPL/'CHPOINT' DISPT/'CHPOINT'
                 LMLump*'LOGIQUE' DECENTR*'LOGIQUE' CHCLIM*'TABLE'
                 OPTRESOL/'TABLE';
* ATTENTION La vitesse est optionnelle, L'ordre est important
* et les types d'arguments qui se suivent aussi pour tester leur
* présence
*
* Attention il faudra transformer les vitesses en débits face
* et sortir le champ
*
*  |-----------------------------------------------------------------|
*  | Phrase d'appel (en GIBIANE)                                     |
*  |-----------------------------------------------------------------|
*  |                                                                 |
*  | SMTr MatrTr CoefDt TbDarTra MassEFMH nomespec                   |
*  | nbespece Difftot Tcini TABRES TABMODI= INITEFMH MoDARCY Porosite|
*  |                      MateDiff ChPSour DeltaT Cini TetaDiff      |
*  |                      TetaConv TetaLin fluimp dircli (QFACE)     |
*  |                      LMLump DECENTR CHCLIM optresol ;           |
*  |                                                                 |
*  |-----------------------------------------------------------------|
*  | Généralités : MATTEFMH construit la matrice de discrétisation   |
*  |               du problème de transport convection-diffusion pour|
*  |               le premier pas de tps d'un algorithme transitoire.|
*  |               Le second membre et les Conditions limites de flux|
*  |               sont pris en compte.                              |
*  |            RESTE TCINI, DECENTR et TERME LIN                    |
*  |-----------------------------------------------------------------|
*  |                                                                 |
*  |-----------------------------------------------------------------|
*  |                         ENTREES                                 |
*  |-----------------------------------------------------------------|
*  | MoDARCY  : modele Darcy.                                        |
*  |                                                                 |
*  | Porosite : champoint de composante 'CK'                         | 
*  |                                                                 |
*  | MateDiff : Tenseur de diffusion  (type iso, ..) champoint       |
*  |             de composante 'K' en isotrope, 'K11', 'K21',        |
*  |            'K22' en anisotrope 2d et  'K11', 'K21', 'K22', 'K31'|
*  |            'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQUE'|
*  |                                                                 |
*  | DeltaT   : Pas de temps.                                        |
*  |                                                                 |
*  | ChPSour  : Champ par points des sources volumiques par unité de |
*  |            temps (support maillage centre). Composante ??????   |
*  |                                                                 |
*  | Cini     : Concentration initiale, CHPOINT centre.              |
*  |            Composante 'H'.                                      |
*  |                                                                 |
*  | Qface    : vitesse aux faces, CHPO face de composantes Vx, Vy   |
*  |            en 2d et Vx, Vy, Vz en 3d. Il s'agit plus exatement  |
*  |            de (V.n)n, c'est à dire de la composante normale de  |
*  |            la vitesse aux faces. ???????? (je pressens que      |
*  |            castem va sortir des flux, cad intégrés sur surfaces)|
*  |                                                                 |
*  | TetaDiff : Valeur de theta pour theta-schéma en temps, opérateur|
*  |            de diffusion. Entre 0 et 1. 0 = explicite, 1 = euler |
*  |            implicite.                                           |
*  |                                                                 |
*  | TetaConv : Valeur de theta pour theta-schéma en temps, opérateur|
*  |            de convection. Entre 0 et 1. 0 = explicite, 1 = Euler|
*  |            implicite.                                           |
*  |                                                                 |
*  | TetaLin  : valeur de theta pour theta-schéma en temps, opérateur|
*  |            linéaire du type coef * C, où C est l'inconnue.      |
*  |            Entre 0 et 1. 0 = explicite, 1 = euler implicite.    |
*  |            ??????????? A voir car peut etre identique à Tetadiff|
*  |                                                                 | 
*  | LMLump   : Logique. Si vrai on effectue une condensation de     |
*  |            masse de la matrice EFMH                             |
*  |                                                                 |
*  | DECENTR  : Logique. Vrai veut dire schémas décentrés et faux    |
*  |            veut dire schéma convectif centré.                   |        
*  |                                                                 |          
*  | CHCLIM   : table d'indice 'NEUMANN' et 'DIRICHLET' contenant les|
*  |            Chpoint à n composantes contenant les conditions aux |
*  |            limites de Neumann et Dirichlet par espece.          |
*  |            L'indice 'FLUXTOT' contient les conditions limites   |
*  |            de flux total et 'FLUMIXTE' concerne une condition   |
*  |            de flux mixte : 'FLUMIXTE' . 'VAL' contient le champ |
*  |            à n composantes indiquant le flux, 'FLUMIXTE' . 'A'  |
*  |            et 'FLUMIXTE' . 'B' les coef (champoints SCAL) tels  |
*  |            que A D grad (C) + B (C) = VAL                       | 
*  |                                                                 |
*  | OPTRESOL : Table dont l'entree est optionnelle définissant      |     
*  |               les options de résolution pour 'KRES'.            |        
*  |                                                                 |
*  |-----------------------------------------------------------------|
*  |                         SORTIES                                 |
*  |-----------------------------------------------------------------|
*  |                                                                 |
*  |                                                                 |
*  | MassEFMH : matrice elementaire EFMH                             |
*  |                                                                 |
*  | MatrTr   : matrice globale sur les traces                       |
*  |                                                                 |
*  | SMTr     : second membre sur les traces                         |
*  |                                                                 |
*  | TbDarTra : table Darcy transitoire utilisée par MHYB, SMTP ...  |
*  |                                                                 |
*  | nomespec : liste des noms de composante des espèces dans Cini   |
*  |                                                                 |
*  | nbespece : nombre de composante de Cini, soit nombre d'especes  |
*  |                                                                 |
*  | nbsource : nombre de composantes du terme source qd X especes   |
*  |                                                                 |
*  | Diffdisp : Dipersivité, tenseur chpoint K11 K22 K33 K21 K31 K32 |
*  |                                                                 |
*  | TABRES   : Table complète définissant les options de résolution |    
*  |            pour 'KRES'.                                         |
*  |                                                                 |
*  | Tcini    : Trace de concentration aux faces (eventuellement à   |
*  |            plusieurs composantes (espèces)                      |
*  |                                                                 |
*  | TABMODI  : table contenant des logiques indiquant la nécessité  |
*  |            ou non de reclalculer certains termes.               |
*  |            'POROSITE' : VRAI si le coefficient devant D/DT      |
*  |                         (porosité) est modifié depuis le dernier|
*  |                         appel                                   |
*  |            'DELTAT'   : VRAI si le pas de tps a changé          |
*  |            'CONVECTI' : VRAI si la vitesse a changé             |
*  |            'COEF_LIN' : VRAI si le coef en facteur de C a changé|
*  |            'DIFFUSI'  : VRAI si les diffusivités ont changé     |
*  |                                                                 |
*  |-----------------------------------------------------------------|
*  |                  VARIABLES INTERNES                             |
*  |-----------------------------------------------------------------|
*  |                                                                 |
*  | CoefDt   : coeff devant dC/dt integre sur les elements          |
*  |                                                                 |
*  | PCONV    : Logique indiquant VRAI si présence de convection     |
*  |                                                                 |
*  | toltheta : 1.D-4  seuil en dessous duquel on considère que la   |
*  |            valeur de theta du theta-schéma est nulle (schéma    |
*  |            explicite) ou au contraire euler-implicite si        |
*  |            theta > 0.9999                                       |
*  |                                                                 |
*  | Tccini   : Trace de concentration aux faces (une composante)    |
*  |                                                                 |
*  | SSource  : Source aux centre (une composante)                   |
*  |                                                                 |
*  | CCini    : concentration aux centres (une composante)           |
*  |                                                                 |
*  | lstcps   : liste des noms de composante des espèces dans Chpsour|
*  |                                                                 |
*  | SSMTr     : second membre sur les traces pour une espèce        |
*  |                                                                 |
*  | MatConv   : matrice globale sur les traces pour la convection   |
*  |                                                                 |
*  | Numdiff  : diffusivité numérique due au décentrement            |
*  |                                                                 |
*  | FLUNEU   : LOGIQUE valant VRAI si conditions de Neumann         |
*  |                                                                 |
*  | CLFLUX   : Chpoint à n composantes contenant les flux imposés   |
*  |            pour chaque espece chimique. nul si pas de flux      |
*  |            OPTIONNEL                                            |
*  |                                                                 |
**********************************************************************
 
*---------------------------------------------------------------------
*---------- On récupere les conditions limites      ------------------
*---------------------------------------------------------------------
 
 
FLUNEU = FAUX                                                         ;
FLUTOT = FAUX                                                         ;
FLUMIX = FAUX                                                         ;
FLUCLI = FAUX                                                         ;
 
* Neumann
'SI' ('EXISTE' CHCLIM 'NEUMANN')                                      ;
   CLFLUX1 = CHCLIM . 'NEUMANN'                                       ;
   FLUNEU = VRAI                                                      ;
'SINON'                                                               ;
*  on crée un champ vide
   CLFLUX1 dum = 'KOPS' MATRIK                                        ;
   'OUBLIER' dum                                                      ;
'FINSI'                                                               ;
 
* Flux total
'SI' ('EXISTE' CHCLIM 'FLUTOTAL')                                     ;
   CLFLUX2 = CHCLIM . 'FLUTOTAL'                                      ;
   FLUTOT = VRAI                                                      ;
'SINON'                                                               ;
*  on crée un champ vide
   CLFLUX2 dum = 'KOPS' MATRIK                                        ;
   'OUBLIER' dum                                                      ;
'FINSI'                                                               ;
 
* Flux mixte
'SI' ('EXISTE' CHCLIM 'FLUMIXTE')                                     ;
*  comme on impose A Dgrad C + B C = flumix, on le traite sous
* la forme D grad C + (B/A) C = flumix/A plus naturelle en EFMH car
* D grad C est le flux diffusif
   COFA   = -1.D0 * CHCLIM . 'FLUMIXTE' . 'COEFA'                     ;
   CLFLUX3 = CHCLIM . 'FLUMIXTE' . 'VAL' '/' COFA                     ;
   CLFLUX3 = CHAN 'ATTRIBUT' CLFLUX3 NATURE DISCRET                   ;
   FLUMIX = VRAI                                                      ;
'SINON'                                                               ;
*  on crée un champ vide
   CLFLUX3 dum = 'KOPS' MATRIK                                        ;
   'OUBLIER' dum                                                      ;
'FINSI'                                                               ;
 
 
* On fabrique le terme de flux complet
'SI' (FLUNEU 'OU'  FLUTOT 'OU' FLUMIX)                                ;
   CLFLUX = CLFLUX1 'ET' CLFLUX2 'ET' CLFLUX3                         ;
   FLUCLI = VRAI                                                      ;
'FINSI'                                                               ; 
 
 
*---------------------------------------------------------------------
*---------- Initialisations de tables, coefficients ------------------
*---------------------------------------------------------------------
 
* Table de logiques indiquant des modifications. Initialisation
TABMODI = TABLE;
TABMODI . 'POROSITE' = FAUX; 
TABMODI . 'CONVECTI' = FAUX; 
TABMODI . 'DELTAT' = FAUX;
TABMODI . 'COEF_LIN' = FAUX;
TABMODI . 'DIFFUSIV' = FAUX;
 
 
'SI' ('EXISTE' QFACE) ;
   PCONV = VRAI;
   'SI' ('EXISTE' DISPL) ;
      DISPERSI = VRAI ;
    'SINON' ;
      DISPERSI = FAUX;
   'FINSI' ;
'SINON' ;
   PCONV = FAUX;
   DISPERSI = FAUX;
'FINSI' ; 
 
 
* tolerance sur theta du theta schéma de discrétisation en temps.
* il faudrait remmettre les theta à 0 ou 1 si nécessaire dans
* procédure amont.
toltheta = 1.D-4;
 
'SI' (TetaConv 'EGA' 0.D0 toltheta) ;
   TetaConv = 0.D0;
'FINSI' ;
'SI' (TetaConv 'EGA' 1.D0 toltheta) ;
   TetaConv = 1.D0;
'FINSI' ;
'SI' (TetaDiff 'EGA' 0.D0 toltheta) ;
   TetaDiff = 0.D0;
'FINSI' ;
'SI' (TetaDiff 'EGA' 1.D0 toltheta) ;
   TetaDiff = 1.D0;
'FINSI' ;
'SI' (TetaLin 'EGA' 0.D0 toltheta) ;
   TetaLin = 0.D0;
'FINSI' ;
'SI' (TetaLin 'EGA' 1.D0 toltheta) ;
   TetaLin = 1.D0;
'FINSI' ;
 
 
* Calcul du terme devant le dC/dt integré sur le volume
* CHAMELEM 'SCAL'. Voir pour VF ?????
CoefDt = ('NOMC' 'SCAL' Porosite);
COefDt = Porosite * ('DOMA' modarcy VOLUME);
CoefDt = 'KCHA' Modarcy CoefDt 'CHAM';
 
*  Creation de la table TbDarTra utilisée par les operateurs MATP, 
*  SMTP, HYBP et HDEBI, dans le cadre des EFMH. invariante si
*  le pas de temps ne bouge pas.
TbDarTra = 'TABLE' 'DARCY_TRANSITOIRE';
TbDarTra . 'THETA' = TetaDiff ;
* essayer si marche avec, sans indice !!!!
'SI' (PCONV) ;
   TbDarTra . 'THETA_CONVECTION'= TetaConv ;
'FINSI' ;
TbDarTra . 'PAS' = DeltaT ;
TbDarTra . 'SURF' = CoefDt ;
 
*---------------------------------------------------------------------
*----------------------- CREATION TABLE POUR RESOLUTION --------------
*---------------------------------------------------------------------
 
 
**************** OPTIONS PAR DEFAUT **************************
 
* création de la table de résolution pour l'opérateur KRES
* On crée la table de résolution avec les options par défaut
* On y remplacera les valeurs définit par l'utilisateur ensuite.
TABRES = 'TABLE' 'METHINV' ;
* type d'inversion
'SI' (PCONV) ;
*  option BCGSTAB par défaut pour une matrice non symétrique
   METHRES = 3;
'SINON' ;
*  option gradient conjugué par défaut pour une matrice symétrique   
   METHRES = 2;
'FINSI' ;
TABRES . 'TYPINV' = METHRES ;
 
* niveau d'impression.
TABRES . 'IMPINV' = 0 ;
 
* Type de renumérotation. Option SLOANE par défaut
*TABRES . 'TYRENU' = 'SLOANE' ;
TABRES . 'TYRENU' = 'SLOA' ;
 
* La gestion des multiplicateurs sera modifiée
* par la suite. Pas d'option pour l'instant
TABRES . 'PCMLAG' = 'APR2' ;
TABRES . 'OUBMAT' = 0 ;
TABRES . 'SCALING' = 0 ;
 
*INDICES SPÉCIFIQUES POUR UNE MÉTHODE ITÉRATIVE
* Nombre maxi d'itérations
TABRES . 'NITMAX' = 1500 ;
* résidu pour la convergence de la méthode
TABRES . 'RESID' = 1.D-15 ;
* valeur minimale du pivot de la méthode
TABRES . 'BCGSBTOL' = 1.D-120 ;
* preconditionnement ILU(0)
TABRES . 'PRECOND' = 3 ;
*relaxation pour MILU0
TABRES . 'MILURELX' = 1.D0 ;
*GMRESTART
TABRES . 'GMRESTRT' = 100 ;
*ILUTLFIL
TABRES . 'ILUTLFIL' = 2;
* drop tolerence pour ILUT2
TABRES . 'ILUTDTOL' = 0.D0;
TABRES . 'ILUTPTOL' = 0.01D0;
TABRES . 'ILUTALPH' = 0.D0;
 
************** OPTIONS UTILISATEUR **************************
 
 
'SI' ('EGA' ('TYPE' OPTRESOL) 'TABLE') ;
* L'utilisateur a défini des options pour la méthode
* de résolution.
 
*  Type d'inversion
   'SI' ('EXISTE' OPTRESOL 'TYPINV') ;
      TABRES . 'TYPINV' = OPTRESOL . 'TYPINV';
   'FINSI' ;
 
*  Niveau d'impression
   'SI' ('EXISTE' OPTRESOL 'IMPINV') ;
      TABRES . 'IMPINV' = OPTRESOL . 'IMPINV';
   'FINSI' ;
 
* Type de renumérotation.
   'SI' ('EXISTE' OPTRESOL 'TYRENU') ;
      TABRES . 'TYRENU' = OPTRESOL . 'TYRENU'; 
   'FINSI' ;
 
* Indices spécifiques aux méthodes itératives
   'SI' ((TABRES . 'TYPINV') > 1); 
*     Nombre maxi d'iterations
      'SI' ('EXISTE' OPTRESOL 'NITMAX') ;
         TABRES . 'NITMAX' = OPTRESOL . 'NITMAX'; 
      'FINSI' ;
*     Valeur du résidu de la méthode
      'SI' ('EXISTE' OPTRESOL 'RESID') ;
         TABRES . 'RESID' = OPTRESOL . 'RESID'; 
      'FINSI' ;
*     valeur minimal du pivot de la méthode
      'SI' ('EXISTE' OPTRESOL 'BCGSBTOL') ;
         TABRES . 'BCGSBTOL' = OPTRESOL . 'BCGSBTOL'; 
      'FINSI' ;
*     precond par diagonale
      'SI' ('EXISTE' OPTRESOL 'PRECOND') ;
         TABRES . 'PRECOND' = OPTRESOL . 'PRECOND';
      'FINSI' ;
*     precon ILUT2
      'SI' ('EXISTE' OPTRESOL 'ILUTLFIL')  ;
         TABRES . 'ILUTLFIL' = OPTRESOL . 'ILUTLFIL' ;
      'FINSI' ;
   'FINSI' ;
*  Pour GMRES
   'SI' ((TABRES . 'TYPINV') EGA 5);
      'SI' ('EXISTE' OPTRESOL 'GMRESTRT') ;
         TABRES . 'GMRESTRT' = OPTRESOL . 'GMRESTRT';
      'SINON' ;
         TABRES . 'GMRESTRT' = 50;
      'FINSI' ;      
   'FINSI' ;
'FINSI' ;
 
 
SI (('EGA' TABRES . 'PRECOND' 8) 'OU' ('EGA' TABRES . 'PRECOND' 7));
   TABRES . 'ILUTDTOL' = 0.1D-2;
FINSI ;
 
*---------------------------------------------------------------------
*--------------------- CALCUL DE LA DISPERSIVITE----------------------
*---------------------------------------------------------------------
 
* On test si le modèle est anisotrope ou non (plus d'une composante si
* anisotrope)
compmat = 'EXTRAIRE' Modarcy 'MATERIAU' ;
ANISO   = ('DIME' compmat) > 1 ;
*'SI' ('NON' ANISO)                                                     ;
*   'SI' (DISPERSI 'OU' DECENTR)                                        ;
*      'MESSAGE' 'ERREUR - le modèle est déclaré isotrope en présence
*      de décentrement ou de dispersivité'                              ;
*      ERREUR 5                                                         ;
*   'SINON'                                                             ;
*      difftot = 'NOMC' 'K' Matediff                                    ;      
*   'FINSI'                                                             ;
*'FINSI'                                                                ;
 
* Remise de la diffusion aux bonnes composantes aniso car + général
 
* la ligne suivante ne marche pas, je dois introduire zozo ????
*Matediff = DIFFANIS Matediff 'EFMH'                                   ;
zozo = DIFFANIS Matediff 'EFMH' ANISO                                  ;
Matediff = zozo;
 
* Seulement si présence de convection et présence de
* dispersivité et disp_l et disp_t
 
'SI' (DISPERSI 'ET' PCONV)                                             ;
   DISPL = 'NOMC' 'SCAL' DISPL                                         ;
   DISPT = 'NOMC' 'SCAL' DISPT                                         ;
*   on calcul la dispersivité
   diffdisp = CALCDISP QELEM DISPL DISPT                               ;
'SINON'                                                                ;
   diffdisp = 0.D0 * Matediff                                          ;
'FINSI'                                                                ;
 
Difftot = Matediff + diffdisp                                          ;
 
*---------------------------------------------------------------------
*--------------------- CALCUL DU DECENTREMENT ------------------------
*---------------------------------------------------------------------
 
* Seulement si présence de convection et option décentrée
* Deuxième matrice masse ????
 
 
'SI' (DECENTR 'ET' PCONV)                                              ;
 
  'MESSAGE' 'on utilise un schéma décentré pour la convection'         ;
 
*  projection sur axe des x de la normale
   normax = 'DOMA' modarcy 'NORMALE'                                   ;
   normax = ('EXCO' 'UX' normax) 'ABS'                                 ;
   normax = 'NOMC' 'FLUX' normax                                       ;
   normax = ('DOMA' modarcy SURFACE) * normax                          ;
   dxmail = 'DIVU' modarcy normax ('ABS' ('DOMA' modarcy ORIENTATION)) ;
*  x 2 car surface volume / (surf haut + surf bas)
   dxmail = 2.D0 * ('DOMA' modarcy 'VOLUME') '/' (dxmail)              ;  
 
*  projection sur axe des x de la normale
   normay = 'DOMA' modarcy 'NORMALE'                                   ;
   normay = ('EXCO' 'UY' normay) 'ABS'                                 ;
   normay = 'NOMC' 'FLUX' normay                                       ;
   normay = ('DOMA' modarcy SURFACE) * normay                          ;
   dymail = 'DIVU' modarcy normay ('ABS' ('DOMA' modarcy ORIENTATION)) ;
   dymail = 2.D0 * ('DOMA' modarcy 'VOLUME') '/' (dymail)              ;   
 
*  porjection sur axe des x de la normale
   'SI' ('EGA' ('VALEUR' 'DIME') 3)                                    ;
      normaz = 'DOMA' modarcy 'NORMALE'                                ;
      normaz = ('EXCO' 'UZ' normaz) 'ABS'                              ;
      normaz = 'NOMC' 'FLUX' normaz                                    ;
      normaz = ('DOMA' modarcy SURFACE) * normaz                       ;
      dzmail = 'DIVU' modarcy normaz ('ABS' ('DOMA' modarcy
                ORIENTATION))                                          ;
      dzmail = 2.D0 * ('DOMA' modarcy 'VOLUME') '/' (dzmail)           ;  
   'FINSI'                                                             ;
 
*  Vitesses extraites par composantes et comme chpo centre
   V1 = 'EXCO' Qelem 'VX' 'SCAL'                                       ;
   V2 = EXCO Qelem 'VY' 'SCAL'                                         ;
   'SI' ('EGA' ('VALEUR' 'DIME') 3)                                    ;
      V3 = EXCO Qelem 'VZ' 'SCAL'                                      ;
   'FINSI'                                                             ;
 
Pe = OPTRESOL .  'PECLET'                                              ; 
*---on incorpore la diffusion numerique suivant X si PECLET > 2 ;
   dum = 'ABS' (dxmail * V1 '/' Pe)                                    ;
   D11P = exco 'K11' difftot 'SCAL'                                    ;
   m1 = masq D11P SUPERIEUR dum                                        ;
   D11P = m1*D11P + ((1 - m1) * dum)                                   ;
   D11P = 'NOMC' 'K11' D11P                                            ;
 
*---on incorpore la diffusion numerique suivant Y si PECLET > 2 ;
   dum = 'ABS' (dymail * V2 '/' Pe)                                    ;
   D22P = exco 'K22' difftot 'SCAL'                                    ;
   m1 = masq D22P SUPERIEUR dum                                        ;
   D22P = m1*D22P + ((1 - m1) * dum)                                   ;
   D22P = 'NOMC' 'K22' D22P                                            ;
 
*---on incorpore la diffusion numerique suivant Z si PECLET > 2 ;
   'SI' ('EGA' ('VALEUR' 'DIME') 3)                                    ;
      dum = 'ABS' (dzmail * V3 '/' Pe)                                 ;
      D33P = exco 'K33' difftot 'SCAL'                                 ;
      m1 = masq D33P SUPERIEUR dum                                     ;
      D33P = m1*D33P + ((1 - m1) * dum)                                ;
      D33P = 'NOMC' 'K33' D33P                                         ;
 
      diffdisp = D11P '+' D22P '+' D33P
               '+' ('EXCO' 'K21' difftot 'K21')
               '+' ('EXCO' 'K31' difftot 'K31')
               '+' ('EXCO' 'K32' difftot 'K32')  
                - Matediff                                             ;
   'SINON'                                                             ;
      diffdisp = D11P '+' D22P '+' ('EXCO' 'K21' difftot 'K21') 
               - Matediff                                              ;
   'FINSI'                                                             ;
 
'FINSI'                                                                ;
 
* diffusion totale (dispers + effective + numérique)
 
difftot = Matediff '+' diffdisp                                        ;
 
* changement en chamelem pour EFMH
 
difftot = 'KCHA' Modarcy difftot 'CHAM'                                ;
difftot = CHANGER difftot TYPE CARACTERISTIQUES                        ;
 
SI (ANISO)                                                             ;
'SI' ('EGA' ('VALEUR' 'DIME') 2)                                       ;
   difftot = 'MATE' Modarcy  'DIRECTION'
      (1. 0. ) 'PARALLELE' 'K11' ('EXCO' 'K11' difftot)
                        'K21' ('EXCO' 'K21' difftot)
                        'K22' ('EXCO' 'K22' difftot)                   ;
'FINSI'                                                                ;
'SI' ('EGA' ('VALEUR' 'DIME') 3)                                       ;
   difftot = 'MATE' Modarcy  'DIRECTION'
     (1. 0. 0.) (0. 1. 0.) 'PARALLELE'
                        'K11' ('EXCO' 'K11' difftot)
                        'K21' ('EXCO' 'K21' difftot)
                        'K22' ('EXCO' 'K22' difftot)
                        'K31' ('EXCO' 'K31' difftot)
                        'K32' ('EXCO' 'K32' difftot)
                        'K33' ('EXCO' 'K33' difftot)                   ;        
'FINSI'                                                                ;   
SINON                                                                  ;
 difftot = 'MATE' Modarcy 'K' difftot                                  ;
FINSI                                                                  ;
 
 
*---------------------------------------------------------------------
*-------------- Matrice masse inverse des EFMH -----------------------
*---------------------------------------------------------------------
 
 
* Calcul des  matrices de masse elementaires inverses
'SI' (LMLump) ;
*   masse lumping
    MassEFMH = 'MHYB' MoDARCY  Difftot 'LUMP' ;
'SINON' ;
    MassEFMH = 'MHYB' MoDARCY  Difftot ;
'FINSI' ;
 
 
 
*----------------------------------------------------------------------
*---------- CALCUL DE LA MATRICE DE DIFFUSION DU PROBL ----------------
*----------------------------------------------------------------------
 
 
* Calcul de la matrice du probleme diffusion transitoire
MatrTr =  'MATP' MoDARCY MassEFMH  TbDarTra ;
 
 
 
*----------------------------------------------------------------------
*--- CALCUL DE LA MATRICE DE CONVECTION ET DES SECONDS MEMBRES --------
*----------------------------------------------------------------------
 
 
************** Nombres d'espèces à gérer ******************************
 
* évaluation du nombre d'espèces.
nomespec = 'EXTRAIRE' Cini 'COMP';
nbespece   =  'DIME' nomespec;
* nb de termes sources, commun (1), ou egal au nombre
* d'especes (nbcompi)
lstcps = 'EXTRAIRE' Chpsour 'COMP';
nbsource   =  'DIME' lstcps;
'SI' ((nbsource 'NEG' 1) 'ET' (nbsource 'NEG' nbespece)) ;
   'MESSAGE' 'La source doit avoir le meme nombre de composantes';
   'MESSAGE' 'que les espèces ou 1 seule composante';
   ERREUR 5;
'FINSI' ;
 
 
* La matrice de convection ne dépend pas de Tcini mais
* est calculée en meme temps que le calcul du second
* membre. On effectue donc un traitement particulier
* dans un cas multiespece pour gagner en temps de calcul.
 
 
********************* Une ou plusieurs espèces ************************
 
'SI' (nbespece > 0) ;
 
*  pour un schéma en temps non euler implicite, il faut
*  la trace à l'instant précédent pour le second membre
*  pour la convection ou la diffusion
   'SI' ((TetaDiff 'NEG' 1.D0 toltheta)
         'OU' (TetaConv 'NEG' 1.D0 toltheta)) ;
*     Calcul de Tcini
       tcini dummy = CALCTRAC MoDARCY Difftot Cini                 
                     nomespec nbespece  LMLump
                     TABRES Tbdartra CHCLIM;
       'OUBLIER'  dummy;
   'SINON' ;
*     La trace de charge n'est pas réellement utilisée car multipliée
*     par 0. On la met à 0 pour simplifier les calculs.
      TCCINI = 'NOMC' 'TH' ('KCHT' Modarcy SCAL 'FACE' 0.D0);
      'REPETER' bloctc nbespece;
         'SI' (&bloctc 'EGA' 1);
            tcini =  ('NOMC' ('EXTRAIRE' &bloctc nomespec)
                   TCCINI);
         'SINON';
            tcini = tcini 'ET' ('NOMC' ('EXTRAIRE'
                    &bloctc nomespec)  ('COPIER' TCCINI));
         'FINSI' ;
      'FIN' bloctc;
   'FINSI';
 
   'REPETER' bloc1 nbespece;
*     On extrait la composante de Cini, Tcini et de la source
      CCini = 'NOMC' 'H' ('EXCO' (extr &bloc1 nomespec) Cini);
      TCCini = 'NOMC' 'TH' ('EXCO' (extr &bloc1 nomespec) Tcini);
      'SI' (nbsource > 1) ;
         SSource = 'NOMC' 'SOUR' ('EXCO' (extr &bloc1 nomespec)
                               Chpsour);
      'SINON' ;
         SSource = 'NOMC' 'SOUR' Chpsour;
      'FINSI' ;    
*     Conditions initiales
      TbDarTra . 'CHARGE' = CCini ;
      TbDarTra . 'TRACE'= TCCini ;
*     Prise en compte du terme source et eventuellement 
*     de la convection avec le schema centre 
      'SI' (PCONV);
*        convection
         'SI' (TetaConv 'NEG' 0.0D0 toltheta);
*           schéma partiellement implicite, matrice MatConv
            MatConv SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource
                                  (NOMC 'FLUX' QFace);
         'SINON' ;
*           schéma explicite, calcul du second membre uniquement
            SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource 
                                      (NOMC 'FLUX' QFace) ;
         'FINSI' ;
      'SINON';
*        pas de convection, calcul du second membre restant
         SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource;   
      'FINSI';
*     On reconstitue un champ de second membre
      'SI' ((&bloc1) 'EGA' 1) ;
         SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) SSMTr;
         SMTR =  SSMTr;
      'SINON' ;
         SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) ('COPIER' SSMTr);
         SMTr = SSMTr 'ET'  SMTr;
      'FINSI' ;
   'FIN' bloc1;
*  on stoque la matrice en assemblant la matrice de convection
*  du dernier calcul (elles sont toutes identiques). On mettra
*  un jour une option pour les sortir si besoin seulement.
   'SI' (PCONV 'ET' (TetaConv 'NEG' 0.0D0 toltheta));
      MatrTr  = MatrTr   'ET' MatConv ;
      'DETRUIT' MatConv;
      'MENAGE' ;
   'FINSI' ;
'FINSI' ;
 
 
 
*---------------------------------------------------------------------
*------ Conditions aux limites de Flux (mixtes et Neumann) -----------
*---------------------------------------------------------------------
 
 
'SI' (FLUCLI) ;  
   SMTr = SMTR 'ET' CLFLUX ;
'FINSI' ;
 
'MENAGE' ;
 
'FINP' SMTr MatrTr TbDarTra MassEFMH nomespec
      nbespece nbsource Diffdisp Tcini TABRES TABMODI;