Download equ_chaleur3Dtet.dgibi
* fichier : equ_chaleur3Dtet.dgibi
************************************************************************
* NOM : equ_chaleur3Dtet.dgibi
* ___
*
* DESCRIPTION : Solution stationnaire de l'équation de la chaleur (3D)
* ___________
*
* GEOMETRIE : Un cube de côté 1 (il sera translaté et
* ---------- tourné) maillé avec des tétraèdres.
*
* EQUATIONS :
* ----------
*
* - Equations :
*
* mu Laplacien(T) = 0 avec mu = 1 ici
*
* - Conditions aux limites :
*
* conditions de Dirichlet sur tout le bord :
* on prend la restriction de la solution exacte au bord
*
* - Solution exacte :
*
* T(x,y)= alpha x + beta y + gamma z + eta xy + delta
*
* pour une discrétisation avec de éléments linéaires,
* on prend eta = 0 (car les triangles linéaires
* disponibles ne discrétisent pas de façon exacte le
* terme bilinéaire en xy).
*
* DISCRETISATION : On teste les types d'éléments continus suivants (voir
* la notice NAVI) :
* * 'QUAF' (température au sommet) ;
* * 'LINE' (température au sommet).
* On devrait également tester les discrétisations
* discontinus (linéaires ou Ctes par élément). Ceci
* n'est pas fait.
* le tértaèdre 'MACRO' ne fonctionne pas non plus.
*
* Le maillage est construit avec l'opérateur VOLU, il est perturbé
* par un bruit blanc gaussien. Il est ensuite translaté d'un
* vecteur arbitraire et tourné d'un angle arbitraire.
*
* Opérateurs utilisés : LAPN (option EF implicite)
*
*
* RESOLUTION : Directe (Crout)
*
*
* TESTS EFFECTUES : La solution du problème discrétisé doit être la
* _______________ même que la solution exacte aux erreurs de
* troncature près.
* On impose donc des valeurs faibles pour l'écart en
* norme L2 entre solution exacte et calculée :
* * 1.D-14
*
* LANGAGE : GIBIANE-CASTEM 2000
* AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
* mél : gounand@semt2.smts.cea.fr
************************************************************************
* VERSION : v1, 30/05/00, version initiale
* HISTORIQUE : v1, 30/05/00, création
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
interact= FAUX ;
complet = FAUX ;
graph = FAUX ;
*
'OPTION' 'DIME' 3 'ELEM' CUB8 ;
'OPTION' 'ISOV' 'SULI' ;
'OPTION' 'ECHO' 0 ;
nbisov = 15 ;
'SI' ('NON' interact) ;
'OPTION' 'TRAC' 'PS' ;
'OPTION' 'ISOV' 'LIGNE' ;
'FINSI' ;
*
** Erreur Linfini entre deux Champoints.
*
DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
err = MAXI (vitp1 - vit) 'ABS' ;
FINPROC err ;
*
** Erreur Pseudo L2 entre deux Champoints.
*
DEBPROC CALCERR2 vitp1*'CHPOINT' vit*'CHPOINT' ;
er0 = 'KOPS' vitp1 '-' vit ;
er2 = 'KOPS' er0 '*' er0 ;
err = 0.D0 ;
suppv = 'EXTRAIRE' er2 'MAIL' ;
compv = 'EXTRAIRE' er2 'COMP' ;
nbcomp = 'DIME' compv ;
nptd = 'NBNO' suppv ;
'REPETER' numcomp nbcomp ;
lacomp = 'EXTRAIRE' compv &numcomp ;
'REPETER' numpo nptd ;
lepoi = suppv 'POIN' &numpo ;
err = err '+' ('EXTRAIRE' er2 lacomp lepoi) ;
'FIN' numpo ;
'FIN' numcomp ;
err = err '/' nptd ;
err = err ** 0.5D0 ;
FINPROC err ;
*
* Procédure paramétrée (type de discrétisation)
* renvoyant l'erreur en norme L2 sur la température.
* On calcule une solution exacte de l'équation de Laplace
* (équation de la chaleur) ;
*
'DEBPROC' CALCUL nraff*'ENTIER' typdis*'MOT' ;
*
* titre global pour les dessins
*
titglob = 'CHAINE' ' ; nraff=' nraff ' typdis=' typdis ;
*
* Paramètres physiques
*
mu = 1. ;
* Conditions aux limites
alpha = (** 2. 0.5) ;
beta = (** 3. 0.5) ;
gamma = (** 5. 0.5) ;
'SI' ('EGA' typdis 'QUAF') ;
eta = (** 7. 0.5) ;
'SINON' ;
eta = 0. ;
'FINSI' ;
delta = (** 1.5 0.5) ;
*
* Géométrie
*
pA = 0. 0. 0. ;
pB = 1. 0. 0. ;
pC = 1. 1. 0. ;
pD = 0. 1. 0. ;
pE = 0. 0. 1. ;
pF = 1. 0. 1. ;
pG = 1. 1. 1. ;
pH = 0. 1. 1. ;
*
* Paramètres de la discrétisation de base
*
'SI' complet ;
nAB = 6 ;
nBC = 6 ;
nCD = 6 ;
nDA = 6 ;
nH = 6 ;
'SINON' ;
nAB = 3 ;
nBC = 3 ;
nCD = 3 ;
nDA = 3 ;
nH = 3 ;
'FINSI' ;
*
* Géométrie discrétisée
*
'OPTION' 'ELEM' 'CUB8' ;
bas = 'DROIT' nAB pA pB ;
droite = 'DROIT' nBC pB pC ;
haut = 'DROIT' nCD pC pD ;
gauche = 'DROIT' nDA pD pA ;
spourt = bas 'ET' droite 'ET' haut 'ET' gauche ;
sbas = 'SURFACE' spourt 'PLAN' ;
mt = 'VOLUME' sbas 'TRAN' nH pE ;
*mt = 'ORIENTER' mt ;
f1 f2 f3 = 'FACE' mt ;
f1 = 'CHANGER' f1 'TRI3' ;
f2 = 'CHANGER' f2 'TRI3' ;
f3 = 'CHANGER' f3 'TRI3' ;
ss = (f1 'ET' f2 'ET' f3) ;
'OPTION' 'ELEM' 'TET4' ;
mt = 'VOLUME' ss ;
pourtour = 'ENVELOPPE' mt ;
*
* Rotation et translation aditionnelle + bruit blanc
* + raffinement
*
vtran = ((** 2 0.5) (* -1 (** 3 0.5)) (** 5 0.5)) ;
orig = (0.D0 0.D0 0.D0) ;
orig2 = (-0.5D0 0.5D0 0.5D0) ;
arot = (** 360. 0.5) ;
*arot = 0.D0 ;
*arot = 90.0D0 ;
ladens = (** ('/' ('MESURE' mt) ('NBEL' mt)) ('/' 1.D0 3.D0)) ;
*'MESSAGE' ('CHAINE' 'ladens=' ladens) ;
ampbruit = (0.02 * ladens) ;
*
* On rajoute le bruit
*
pmt = 'CHANGER' mt 'POI1' ;
pcnt= 'CHANGER' pourtour 'POI1' ;
inmt= 'DIFF' pmt pcnt;
brbl1 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
brbl2 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
brbl = 'ET' ('NOMC' 'UX' brbl1) ('NOMC' 'UY' brbl2) ;
mt = 'PLUS' mt brbl;
*
* On raffine nraff fois
*
'SI' (nraff > 0) ;
'REPETER' bcl nraff ;
mt = 'CHANGER' mt 'QUADRATIQUE' ;
$mt = 'DOMA' mt 'MACRO' ;
mt = ($mt . 'MAILLAGE') ;
pourtour = 'CHANGER' pourtour 'QUADRATIQUE' ;
$pourtou = 'DOMA' pourtour 'MACRO' ;
pourtour = ($pourtou . 'MAILLAGE') ;
'MENAGE' ;
'FIN' bcl ;
'FINSI' ;
*
* Eventuellement, on trace le résultat
*
'SI' graph ;
titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
' NBELEM=' ('NBEL' mt) titglob ;
'TRACER' mt 'NOEUD' 'TITRE' titgeo ;
'FINSI' ;
*
* Modèle
*
_mt = 'CHANGER' mt 'QUAF' ;
_pourt = 'CHANGER' pourtour 'QUAF' ;
'ELIMINATION' (_mt 'ET' _pourt) 1.D-6 ;
$mt = 'MODELISER' _mt 'NAVIER_STOKES' typdis ;
$pourt = 'MODELISER' _pourt 'NAVIER_STOKES' typdis ;
*
* Solution exacte trilinéaire
*
xx yy zz = 'COORDONNEE' ('DOMA' $mt 'SOMMET') ;
solex = '+' ('+' ('*' xx alpha) ('*' yy beta)) ('*' zz gamma) ;
*solex = '+' solex ('*' ('*' ('*' xx yy) zz) eta) ;
solex = '+' solex ('*' ('*' xx yy) eta) ;
solex = '+' solex delta ;
solex = 'KCHT' $mt 'SCAL' 'SOMMET' solex ;
*
* On tourne et on translate
*
tmt tpourt _tmt _tpourt solex =
'TOURNER' mt pourtour _pourt _mt solex arot orig orig2 ;
tmt tpourt _tmt _tpourt solex =
'PLUS' tmt tpourt _tpourt _tmt solex vtran ;
'ELIMINATION' 1.D-6 (_tmt 'ET' _tpourt 'ET' tmt 'ET' tpourt
'ET' ('EXTRAIRE' solex 'MAIL')) ;
$tmt = 'MODELISER' _tmt 'NAVIER_STOKES' typdis ;
$tpourt = 'MODELISER' _tpourt 'NAVIER_STOKES' typdis ;
solex = 'KCHT' $tmt 'SCAL' 'SOMMET' solex ;
mtmt = 'DOMA' $tmt 'MAILLAGE' ;
cmtmt = 'DOMA' $tpourt 'MAILLAGE' ;
*
* Mise en place du calcul numérique
*
* équation de Laplace
*
rv = 'EQEX' $tmt 'ITMA' 1
'OPTI' 'EF' 'IMPL'
'ZONE' $tmt 'OPER' 'LAPN' mu 'INCO' 'TN' 'TN' ;
*
* conditions aux limites
*
csolex = 'REDU' solex cmtmt ;
rv = 'EQEX' rv
'CLIM' 'TN' 'TIMP' cmtmt csolex ;
* ________________
*
* INITIALISATION
* ________________
*
rv . 'INCO' = 'TABLE' 'INCO' ;
rv . 'INCO' . 'TN' = 'KCHT' $tmt 'SCAL' 'SOMMET' 0. ;
*
rv . 'NITER' = 1 ;
rv . 'METHINV' . 'TYPINV' = 1 ;
rv . 'METHINV' . 'IMPINV' = 2 ;
* __________
*
* CALCUL
* __________
*
EXEC rv ;
*
* Résultats
*
'SI' graph ;
*
* solutions exactes
*
tn = solex ;
titt = 'CHAINE' 'Température exacte' titglob ;
'TRACER' tn tmt tpourt nbisov 'TITRE' titt ;
*
* solutions calculées
*
* tn = rv . 'INCO' . 'TN' ;
* titt = 'CHAINE' 'Température calculée' titglob ;
* 'TRACER' tn tmt tpourt nbisov 'TITRE' titt ;
'FINSI' ;
*
* Calcul des erreurs par rapport à la solution analytique
*
tn = rv . 'INCO' . 'TN' ;
errlit = CALCERR tn solex ;
errl2t = CALCERR2 tn solex ;
'MESSAGE' '-------------------------------------------------' ;
'MESSAGE' ('CHAINE' 'Erreur en norme LINF' titglob ' : ' errlit) ;
'MESSAGE' ('CHAINE' 'Erreur en norme L2 ' titglob ' : ' errl2t) ;
'MESSAGE' '-------------------------------------------------' ;
'FINPROC' errl2t ;
* Fin de la procédure de calcul *
*___________________________________________________________*
*
*-----------------------------------------------------------
* Paramètres du calcul
*
* ldiscr : type d'éléments à tester.
'SI' complet ;
ldiscr = 'MOTS' 'QUAF' 'LINE' ;
'SINON' ;
ldiscr = 'MOTS' 'QUAF' 'LINE' ;
'FINSI' ;
*
*-----------------------------------------------------------*
* Boucles sur les différents paramètres du calcul *
ok = VRAI ;
*
* Boucle sur les discrétisations
*
'REPETER' idiscr ('DIME' ldiscr) ;
discr = 'EXTRAIRE' ldiscr &idiscr ;
errl2t = CALCUL 0 discr ;
ok = ('ET' ok (errl2t '<' 1.D-14)) ;
'SI' ('ET' ('NON' complet) ('NON' ok)) ;
'MESSAGE' 'On a dépassé les marges d_erreur :' ;
'MESSAGE' ' * 1.D-14 pour une méthode directe' ;
'MESSAGE' ' * 1.D-12 pour une méthode itérative' ;
* 'ERREUR' 5 ;
'FINSI' ;
'FIN' idiscr ;
'SI' interact ;
'OPTION' 'DONN' 5 ;
'FINSI' ;
'MESSAGE' 'Tout s_est bien passé' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer