* 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' ;
© Cast3M 2003 - All rights reserved.
Disclaimer