Accueil | Cast3M

* fichier :  darcy1.dgibi
************************************************************************
************************************************************************
*****************  CAS TEST : darcy1.dgibi  *************************
*
GRAPH = faux   ;
'SAUT' 'PAGE' ;
*
*-------------------------------------------------------------------
* TEST DARCY9
*
* C'est une reprise de DARCY1
* avec utilisation de conditions aux limites d'inégalité
*
*-  1 premier calcul comme DARCY1 pour vérifier que les modifs conservent
*     le résultat
*-  1 deuxième calcul avec CL d'inégalité et terme source pour l'exemple.
*
* Les modifications sont repérées par ------>
*
*  . multiplication par -1 de la matrice issue de MATP,
*                          du terme source issu de SMTP
*                          des débits au second membre,
*  . on ne touche pas aux multiplicateurs de Lagrange (BLOQ et DEPI)
*
*-------------------------------------------------------------------
*
'SAUT' 'PAGE' ;
*
*- Options générales de calcul.
*
'TITR' 'EFMH DARCY CL inégalité : darcy9.dgibi' ;
* Dimension 2, éléments à générer quadrangulaires ou triangulaires :
'OPTI' 'DIME' 2 'ELEM' 'TRI3'       ;
* Pas d'affichage à l'écran des lignes de commandes lues :
'OPTI' 'ECHO' 0           ;
*                                                     ------------
*                                                     = 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 ;
 
* ==============
* Premier Calcul
* ==============
 
* Cas test 1D classique pour valider les modifications sur solution
* analytique
 
*                                          = Solution analytique =
*                                          -----------------------
XX YY = 'COOR' (DOMA MODHYB 'FACE')    ;
HG    = 100.D0     ;
HD    =  20.D0     ;
DH    = HD - HG    ;
PANAF = (XX - L) * DH / L + HD ;
 
XXC YYC = 'COOR' (DOMA MODHYB 'CENTRE' )  ;
PANAC   = (XXC - L) * DH / L + HD ;
 
VXANA = ( -1.D0 ) * VK * DH / L ;
VYANA = 0.D0        ;
 
VANAF = 'MANU' 'CHPO' (DOMA MODHYB 'FACE')   2 'VX' VXANA 'VY' VYANA ;
VANAC = 'MANU' 'CHPO' (DOMA MODHYB 'CENTRE') 2 'VX' VXANA 'VY' VYANA ;
*
*                                             = Modèle numérique =
*                                             --------------------
CND1A = 'MHYB' MODHYB  MATI3 ;
*------> on multiplie par -1. pour qu'elle soit définie positive
HND1A = -1.*('MATP' MODHYB  CND1A) ;
 
* CLs de Dirichlet
* 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          ;
 
* CL de Neuman
FG    = VK * DH / L     ;
 
FGAU  = 'MANU' 'CHPO' CEGAU 1 'FLUX' FG   ;
LGAU  = 'REDU' HYSUR  CEGAU   ;
FLGAU = FGAU * LGAU           ;
FLGAU = 'CHAN' FLGAU 'ATTRIBUT' 'NATURE' 'DISCRET' ;
 
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
*------> les flux sont multipliés par -1. pour correspondre à HND1A
CCC3  = HND1A 'ET' BBHAU 'ET' BBBAS  ;
FFF3  =            EEHAU 'ET' EEBAS 'ET' (-1.* (FLDRO 'ET' FLGAU)) ;
 
CHTER3 = 'RESO' CCC3 FFF3 ;
 
PCEN3  = 'HYBP' MODHYB CND1A CHTER3 ;
QFACE3 = 'HDEB' MODHYB CND1A PCEN3 CHTER3 ;
VCENT3 = 'HVIT' MODHYB QFACE3 ;
VFACE3 = ('EXCO' QFACE3 'FLUX' 'SCAL') * HYNOR / HYSUR;
 
* Calcul ERREUR
ERRTP3 = 'EXCO' CHTER3 'TH' 'SCAL' ;
ERRTP3 = ERRTP3 - PANAF / PANAF    ;
ERRTP3 = 'ABS' ERRTP3  ;
MAXTP3 = 'MAXI' ERRTP3 ;
 
ERRP3  = 'EXCO' PCEN3 'H' 'SCAL'   ;
ERRP3  = ERRP3 - PANAC / PANAC     ;
ERRP3  = 'ABS' ERRP3   ;
MAXP3  = 'MAXI' ERRP3  ;
 
MOT1   = 'MOTS' 'VX' 'VY' ;
VDVD   = 'PSCA' VANAC VANAC MOT1 MOT1 ;
VD3    = VANAC - VCENT3           ;
VC3    = 'PSCA' VD3 VD3 MOT1 MOT1 ;
SDC3   = 'ABS' ( VC3 / VDVD )     ;
SDC3   = SDC3 '**' 0.5D0          ;
MAXV3  = 'MAXI' SDC3   ;
 
'SI' (GRAPH) ;
  'TITRE' 'Charge pour pb uniforme 1D' ;
  'TRACER' ('KCHA' PCEN3 modhyb 'CHAM') modhyb ;
'FINSI' ;
 
* ===============
* Deuxième Calcul
* ===============
 
* Cas avec CL d'inégalité.
 
* Nouvelles CL
* ------------
* on avait Hgauche=HG=100, on met Hgauche >= HG
BBGAU = 'BLOQ' CEGAU 'MINIMUM' 'TH' ;
TTIMP = 'REDU' PANAF CEGAU ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
EEGAU = 'DEPI' BBGAU TTIM2          ;
*
* on avait Hdroite=HD=20 , on met Hdroite <= HD
BBDRO = 'BLOQ' CEDRO 'MAXIMUM' 'TH' ;
TTIMP = 'REDU' PANAF CEDRO ;
TTIM2 = 'EXCO' TTIMP 'SCAL' 'TH'    ;
EEDRO = 'DEPI' BBDRO TTIM2          ;
 
* Résolution et post-traitement
CCC4   = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBGAU 'ET' BBDRO ;
FFF4   =            EEHAU 'ET' EEBAS 'ET' EEGAU 'ET' EEDRO ;
CHTER4 = 'RESO' CCC4 FFF4 ;
 
PCEN4  = 'HYBP' MODHYB CND1A CHTER4 ;
QFACE4 = 'HDEB' MODHYB CND1A PCEN4 CHTER4 ;
VCENT4 = 'HVIT' MODHYB QFACE4 ;
VFACE4 = ('EXCO' QFACE4 'FLUX' 'SCAL') * HYNOR / HYSUR;
 
'SI' (GRAPH) ;
  'TITRE' 'Charge pour pb uniforme 1D avec CL unilatérales' ;
  'TRACER' ('KCHA' PCEN4 modhyb 'CHAM') modhyb ;
'FINSI' ;
 
 
* Calcul ERREUR
ERRTP4 = 'EXCO' CHTER4 'TH' 'SCAL' ;
ERRTP4 = ERRTP4 - PANAF / PANAF    ;
ERRTP4 = 'ABS' ERRTP4  ;
MAXTP4 = 'MAXI' ERRTP4 ;
 
ERRP4  = 'EXCO' PCEN4 'H' 'SCAL'   ;
ERRP4  = ERRP4 - PANAC / PANAC     ;
ERRP4  = 'ABS' ERRP4   ;
MAXP4  = 'MAXI' ERRP4  ;
 
MOT1   = 'MOTS' 'VX' 'VY' ;
VDVD   = 'PSCA' VANAC VANAC MOT1 MOT1 ;
VD4    = VANAC - VCENT4           ;
VC4    = 'PSCA' VD4 VD4 MOT1 MOT1 ;
SDC4   = 'ABS' ( VC4 / VDVD )     ;
SDC4   = SDC4 '**' 0.5D0          ;
MAXV4  = 'MAXI' SDC4   ;
 
* Troisème Calcul
* ===============
 
* Cas avec CL d'inégalité.
* et un terme source pour illustrer son emploi.
 
* on l'impose au point centre proche du barycentre
pt     = 'BARYCENTRE' carre ;
pts    = 'POINT' hycen 'PROC' pt ;
cs1    = 'MANU' 'CHPO' ('MANU' pts 'POI1') 1 'SOUR' 100. ;
tersc1 = 'KCHT' modhyb 'SCAL' 'CENTRE' 'COMP' 'SOUR' 0. cs1 ;
*------> les termes sources sont multipliés par -1. pour correspondre à HND1A
SMTR   = -1.*('SMTP' modhyb CND1A tersc1) ;
 
* Résolution et post-traitement
CCC5   = HND1A 'ET' BBHAU 'ET' BBBAS 'ET' BBGAU 'ET' BBDRO ;
FFF5   =            EEHAU 'ET' EEBAS 'ET' EEGAU 'ET' EEDRO 'ET' SMTR ;
CHTER5 = 'RESO' CCC5 FFF5 ;
 
PCEN5  = 'HYBP' MODHYB CND1A CHTER5 ;
QFACE5 = 'HDEB' MODHYB CND1A PCEN5 CHTER5 ;
VCENT5 = 'HVIT' MODHYB QFACE5 ;
VFACE5 = ('EXCO' QFACE5 'FLUX' 'SCAL') * HYNOR / HYSUR;
 
'SI' (GRAPH) ;
  'TITRE' 'Charge pour pb uniforme 1D avec CL unilatérales' ;
  'TRACER' ('KCHA' PCEN5 modhyb 'CHAM') modhyb ;
'FINSI' ;
 
* On vérifie que la charge a bien varié dans le bon sens (positif)
vh4    = 'EXTR' pcen4 'H' pts ;
vh5    = 'EXTR' pcen5 'H' pts ;
 
LOG3   = 'NON' (vh5 > vh4) ;
*
*                                                 -------------------
*                                                 = Gestion ERREURS =
*                                                 -------------------
*
'SAUT' 'PAGE'                                                    ;
'SAUT' 2 'LIGNE'                                                 ;
'MESS' '                           ERREURS RELATIVES           ' ;
'SAUT' 1 'LIGNE'                                                 ;
'MESS' ' cas test        TH                  H                V' ;
'SAUT' 1 'LIGNE'                                                 ;
'MESS' ' numero3   ' maxtp3  '   ' maxp3 '   ' maxv3             ;
'SAUT' 1 'LIGNE'                                                 ;
'MESS' ' numero4   ' maxtp4  '   ' maxp4 '   ' maxv4             ;
'SAUT' 2 'LIGNE'                                                 ;
*
EPS0 = 1.E-13 ;
LOG1 = MAXTP3 > EPS0 ; LOG2 = MAXTP4 > EPS0  ;
LOG4 = MAXP3  > EPS0 ; LOG5 = MAXP4  > EPS0  ;
LOG7 = MAXV3  > EPS0 ; LOG8 = MAXV4  > EPS0  ;
LTP0 = LOG1 'OU' LOG2 ;
LP0  = LOG4 'OU' LOG5 ;
LV0  = LOG7 'OU' LOG8 ;
L0   = LTP0 'OU' LP0  'OU' LV0 'OU' LOG3 ;
 
* Compte-rendu d'erreur
'SI' L0        ;
   'ERREUR' 5  ;
'SINON'        ;
   'ERREUR' 0  ;
'FINSI'        ;
 
'FIN' ;