Accueil | Cast3M

* fichier :  darcy1.dgibi
************************************************************************
************************************************************************
*****************  CAS TEST : darcy1.dgibi  *************************
*
GRAPH = 'N'   ;
'SAUT' 'PAGE' ;
*
*-------------------------------------------------------------------
* TEST DARCY1
*
*                  CALCUL DARCY ISOTROPE
*    Résolution par une méthode d'éléments finis mixtes hybrides.
*
* Ce test permet de vérifier le bon fonctionnement des opérateurs
* utilisés pour la résolution des équations de DARCY par une méthode
* d'éléments finis mixtes hybrides dans CASTEM2000.
*
* On effectue trois calculs sur un domaine carré, maillé par des
* quadrangles réguliers. Le phénomène étudié est monodimensionnel.
* Les conditions aux limites varient suivant le cas considéré :
* on impose le flux ou la charge sur les cotés du domaine.
*
* La solution analytique en charge étant un polynome de degré un
* en x, la vitesse est constante et horizontale.
* Si on note L la longueur du domaine, HD (resp. HG) la charge à
* droite du domaine (resp. à gauche) et dH=HD-HG on a :
*        h(x) = (x-L) * dH/L + HD
* Et par suite
*        V(x) = (-K*dH/L ; 0.)
*
* On s'attend à une précision de l'ordre de la précision machine.
*
*-------------------------------------------------------------------
*
'SAUT' 'PAGE' ;
*
*- Options générales de calcul.
*
'TITR' 'EFMH DARCY ISOTROPE 1D : darcy1.dgibi' ;
* Dimension 2, éléments à générer quadrangulaires ou triangulaires :
'OPTI' 'DIME' 2 'ELEM' 'QUA4'                  ;
* Pas d'affichage à l'écran des lignes de commandes lues :
'OPTI' 'ECHO' 1                                ;
*                                                     ------------
*                                                     = MAILLAGE =
*                                                     ------------
*
*- Création des points supports des DROITES
*
* Dimensions :
L  = 1.D0    ;
HA = 1.D0    ;
 
* Coordonnées :
XG = 0.D0    ;
XD = XG + L  ;
YG = 0.D0    ;
YD = YG + HA ;
 
* Points
A1 = XG YG   ;
A3 = XD YG   ;
D1 = XG YD   ;
D3 = XD YD   ;
*
*- Création des DROITES frontières
*
* Discrétisation dans les deux directions :
INUMX = 9                  ;
INUMY = 4                  ;
* droites maillées :
DRBAS = A3 'DROI' INUMX A1 ;
DRGAU = A1 'DROI' INUMY D1 ;
DRHAU = D1 'DROI' INUMX D3 ;
DRDRO = D3 'DROI' INUMY A3 ;
*
*- Création maillage GEOMETRIQUE
*
* on va faire des parallèlépipèdes réguliers grâce à DALLER
* Les points définissant les quadrangles sont appellés 'SOMMET'
CARRE = 'DALLER' DRBAS DRGAU DRHAU DRDRO 'PLAN' ;
*
*- Création maillage HYBRIDE y compris sous-objets (cond. limites)
* ce sont les maillages avec les points sommets + les points faces +
* les points centres.
 
* on crée celui du maillage total
DOMHYB = 'CHANGER' CARRE 'QUAF' ;
* et ceux des droites sur lesquelles on voudra mettre des conditions
* aux limites (en fait les quatre côtés, ici) :
DOMGAU = 'CHANGER' DRGAU 'QUAF' ;
DOMDRO = 'CHANGER' DRDRO 'QUAF' ;
DOMHAU = 'CHANGER' DRHAU 'QUAF' ;
DOMBAS = 'CHANGER' DRBAS 'QUAF' ;
 
* On a créé à chaque fois les points faces et centres.
* Comme CARRE, DRGAU, DRDRO, DRHAU et DRBAS avaient des segments en
* commun, les opérations précédentes ont créé beaucoup de points en
* double. Il faut les éliminer pour qu'on parle bien de la même chose.
 
* critère d'élimination (inférieur à la plus petite dimension d'arête
ELI0  = L / INUMX / 10.D0   ;
* élimination proprement dite :
'ELIMINATION' ( DOMHYB ET DOMGAU ET DOMDRO ET DOMHAU ET DOMBAS) ELI0 ;
*                                                 ----------------
*                                                 = MODELISATION =
*                                                 ----------------
* Définition du modèle physique attaché au domaine d'étude :
* Il servira d'une part pour la résolution,
* d'autre part pour extraire les points faces, centres, et autres
* informations utiles
MODHYB = 'MODELE' DOMHYB 'DARCY' 'ISOTROPE'  ;
 
* On voudra aussi disposer des points faces des frontières,
* c'est-à-dire des points centres des lignes constituant ces
* frontières. Pour cela, on attache à ces lignes une modélisation,
* inutile en soi, mais qui dit implicitement (comme c'est un modèle
* 'DARCY'), que l'on travaille en éléments finis mixte hybrides,
* et donc qu'on aura besoin, à un moment ou à un autre, des points
* centres, faces, et autres :
MODGAU = 'MODELE' DOMGAU 'DARCY' 'ISOTROPE'  ;
MODDRO = 'MODELE' DOMDRO 'DARCY' 'ISOTROPE'  ;
MODHAU = 'MODELE' DOMHAU 'DARCY' 'ISOTROPE'  ;
MODBAS = 'MODELE' DOMBAS 'DARCY' 'ISOTROPE'  ;
 
* On peut maintenant extraire via ces modèles, les points centre des
* maillages 'QUAF' des frontières, qui sont en fait les points faces
* en limite du domaine :
CEDRO  = 'DOMA' MODDRO 'CENTRE' ;
CEHAU  = 'DOMA' MODHAU 'CENTRE' ;
CEBAS  = 'DOMA' MODBAS 'CENTRE' ;
CEGAU  = 'DOMA' MODGAU 'CENTRE' ;
 
* Les autres petites informations, sur le domaine entier :
* Maillage des points centres :
HYCEN = 'DOMA' MODHYB 'CENTRE' ;
* Maillage des points faces :
HYFAC = 'DOMA' MODHYB 'FACE' ;
* Longueur des arêtes (on est en 2D) :
HYSUR = 'DOMA' MODHYB 'SURFACE'  ; 
* Surface des mailles :
HYVOL = 'DOMA' MODHYB 'VOLUME'  ; 
* Vecteurs normaux aux faces, champ-points à autant de composantes que
* la dimension de l'espace : 'UX' , 'UY' (, 'UZ').
HYNOR = 'DOMA' MODHYB 'NORMALE'  ;
*
*                             ------------------------------------
*                             = Tenseur de perméabilité isotrope =
*                             ------------------------------------
* La grandeur physique correspondant au laplacien. C'est ici la
* perméabilité hydraulique.
VK    = 0.75D0               ;
* On la stoque dans un objet de type matériau
MATI3 = 'MATERIAU' MODHYB 'K' VK ;
*                                          -----------------------
*                                          = Solution analytique =
*                                          -----------------------
* Elle va servir à deux choses : définir les conditions aux limites,
* et vérifier la précision de la résolution.
*
*                                        Trace de charge aux faces
*      La charge étant linéaire, Th moyen est au centre de la face
*
XX YY = 'COOR' (DOMA MODHYB 'FACE')    ;
HG    = 100.D0                 ;
HD    =  20.D0                 ;
DH    = HD - HG                ;
PANAF = (XX - L) * DH / L + HD ;
*
*                                  Charge aux centres des éléments
*  La charge étant linéaire, la moyenne est au centre de l'élément
*
XXC YYC = 'COOR' (DOMA MODHYB 'CENTRE' )  ;
PANAC   = (XXC - L) * DH / L + HD ;
*
*                              Composantes scalaires de la vitesse
*
VXANA = ( -1.D0 ) * VK * DH / L ;
VYANA = 0.D0                    ;
*
*                                 Vitesse aux faces et aux centres
*
VANAF = 'MANU' 'CHPO' (DOMA MODHYB 'FACE')   2 'VX' VXANA 'VY' VYANA ;
VANAC = 'MANU' 'CHPO' (DOMA MODHYB 'CENTRE') 2 'VX' VXANA 'VY' VYANA ;
*
*                                        -------------------------
*                                        = matrice Masse HYBride =
*                                        -------------------------
* C'est l'inverse du tenseur de perméabilité. Il ne faut donc jamais
* avoir mis une valeur de perméabilité nulle.
CND1A = 'MHYB' MODHYB  MATI3 ;
*                                        -------------------------
*                                        = MAtrice globale en TH =
*                                        -------------------------
* La matrice correspondant à div (K 'GRAD').
* C'est un objet de type RIGIDITE.
HND1A = 'MATP' MODHYB  CND1A ;
*                                                 ----------------
*                                                 =  TH imposée  =
*                                                 ----------------
* Les conditions aux limites de Dirichlet portent sur la charge.
* On est en EFMH, donc ce sont les valeurs aux faces qu'on utilise.
* On les appelle trace de charge. C'est un champ-point appuyé aux
* face, de nom de composante 'TH'.
 
* Contribution au membre de gauche pour chaque sous-partie de la frontière
BBGAU = 'BLOQ' CEGAU 'TH' ;
BBDRO = 'BLOQ' CEDRO 'TH' ;
BBHAU = 'BLOQ' CEHAU 'TH' ;
BBBAS = 'BLOQ' CEBAS 'TH' ;
* Contribution au membre de droite pour chaque sous-partie de la frontière
TTIMP = 'REDU' PANAF CEGAU ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
EEGAU = 'DEPI' BBGAU TTIM2          ;
*
TTIMP = 'REDU' PANAF CEDRO ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
EEDRO = 'DEPI' BBDRO TTIM2          ;
*
TTIMP = 'REDU' PANAF CEBAS ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
EEBAS = 'DEPI' BBBAS TTIM2          ;
*
TTIMP = 'REDU' PANAF CEHAU ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
EEHAU = 'DEPI' BBHAU TTIM2          ;
*                                                 ----------------
*                                                 = Flux imposé  =
*                                                 ----------------
* Conditions aux limites de Neuman.
* En discrétisation EFMH comme ici, ce sont des débits intégrés aux
* faces. C'est un champ-point de composante 'FLUX'
 
* Il faut tenir compte de l'orientation des normales aux faces,
* toujours sortantes aux frontières du domaine. Une valeur de débit
* positive correspond donc à un flux sortant.
 
* flux rentrant à gauche - (- K 'GRAD' h) :
FG = VK * DH / L     ;
* champ de débit imposé à gauche, constant sur les faces de la
* frontière de gauche (il faut multiplier le flux par la longueur des
* arêtes pour obtenir le débit intégré) :
FGAU  = 'MANU' 'CHPO' CEGAU 1 'FLUX' FG   ;
LGAU  = 'REDU' HYSUR  CEGAU               ;
FLGAU = FGAU * LGAU                                ;
* on en fait un champ de nature discret (diffus conviendrait tout aussi bien) :
* pour que sa nature ne reste pas 'INDETERMINE'. Cela permettra de le
* concaténer avec d'autres champs par la suite.
FLGAU = 'CHAN' FLGAU 'ATTRIBUT' 'NATURE' 'DISCRET' ;
*
* idem pour le flux sortant à droite :
FD = ( -1.D0 ) * FG  ;
FDRO  = 'MANU' 'CHPO' CEDRO 1 'FLUX' FD   ;
LDRO  = 'REDU' HYSUR  CEDRO               ;
FLDRO = FDRO * LDRO                                ;
FLDRO = 'CHAN' FLDRO 'ATTRIBUT' 'NATURE' 'DISCRET' ;
*                                                 ----------------
*                                                 =  Assemblage  =
*                                                 ----------------
* Création des membres de gauche et de droite de l'équation "A TH = B"
* pour les trois jeux de conditions aux limites choisis.
* CCC# correspond à A, FFF# correspond à B.
* Les CL de Neuman n'ont pas de contribution au terme de gauche.
 
* CL de Dirichlet partout :
CCC1  = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBDRO 'ET' BBGAU ;
FFF1  =            EEHAU 'ET' EEBAS 'ET' EEDRO 'ET' EEGAU ;
* CL de Dirichlet à gauche et à droite,
* et flux nul imposé en haut et en bas (c'est ce qui se passe par
* défaut, donc il suffit de ne rien mettre) :
CCC2  = HND1A 'ET' BBDRO 'ET' BBGAU ;
FFF2  =            EEDRO 'ET' EEGAU ;
* CL de Dirichlet en hau et en bas,
* et flux imposés à gauche et à droite :
CCC3  = HND1A 'ET' BBHAU 'ET' BBBAS                       ;
FFF3  =            EEHAU 'ET' EEBAS 'ET' FLDRO 'ET' FLGAU ;
*                                                -----------------
*                                                = Résolution TH =
*                                                -----------------
* On résoud les trois systèmes d'équations (CL comprises)
* On obtient en sortie la trace de charge, valeur moyenne de la charge
* sur chaque arête.
CHTER1 = 'RESO' CCC1 FFF1 ;
CHTER2 = 'RESO' CCC2 FFF2 ;
CHTER3 = 'RESO' CCC3 FFF3 ;
*                                                 ----------------
*                                                 = Résolution H =
*                                                 ----------------
* On en déduit la charge moyenne par élément, sous la forme d'un
* champ-point de composante 'H' appuyé aux centres des éléments.
PCEN1 = 'HYBP' MODHYB CND1A CHTER1 ;
PCEN2 = 'HYBP' MODHYB CND1A CHTER2 ;
PCEN3 = 'HYBP' MODHYB CND1A CHTER3 ;
*                                                 ----------------
*                                                 = Résolution V =
*                                                 ----------------
* On en déduit aussi pour chaque jeu de conditions aux limites,
* les débits intégrés sur chaque arête, sous la forme d'un champ-point
* de composante 'FLUX' appuyé aux centre des faces :
QFACE1 = 'HDEB' MODHYB CND1A PCEN1 CHTER1 ;
* La vitesse moyenne dans chaque élément, sous la forme d'un champ-point
* de composantes 'VX', 'VY' (,'VZ') appuyé aux centres des éléments :
VCENT1 = 'HVIT' MODHYB QFACE1 ;
* La vitesse moyenne à travers chaque face, sous la forme d'un champ-point
* de composante 'SCAL' appuyé aux centres des faces :
VFACE1 = ('EXCO' QFACE1 'FLUX' 'SCAL') * HYNOR / HYSUR;
 
* de même pour le 2e jeux de CLs
QFACE2 = 'HDEB' MODHYB  CND1A  PCEN2 CHTER2 ;
VCENT2 = 'HVIT' MODHYB  QFACE2         ;
VFACE2 = ('EXCO' QFACE2 'FLUX' 'SCAL') * HYNOR / HYSUR;
 
* de même pour le 3e jeux de CLs
QFACE3 = 'HDEB' MODHYB  CND1A PCEN3 CHTER3  ;
VCENT3 = 'HVIT' MODHYB  QFACE3       ;
VFACE3 = ('EXCO' QFACE3 'FLUX' 'SCAL') * HYNOR / HYSUR;
 
*                                                -----------------
*                                                = Calcul ERREUR =
*                                                -----------------
* Comparaison calcul analytique, calcul numérique :
 
* Erreur relative en trace de charge (moyenne aux faces)
* Erreur relative en charge (moyenne aux centres des éléments)
*
* Jeu de CL n°1
ERRTP1 = 'EXCO' CHTER1 'TH' 'SCAL' ;
ERRTP1 = ERRTP1 - PANAF / PANAF    ;
ERRTP1 = 'ABS' ERRTP1              ;
ERRP1  = 'EXCO' PCEN1 'H' 'SCAL'   ;
ERRP1  = ERRP1 - PANAC / PANAC     ;
ERRP1  = 'ABS' ERRP1               ;
* Jeu de CL n°2
ERRTP2 = 'EXCO' CHTER2 'TH' 'SCAL' ;
ERRTP2 = ERRTP2 - PANAF / PANAF    ;
ERRTP2 = 'ABS' ERRTP2              ;
ERRP2  = 'EXCO' PCEN2 'H' 'SCAL'   ;
ERRP2  = ERRP2 - PANAC / PANAC     ;
ERRP2  = 'ABS' ERRP2               ;
* Jeu de CL n°3
ERRTP3 = 'EXCO' CHTER3 'TH' 'SCAL' ;
ERRTP3 = ERRTP3 - PANAF / PANAF    ;
ERRTP3 = 'ABS' ERRTP3              ;
ERRP3  = 'EXCO' PCEN3 'H' 'SCAL'   ;
ERRP3  = ERRP3 - PANAC / PANAC     ;
ERRP3  = 'ABS' ERRP3               ;
*
* Erreur relative en vitesse aux centres des éléments
MOT1 = 'MOTS' 'VX' 'VY'             ;
VDVD = 'PSCA' VANAC VANAC MOT1 MOT1 ;
* Jeu de CL n°1
VD1  = VANAC - VCENT1           ;
VC1  = 'PSCA' VD1 VD1 MOT1 MOT1 ;
SDC1 = 'ABS' ( VC1 / VDVD )     ;
SDC1 = SDC1 '**' 0.5D0          ;
* Jeu de CL n°2
VD2  = VANAC - VCENT2           ;
VC2  = 'PSCA' VD2 VD2 MOT1 MOT1 ;
SDC2 = 'ABS' ( VC2 / VDVD )     ;
SDC2 = SDC2 '**' 0.5D0          ;
* Jeu de CL n°3
VD3  = VANAC - VCENT3           ;
VC3  = 'PSCA' VD3 VD3 MOT1 MOT1 ;
SDC3 = 'ABS' ( VC3 / VDVD )     ;
SDC3 = SDC3 '**' 0.5D0          ;
*
*                                                 -------------------
*                                                 = Tracé resultats =
*                                                 -------------------
'SI' ('NEG' GRAPH 'N') ;
*
*- Transformation des quantités aux centres en MCHAML
*
ERRP1 = 'KCHA' MODHYB 'CHAM' ERRP1 ;
ERRP2 = 'KCHA' MODHYB 'CHAM' ERRP2 ;
ERRP3 = 'KCHA' MODHYB 'CHAM' ERRP3 ;
SDS1  = 'KCHA' MODHYB 'CHAM' SDC1  ;
SDS2  = 'KCHA' MODHYB 'CHAM' SDC2  ;
SDS3  = 'KCHA' MODHYB 'CHAM' SDC3  ;
*
* Dans chaque cas on trace
*   L'erreur relative sur la charge moyenne aux centres des éléments
*   L'erreur relative sur la Vitesse aux centres des éléments
*
'TITR' 'darcy1/1 : Erreur relative sur la charge'  ;
'TRAC' MODHYB ERRP1                                ;
'TITR' 'darcy1/1 : Erreur relative sur la vitesse' ;
'TRAC' MODHYB SDS1                                 ;
*
'TITR' 'darcy1/2 : Erreur relative sur la charge'  ;
'TRAC' MODHYB ERRP2                                ;
'TITR' 'darcy1/2 : Erreur relative sur la vitesse' ;
'TRAC' MODHYB SDS2                                 ;
*
'TITR' 'darcy1/3 : Erreur relative sur la charge'  ;
'TRAC' MODHYB ERRP3                                ;
'TITR' 'darcy1/3 : Erreur relative sur la vitesse' ;
'TRAC' MODHYB SDS3                                 ;
*
'FINSI' ;
*
*                                                 -------------------
*                                                 = Gestion ERREURS =
*                                                 -------------------
MAXTP1 = 'MAXI' ERRTP1 ;
MAXTP2 = 'MAXI' ERRTP2 ;
MAXTP3 = 'MAXI' ERRTP3 ;
MAXP1  = 'MAXI' ERRP1  ;
MAXP2  = 'MAXI' ERRP2  ;
MAXP3  = 'MAXI' ERRP3  ;
MAXV1  = 'MAXI' SDC1   ;
MAXV2  = 'MAXI' SDC2   ;
MAXV3  = 'MAXI' SDC3   ;
*
'SAUT' 'PAGE'                                                    ;
'SAUT' 2 'LIGNE'                                                 ;
'MESS' '                           ERREURS RELATIVES           ' ;
'SAUT' 1 'LIGNE'                                                 ;
'MESS' ' cas test        TH                  H                V' ;
'SAUT' 1 'LIGNE'                                                 ;
'MESS' ' numero1   ' maxtp1  '   ' maxp1 '   ' maxv1             ;
'SAUT' 1 'LIGNE'                                                 ;
'MESS' ' numero2   ' maxtp2  '   ' maxp2 '   ' maxv2             ;
'SAUT' 1 'LIGNE'                                                 ;
'MESS' ' numero3   ' maxtp3  '   ' maxp3 '   ' maxv3             ;
'SAUT' 2 'LIGNE'                                                 ;
*
EPS0 = 1.E-13 ;
LOG1 = MAXTP1 > EPS0 ; LOG2 = MAXTP2 > EPS0 ; LOG3 = MAXTP3 > EPS0 ;
LOG4 = MAXP1  > EPS0 ; LOG5 = MAXP2  > EPS0 ; LOG6 = MAXP3  > EPS0 ;
LOG7 = MAXV1  > EPS0 ; LOG8 = MAXV2  > EPS0 ; LOG9 = MAXV3  > EPS0 ;
LTP0 = LOG1 'OU' LOG2 'OU' LOG3 ;
LP0  = LOG4 'OU' LOG5 'OU' LOG6 ;
LV0  = LOG7 'OU' LOG8 'OU' LOG9 ;
L0   = LTP0 'OU' LP0  'OU' LV0  ;
 
* Compte-rendu de fin et sortie
'SI' L0        ;
   'ERREUR' 5  ;
'SINON'        ;
   'ERREUR' 0  ;
'FINSI'        ;
*
'FIN' ;