* 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' ;