Télécharger la plaquette Cast3M


fr en

Télécharger equ_chaleur3Dtet.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  equ_chaleur3Dtet.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ************************************************************************
  5. * NOM         : equ_chaleur3Dtet.dgibi
  6. * ___
  7. *
  8. * DESCRIPTION : Solution stationnaire de l'équation de la chaleur  (3D)
  9. * ___________  
  10. *
  11. *               GEOMETRIE : Un cube de côté 1 (il sera translaté et
  12. *               ----------  tourné) maillé avec des tétraèdres.
  13. *
  14. *               EQUATIONS :
  15. *               ----------
  16. *
  17. *               - Equations :
  18. *
  19. *                  mu Laplacien(T) = 0 avec mu = 1 ici
  20. *                  
  21. *               - Conditions aux limites :
  22. *
  23. *                  conditions de Dirichlet sur tout le bord :
  24. *                  on prend la restriction de la solution exacte au bord
  25. *
  26. *               - Solution exacte :
  27. *
  28. *                  T(x,y)= alpha x + beta y + gamma z + eta xy + delta
  29. *
  30. *                  pour une discrétisation avec de éléments linéaires,
  31. *                  on prend eta = 0 (car les triangles linéaires
  32. *                  disponibles ne discrétisent pas de façon exacte le
  33. *                  terme bilinéaire en xy).
  34. *
  35. * DISCRETISATION : On teste les types d'éléments continus suivants (voir
  36. *                  la notice NAVI) : 
  37. *                  * 'QUAF' (température au sommet) ;
  38. *                  * 'LINE' (température au sommet).
  39. *                  On devrait également tester les discrétisations
  40. *                  discontinus (linéaires ou Ctes par élément). Ceci
  41. *                  n'est pas fait.
  42. *                  le tértaèdre 'MACRO' ne fonctionne pas non plus.
  43. *
  44. *       Le maillage est construit avec l'opérateur VOLU, il est perturbé
  45. *       par un bruit blanc gaussien. Il est ensuite translaté d'un
  46. *       vecteur arbitraire et tourné d'un angle arbitraire.
  47. *
  48. * Opérateurs utilisés : LAPN (option EF implicite)
  49. *
  50. *
  51. * RESOLUTION :     Directe (Crout)
  52. *
  53. *
  54. * TESTS EFFECTUES : La solution du problème discrétisé doit être la 
  55. * _______________   même que la solution exacte aux erreurs de
  56. *                   troncature près.
  57. *                   On impose donc des valeurs faibles pour l'écart en
  58. *                   norme L2 entre solution exacte et calculée :
  59. *                   * 1.D-14
  60. *
  61. * LANGAGE     : GIBIANE-CASTEM 2000
  62. * AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  63. *               mél : gounand@semt2.smts.cea.fr
  64. ************************************************************************
  65. * VERSION    : v1, 30/05/00, version initiale
  66. * HISTORIQUE : v1, 30/05/00, création
  67. * HISTORIQUE :
  68. ************************************************************************
  69. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  70. * en cas de modification de ce sous-programme afin de faciliter
  71. * la maintenance !
  72. ************************************************************************
  73. interact= FAUX ;
  74. complet = FAUX ;
  75. graph   = FAUX ;
  76. *
  77. 'OPTION' 'DIME' 3 'ELEM' CUB8  ;
  78. 'OPTION' 'ISOV' 'SULI' ;
  79. 'OPTION' 'ECHO' 0 ;
  80. nbisov = 15 ;
  81. 'SI' ('NON' interact) ;
  82.   'OPTION' 'TRAC' 'PS' ;
  83.   'OPTION' 'ISOV' 'LIGNE' ;
  84. 'FINSI' ;  
  85. *
  86. ** Erreur Linfini entre deux Champoints.
  87. *
  88. DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
  89.   err = MAXI (vitp1 - vit) 'ABS' ;
  90. FINPROC err ;
  91. *
  92. ** Erreur Pseudo L2 entre deux Champoints.
  93. *
  94. DEBPROC CALCERR2 vitp1*'CHPOINT' vit*'CHPOINT' ;
  95.   er0 = 'KOPS' vitp1 '-' vit ;
  96.   er2 = 'KOPS' er0 '*' er0 ;
  97.   err = 0.D0 ;
  98.   suppv  = 'EXTRAIRE' er2 'MAIL' ;
  99.   compv  = 'EXTRAIRE' er2 'COMP' ;
  100.   nbcomp = 'DIME' compv ;
  101.   nptd   = 'NBNO' suppv ;
  102.   'REPETER' numcomp nbcomp ;
  103.   lacomp = 'EXTRAIRE' compv &numcomp ;
  104.     'REPETER' numpo nptd ;
  105.       lepoi = suppv 'POIN' &numpo ;
  106.       err = err '+' ('EXTRAIRE' er2 lacomp lepoi) ;
  107.     'FIN' numpo ;
  108.   'FIN' numcomp ;
  109.   err = err '/' nptd ;
  110.   err = err ** 0.5D0 ;  
  111. FINPROC err ;
  112. *
  113. * Procédure paramétrée (type de discrétisation)
  114. * renvoyant l'erreur en norme L2 sur la température.
  115. * On calcule une solution exacte de l'équation de Laplace
  116. * (équation de la chaleur) ;
  117. *
  118. 'DEBPROC' CALCUL nraff*'ENTIER' typdis*'MOT' ;
  119. *
  120. * titre global pour les dessins
  121. *
  122. titglob = 'CHAINE' ' ; nraff=' nraff ' typdis=' typdis ;
  123. *
  124. * Paramètres physiques
  125. *
  126. mu = 1. ;
  127. * Conditions aux limites
  128. alpha = (** 2. 0.5) ;
  129. beta  = (** 3. 0.5) ;
  130. gamma = (** 5. 0.5) ;
  131. 'SI' ('EGA' typdis 'QUAF') ;
  132.   eta = (** 7. 0.5)  ;
  133. 'SINON' ;
  134.   eta = 0. ;
  135. 'FINSI' ;    
  136. delta = (** 1.5 0.5) ;
  137. *
  138. *  Géométrie
  139. *
  140. pA = 0.   0.  0. ;
  141. pB = 1.   0.  0. ;
  142. pC = 1.   1.  0. ;
  143. pD = 0.   1.  0. ;
  144. pE = 0.   0.  1. ;
  145. pF = 1.   0.  1. ;
  146. pG = 1.   1.  1. ;
  147. pH = 0.   1.  1. ;
  148. *
  149. * Paramètres de la discrétisation de base
  150. *
  151. 'SI' complet ;
  152.   nAB = 6 ;
  153.   nBC = 6 ;
  154.   nCD = 6 ;
  155.   nDA = 6 ;
  156.   nH  = 6 ;
  157. 'SINON' ;
  158.   nAB = 3 ;
  159.   nBC = 3 ;
  160.   nCD = 3 ;
  161.   nDA = 3 ;
  162.   nH  = 3 ;
  163. 'FINSI' ;  
  164. *
  165. * Géométrie discrétisée
  166. *
  167. 'OPTION' 'ELEM' 'CUB8' ;
  168. bas    = 'DROIT' nAB pA pB ;
  169. droite = 'DROIT' nBC pB pC ;
  170. haut   = 'DROIT' nCD pC pD ;
  171. gauche = 'DROIT' nDA pD pA ;
  172. spourt = bas 'ET' droite 'ET' haut 'ET' gauche ;
  173. sbas   = 'SURFACE' spourt 'PLAN' ;
  174. mt = 'VOLUME' sbas 'TRAN' nH  pE ;
  175. *mt = 'ORIENTER' mt ;
  176. f1 f2 f3 = 'FACE' mt  ;
  177. f1 = 'CHANGER' f1 'TRI3' ;
  178. f2 = 'CHANGER' f2 'TRI3' ;
  179. f3 = 'CHANGER' f3 'TRI3' ;
  180. ss = (f1 'ET' f2 'ET' f3) ;
  181. 'OPTION' 'ELEM' 'TET4' ;
  182. mt = 'VOLUME' ss ;
  183. pourtour = 'ENVELOPPE' mt ;
  184. *
  185. * Rotation et translation aditionnelle + bruit blanc
  186. * + raffinement
  187. *
  188. vtran = ((** 2 0.5) (* -1 (** 3 0.5)) (** 5 0.5)) ;
  189. orig  = (0.D0 0.D0 0.D0) ;
  190. orig2  = (-0.5D0 0.5D0 0.5D0) ;
  191. arot  = (** 360. 0.5) ;
  192. *arot  = 0.D0 ;
  193. *arot = 90.0D0 ;
  194. ladens = (** ('/' ('MESURE' mt) ('NBEL' mt)) ('/' 1.D0 3.D0)) ;
  195. *'MESSAGE' ('CHAINE' 'ladens=' ladens) ;
  196. ampbruit = (0.02 * ladens) ;
  197. *
  198. * On rajoute le bruit
  199. *
  200. pmt = 'CHANGER' mt 'POI1' ;
  201. pcnt= 'CHANGER' pourtour 'POI1' ;
  202. inmt= 'DIFF' pmt pcnt;
  203. brbl1 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
  204. brbl2 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
  205. brbl = 'ET' ('NOMC' 'UX' brbl1) ('NOMC' 'UY' brbl2) ;
  206. mt = 'PLUS' mt brbl;
  207. *
  208. * On raffine nraff fois 
  209. *
  210. 'SI' (nraff > 0) ;
  211.   'REPETER' bcl nraff ;
  212.     mt  = 'CHANGER' mt 'QUADRATIQUE' ;
  213.     $mt = 'DOMA' mt 'MACRO' ;
  214.     mt  = ($mt . 'MAILLAGE') ;
  215.     pourtour  = 'CHANGER' pourtour 'QUADRATIQUE' ;
  216.     $pourtou = 'DOMA' pourtour 'MACRO' ;
  217.     pourtour  = ($pourtou . 'MAILLAGE') ;
  218.     'MENAGE' ;
  219.   'FIN' bcl ;
  220. 'FINSI' ;
  221. *
  222. * Eventuellement, on trace le résultat
  223. *
  224. 'SI' graph ;
  225.    titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
  226.                       ' NBELEM=' ('NBEL' mt) titglob ;
  227.    'TRACER' mt 'NOEUD' 'TITRE' titgeo ;
  228. 'FINSI' ;
  229. *
  230. * Modèle
  231. *
  232. _mt     = 'CHANGER' mt 'QUAF' ;
  233. _pourt  = 'CHANGER' pourtour 'QUAF' ;
  234. 'ELIMINATION' (_mt 'ET' _pourt) 1.D-6 ;
  235. $mt = 'MODELISER' _mt 'NAVIER_STOKES' typdis ;
  236. $pourt = 'MODELISER' _pourt 'NAVIER_STOKES' typdis ;
  237. *
  238. * Solution exacte trilinéaire
  239. *
  240. xx yy zz = 'COORDONNEE' ('DOMA' $mt 'SOMMET') ;
  241. solex = '+' ('+' ('*' xx alpha) ('*' yy beta)) ('*' zz gamma) ;
  242. *solex = '+' solex ('*' ('*' ('*' xx yy) zz) eta) ;
  243. solex = '+' solex ('*' ('*' xx yy) eta) ;
  244. solex = '+' solex delta ;
  245. solex = 'KCHT' $mt 'SCAL' 'SOMMET' solex ;
  246. *
  247. * On tourne et on translate
  248. *
  249. tmt tpourt _tmt _tpourt solex =
  250. 'TOURNER' mt pourtour _pourt _mt solex arot orig orig2 ;
  251. tmt tpourt _tmt _tpourt solex =
  252. 'PLUS' tmt tpourt _tpourt _tmt solex vtran ;
  253. 'ELIMINATION' 1.D-6 (_tmt 'ET' _tpourt 'ET' tmt 'ET' tpourt
  254.                      'ET' ('EXTRAIRE' solex 'MAIL')) ;
  255. $tmt = 'MODELISER' _tmt 'NAVIER_STOKES' typdis ;
  256. $tpourt = 'MODELISER' _tpourt 'NAVIER_STOKES' typdis ;
  257. solex = 'KCHT' $tmt 'SCAL' 'SOMMET' solex ;
  258. mtmt = 'DOMA' $tmt 'MAILLAGE' ;
  259. cmtmt = 'DOMA' $tpourt 'MAILLAGE' ;
  260. *
  261. * Mise en place du calcul numérique
  262. *
  263. * équation de Laplace
  264. *
  265. rv = 'EQEX' $tmt 'ITMA' 1
  266. 'OPTI' 'EF' 'IMPL'
  267. 'ZONE' $tmt 'OPER' 'LAPN' mu 'INCO' 'TN' 'TN' ;
  268. *
  269. * conditions aux limites
  270. *
  271. csolex = 'REDU' solex cmtmt ;
  272. rv = 'EQEX' rv
  273.  'CLIM' 'TN' 'TIMP' cmtmt csolex ;
  274. *                      ________________
  275. *
  276. *                       INITIALISATION
  277. *                      ________________
  278. *
  279. rv . 'INCO' = 'TABLE' 'INCO' ;
  280. rv . 'INCO' . 'TN' = 'KCHT' $tmt 'SCAL' 'SOMMET' 0. ; 
  281. *
  282. rv . 'NITER' = 1 ;
  283. rv . 'METHINV' . 'TYPINV'   = 1 ;
  284. rv . 'METHINV' . 'IMPINV'   = 2 ;
  285. *                         __________
  286. *
  287. *                           CALCUL
  288. *                         __________
  289. *
  290. EXEC rv ;
  291. *
  292. * Résultats
  293. *
  294. 'SI' graph ;
  295. *
  296. * solutions exactes
  297. *
  298.   tn = solex ;
  299.   titt = 'CHAINE' 'Température exacte' titglob ;
  300.   'TRACER' tn tmt tpourt nbisov 'TITRE' titt ;
  301. *
  302. * solutions calculées
  303. *
  304. *  tn = rv . 'INCO' . 'TN' ;
  305. *  titt = 'CHAINE' 'Température calculée' titglob ;
  306. *  'TRACER' tn tmt tpourt nbisov 'TITRE' titt ;
  307. 'FINSI' ;
  308. *
  309. * Calcul des erreurs par rapport à la solution analytique
  310. *
  311. tn = rv . 'INCO' . 'TN' ;
  312. errlit = CALCERR  tn solex ;
  313. errl2t = CALCERR2 tn solex ;
  314. 'MESSAGE' '-------------------------------------------------' ;
  315. 'MESSAGE' ('CHAINE' 'Erreur en norme LINF' titglob ' : ' errlit) ;
  316. 'MESSAGE' ('CHAINE' 'Erreur en norme L2  ' titglob ' : ' errl2t) ;
  317. 'MESSAGE' '-------------------------------------------------' ;
  318. 'FINPROC' errl2t ;
  319. * Fin de la procédure de calcul                             *
  320. *___________________________________________________________*
  321. *
  322. *-----------------------------------------------------------
  323. * Paramètres du calcul
  324. *
  325. * ldiscr : type d'éléments à tester.
  326. 'SI' complet ;
  327.   ldiscr = 'MOTS' 'QUAF' 'LINE' ;
  328. 'SINON' ;  
  329.   ldiscr = 'MOTS' 'QUAF' 'LINE' ;
  330. 'FINSI' ;  
  331. *
  332. *-----------------------------------------------------------*
  333. * Boucles sur les différents paramètres du calcul           *
  334. ok = VRAI ;
  335. *
  336. *  Boucle sur les discrétisations
  337. *
  338. 'REPETER' idiscr ('DIME' ldiscr) ;
  339.   discr = 'EXTRAIRE' ldiscr &idiscr ;
  340.   errl2t = CALCUL 0 discr ;
  341.   ok = ('ET' ok (errl2t '<' 1.D-14)) ;
  342.   'SI' ('ET' ('NON' complet) ('NON' ok)) ;
  343.     'MESSAGE' 'On a dépassé les marges d_erreur :' ;
  344.     'MESSAGE' ' * 1.D-14  pour une méthode directe' ;
  345.     'MESSAGE' ' * 1.D-12  pour une méthode itérative' ;
  346. *    'ERREUR' 5 ;
  347.   'FINSI' ;
  348. 'FIN' idiscr ;
  349.  
  350. 'SI' interact ;
  351.    'OPTION' 'DONN' 5 ;
  352. 'FINSI' ;
  353. 'MESSAGE' 'Tout s_est bien passé' ;
  354. 'FIN' ;
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales