* fichier :  darcy5.dgibi
************************************************************************
* Section : Fluides Darcy
************************************************************************
*
********************  CAS TEST : darcy5.dgibi  **********************
*
 GRAPH = FAUX   ;
*GRAPH = VRAI   ;
'SAUT' 'PAGE' ;
*
*-------------------------------------------------------------------
* TEST DARCY5
*
*    Validation de la résolution des équations de DARCY formulées en
*             (VITESSE - PRESSION) par comparaison  avec 
* la résolution des équations de DARCY formulées en (VITESSE - CHARGE)
*
*
*    (Résolution par une méthode d'éléments finis mixtes hybrides)
*
* Le maillage est composés de triangles et de quadrangles et contient
* des lignes de maillages courbes 
*
* la massif est supposé isotrope, indéformable et saturé
*
* Pour permettre la comparaison, la densité de l'eau reste constante 
*
*-------------------------------------------------------------------
*
 OPTI ECHO 0                                              ;
*
'SAUT' 'PAGE' ;
*
*- Options générales de calcul
*
TITRE 
'EFMY DARCY FORMULATION (VITESSE - PRESSION) 2D : darcy5.dgibi' ;

OPTI DIME 2 ELEM QUA4                                    ;
OPTI ISOV SURF                                           ;

*
*- Données physiques
*
*- Densité de l'eau
RHO0 = 1000. ;
*
*- Pression athmosphérique
P0 = 1.013E5 ;
*
*- Composante du vecteur gravité
Gx0 = 0. ;
Gy0 = -9.81 ;
*
*- Perméabilité intrinsèque
LA0 = 1.e-9 ;
*
*======================================================================
*
*- Construction du maillage
*
*======================================================================
*
*- Discrétisation spatiale
*
ENX = 15          ;
ENY0 = 5         ;
ENY1 = 5        ;
ENY2 = 5         ;
*
*- Création des points
*
A0 = 0.d0 -2.d3 ; B0 = 5.d3 -2.d3 ;
*
A1 = 0.d0 -1.d3 ; B1 = 5.d3 -1.d3 ; C1 = 2.d3 -1.25d3 ;
*
A2 = 0.d0 0.0d3 ; B2 = 5.d3 0.0d3 ; C2 = 2.d3 0.25d3 ;
*
A3 = 0.d0 1.5d3 ; B3 = 5.d3 1.d3 ; C3 = 2.5d3 1.5d3 ;
*
*- Création des droites
*
AB0 = DROI ENX A0 B0 ;
AB1 = A1 cer3 C1 ENX B1 ;
AB2 = A2 cer3 C2 ENX B2 ;
AB3 = A3 cer3 C3 ENX B3 ;
*
*- Creation des SURFACES
*
NIVEAU0 = AB0 regle ENY0  AB1 ;
OPTI ELEM TRI3 ;
NIVEAU1 = AB1 regle ENY1  AB2 ;
OPTI ELEM QUA4 ;
NIVEAU2 = AB2 regle ENY2  AB3 ;  
MASSIF0 = NIVEAU0 ET NIVEAU1 ET NIVEAU2  ;
HAUT  = INVE AB3 ;
MASSIF0 = ORIENT MASSIF0        ;
QMASSIF = CHANGE MASSIF0 QUAF ;
QFHAUT  = CHANGE HAUT QUAF ;
 ELIM 0.01 (QMASSIF ET QFHAUT ) ;
*
*- Domaines
*
*HYTOT = DOMA MASSIF0 0.01   ;
*CHYB1 = DOMA HYTOT 'SURFACE' ;
*CHYB2 = DOMA HYTOT 'NORMALE' ;
*MCHYB = DOMA HYTOT 'ORIENTAT' ;
*HYHAUT = DOMA HAUT INCL HYTOT 0.01 ;
*
*- Modélisation
*
*MODHYB = MODE HYTOT DARCY ISOTROPE HYQ4  HYT3                 ;
 MODHYB = MODE QMASSIF DARCY ISOTROPE                 ;
 MODHAU = MODE QFHAUT DARCY ISOTROPE                 ;
 CHYB1 = DOMA MODHYB 'SURFACE' ;
 CHYB2 = DOMA  MODHYB 'NORMALE' ;
 CEHAUT= DOMA MODHAU 'CENTRE' ;
 CETOT= DOMA MODHYB 'CENTRE' ;
*
*- perméabilité intrinsèque
*
LMME =  MANU CHPO CETOT 1 'K' la0 'NATURE' 'DIFFU' ;
*
*======================================================================
*
*- Résolution en  en (VITESSE ; PRESSION)
*
*======================================================================
*
LMMA = KCHA MODHYB 'CHAM' LMME                       ;
MATI_P  = MATERIAU  MODHYB 'K' LMMA                 ;
*
****  Masse HYBride
CND1A_P  = MHYBR MODHYB  MATI_P                ;
*
****  Masse FORce
M_P  = MHYBR MODHYB 'MASSE'                ;
*
****  MAtrice en Trace de pression TP
HND1A_P  = MATP MODHYB  CND1A_P                ;
*
*- Conditions aux limites
*
BHAUT_P = BLOQUE CEHAUT 'TH' ;
*
*- TP imposee
*
PIMP = MANU CHPO CEHAUT 1 'TH' P0 'NATURE' 'DISCRET' ;
EHAUT_P = DEPI BHAUT_P PIMP  ; 
*
*-  contribution des forces de volume au second membre 
*
RHO = MANU CHPO CETOT 1 'SCAL' RHO0 'NATURE' 'DISCRET' ;
RCH =  MANU CHPO CETOT 2 'FX' gx0 'FY' gy0  'NATURE' 'DISCRET' ;
RCH =  RHO * RCH ;
GRAV2TP = SQTP MODHYB  HND1A_P M_P RCH ;
*
*- Assemblage matrice et second membre
*
CCC1_P = HND1A_P ET BHAUT_P ;
* second membre du système matriciel en trace de charge
FFF1_P = GRAV2TP et EHAUT_P ;
*
*- Resolution en trace de pression
*
CHTER1_P = RESOUDRE CCC1_P FFF1_P ;
CHTER1_P = exco CHTER1_P 'TH' 'TH';
*
*- Calcul des pressions élémentaires
*
PCEN1 = HYBP MODHYB  CND1A_P CHTER1_P M_P RCH  ;
*
*- Calcul des débits
*
QFACE1_P = HDEBI MODHYB  CND1A_P PCEN1 CHTER1_P  M_P RCH;
VCENT1_P = HVIT  MODHYB  QFACE1_P  ;
*VV = VECT VCENT1_P 0.5e9 'VX' 'VY'  VERT ;
*titre 'VITESSES AUX CENTRES DES ELEMENTS : ECHL=1e9:2' ;
*TRAC VV CETOT ;
*
*=======================================================================
*
*- Résolution en (VITESSE ; CHARGE)
*
*=======================================================================
*
*
*- conductivité hydraulique
*
LMME = LMME * (ABS Gy0) * RHO0  ;
LMMA = KCHA MODHYB 'CHAM' LMME                       ;
MATI_H  = MATERIAU  MODHYB 'K' LMMA                 ;
*
****  Masse HYBride
CND1A_H  = MHYBR MODHYB  MATI_H                ;
*
****  MAtrice en Trace de charge TH
HND1A_H  = MATP MODHYB  CND1A_H                ;
*
*- Conditions aux limites
*
BHAUT_H = BLOQUE CEHAUT 'TH' ;
*
*- TH imposee
*
PIMP = NOMC 'TH' (CEHAUT COOR 2) ;
EHAUT_H = DEPI BHAUT_H PIMP  ;   
*
*- Assemblage matrice et second membre
*
CCC1_H = HND1A_H ET BHAUT_H ;
* second membre du système matriciel en trace de charge
FFF1_H = EHAUT_H ;
*
*- Resolution en trace de charge
*
CHTER1_H = RESOUDRE CCC1_H FFF1_H ;
CHTER1_H = exco CHTER1_H 'TH' 'TH';
*
*- Calcul de la charge
*
TCEN1 = HYBP MODHYB  CND1A_H CHTER1_H   ;
*
*- Calcul de la vitesse
*
QFACE1_H = HDEBI MODHYB  CND1A_H TCEN1 CHTER1_H  ;
VCENT1_H = HVIT  MODHYB  QFACE1_H  ;
*VV = VECT VCENT1_H 0.5e9 'VX' 'VY'  VERT ;
*titre 'CALCUL EN CHARGE : VITESSES AUX CENTRES DES ELTS : ECHL=1e9:2';
*TRAC VV CETOT ;
*
*=======================================================================
*
*- Comparaison des deux résolutions
*
*=======================================================================
*
*
*- comparaison des traces de charges
*
TP2 =  (NOMC CHTER1_P SCAL) - P0 / (rho0 * gy0)  ;
TP2 = ((DOMA MODHYB 'FACE') COOR 2) - TP2 ;
ER2 = (NOMC CHTER1_H SCAL)/TP2 - 1. ;
maxtp1 = ER2 ABS MAXI ;
*
*
*- comparaison des charges
*
PCEN2 = (NOMC PCEN1 SCAL) - P0 / (rho0 * gy0)  ;
PCEN2 = CETOT COOR 2 - PCEN2 ;
ER2 = (NOMC TCEN1 SCAL)/PCEN2 - 1. ;
maxp1 = ER2 ABS MAXI ;

'SI' GRAPH  ;
TITRE 'VARIATIONS RELATIVES DES CHARGES DH/H' ; 
TRAC (KCHA MODHYB 'CHAM' er2) MODHYB (CONTOUR
 (DOMA MODHYB 'MAILLAGE')) ;
'FINSI' ;
*
*- comparaison des vitesse
*
QFACE1_P = NOMC SCAL QFACE1_P                ;
VFACE1_P = QFACE1_P * CHYB2 / CHYB1          ;
QFACE1_H = NOMC SCAL QFACE1_H                ;
VFACE1_H = QFACE1_H * CHYB2 / CHYB1          ;
maxvf1 = (kops QFACE1_P '-' QFACE1_H) / (maxi QFACE1_H) ABS MAXI ;

MOT1 = 'MOTS' 'VX' 'VY'             ;
VDVD = 'PSCA' VCENT1_H VCENT1_H MOT1 MOT1 ;
VD1  = VCENT1_P - VCENT1_H           ;
VC1  = 'PSCA' VD1 VD1 MOT1 MOT1 ;
SDC1 = 'ABS' ( VC1 / VDVD )     ;
SDC1 = SDC1 '**' 0.5D0          ;
maxvc1 = MAXI SDC1 ;

'SI' GRAPH  ;
TITRE ' VARIATIONS RELATIVES DES VITESSES |DV|/|V|' ; 
TRAC ('KCHA' MODHYB 'CHAM' SDC1) MODHYB (CONTOUR
 (DOMA MODHYB 'MAILLAGE')) ;
'FINSI' ;

*
'SAUT' 'PAGE'                                                        ;
'SAUT' 2 'LIGNE'                                                     ;
'MESS' '                       VARIATIONS RELATIVES                '   ;
'SAUT' 1 'LIGNE'                                                     ;
'MESS' '      TH                H                Vf              Vc' ;
'SAUT' 1 'LIGNE'                                                     ;
'MESS' ' ' maxtp1  ' ' maxp1 ' ' maxvf1 ' ' maxvc1                   ;
'SAUT' 2 'LIGNE'                                                     ;
*
EPS0 = 1.E-7 ; 
LOG1 = MAXTP1 > EPS0                      ;
LOG2 = MAXP1  > EPS0                      ;
LOG3 = MAXVF1 > EPS0                      ;
LOG4 = MAXVC1 > EPS0                      ;
L0   = LOG1 'OU' LOG2 'OU' LOG3 'OU' LOG4 ;
*
'SI' ( L0 ) ;
   'ERRE' 5 ;
'SINO'      ;
   'ERRE' 0 ;
'FINS'      ;

FIN ;