Télécharger cmct1.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : cmct1.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. *
  5. * CAST TEST PORTANT SUR L'OPÉRATEUR CMCT et l'opérateur *
  6. *
  7. * Principe :
  8. * On considère une matrice formée de l'assemblage
  9. * de matrice de blocage nuls cl1 et cl2
  10. * d'une matrie de masse diagonale massc1 (LUMP)
  11. * et un chargement formé
  12. * de forces reparties : chpfi
  13. * de second membre sur les relations
  14. *
  15. * On resoud indirectement par condensation sur les
  16. * multiplicateurs de Lagrange de cl3 puis par
  17. * une redescente sur tous les noeuds.
  18. *
  19. * On résoud directement
  20. *
  21. * On compare
  22. *
  23. *
  24. * degay 22 04 97
  25. *
  26.  
  27.  
  28. *========== taille du maillage ======================
  29. m = 30 ;
  30. n = m ;
  31. E1 = 1.d15 ;
  32. r0 = 7.d35 ;
  33. graph = faux ;
  34. * coefficient qui sert à abimer le conditionnement du système
  35. alpha = 1.d6 ;
  36. *========== maillage ================================
  37. 'OPTI' 'DIME' 2 'MODE' 'PLAN' 'DEFO';
  38. 'OPTI' 'ELEM' 'SEG2';
  39. p1 = 0. 0. ;
  40. p2 = 10. 0. ;
  41. li1 = d p1 n p2 ;
  42. 'OPTI' 'ELEM' 'QUA4' ;
  43. su1 = 'TRAN' li1 m ( 0. 1. ) ;
  44. ls1 = 'COTE' 3 su1 ;
  45. P3 = 'POINT' su1 'PROC' ( 10. 1. ) ;
  46. P4 = 'POINT' su1 'PROC' ( 0. 1. ) ;
  47. lr0 = 'COTE' 4 su1 ;
  48. lr10 = 'COTE' 2 su1 ;
  49.  
  50. * rotation du maillage de 45 pour avoir des conditions
  51. * aux limites composées
  52. 'DEPL' su1 'TOUR' 45. ( 0. 5. ) ;
  53.  
  54. *============== modele et matériau ===================
  55. mod1 = 'MODE' su1 mecanique elastique ;
  56. mod2 = 'MODE' ( li1 et ls1 ) mecanique elastique coq2 ;
  57. mat1 = 'MATE' mod1 'YOUN' E1 'NU' 0.3 'RHO' r0 ;
  58. mat2 = 'MATE' mod2 'YOUN' E1 'NU' 0.3 'RHO' r0 'EPAI' 0.05 ;
  59.  
  60. *=============== matrices classique ==================
  61. mas1 = 'LUMP' mod1 mat1 ;
  62. mas2 = 'LUMP' mod2 mat2 ;
  63. massc1 = mas1 'ET' mas2 ;
  64.  
  65. * stockage sous la forme d'un champ point
  66. chpu1 = 'MANU' 'CHPO' su1 4 'UX' 1. 'UY' 1. 'UZ' 1.
  67. 'RZ' 1. ;
  68. chpm1 = massc1 '*' chpu1 ;
  69. chpm1 = 'NOMC' chpm1
  70. ( 'MOTS' 'FX' 'FY' 'MZ' )
  71. ( 'MOTS' 'UX' 'UY' 'RZ' )
  72. 'NATU' 'DISCRET' ;
  73. * inversion du champ point
  74. chpm1 = 'INVE' chpm1 ;
  75.  
  76. *================= conditions aux limites =============
  77. cl0 = 'SYMT' lr10 'DEPL' P2 P3 0.01 ;
  78. * on multiplie la matrice par alpha pour tester le conditionnement
  79. cl1 = ( 'BLOQ' p1 'RZ' ) * alpha ;
  80. cl2 = 'SYMT' lr0 'DEPL' P1 P4 0.01 ;
  81. cl3 = cl0 'ET' cl1 'ET' cl2 ;
  82.  
  83. * matrice globale du systeme
  84. rigt1 = cl3 'ET' massc1 ;
  85.  
  86. * maillage support des multiplicateurs
  87. mailmult = 'EXTR' cl3 'MAIL' 'MULT' ;
  88.  
  89.  
  90. *================= Chargement ===========================
  91. * on impose un deplacement unité de lr10 selon l'axe de la plaque
  92. chpd0 = ( 'DEPI' cl0 1.d0 ) 'ET' ( 'DEPI' cl1 ALPHA ) ;
  93. * on applique des forces sur le maillage
  94. chpfi = 'CHPOINT' 'ALEATOIRE' ( mas1 et mas2 ) ;
  95. chpfi = 'NOMC' chpfi ( 'MOTS' 'UX ' 'UY ' 'RZ ' )
  96. ( 'MOTS' 'FX ' 'FY ' 'MZ ' )
  97. 'NATU' 'DISCRET' ;
  98.  
  99.  
  100. *============ résolution indirecte ====================
  101. temps zero ;
  102.  
  103. *-- avec SUPER
  104. sup1 = 'SUPE' 'RIGI' rigt1 mailmult ;
  105. rigc1 = 'EXTR' sup1 'RIGI' ;
  106. **chpfci1 = 'SUPE' sup1 'CHAR' chpfi ;
  107. **chpd01 = 'SUPE' sup1 'CHAR' chpd0;
  108. **depfl1 = 'RESO' rigc1 ( chpd01 'ET' chpfci1 ) ;
  109. * on redescend
  110. **chpf1 = 'REAC' cl3 depfl1 ;
  111. **dep1 = 'RESO' massc1 ( chpf1 'ET' chpfi );
  112. chpfci1 = 'SUPE' sup1 'CHAR' (chpd0 'ET' chpfi);
  113. depfl1 = 'RESO' rigc1 chpfci1 ;
  114. dep1 = super depl sup1 depfl1 (chpd0 'ET' chpfi);
  115. chpf1 = 'REAC' cl3 dep1 ;
  116. dep1 = 'ENLE' dep1 'LX';
  117.  
  118. TABTPS = TEMP 'NOEC';
  119. t0 = TABTPS.'TEMPS_CPU'.'INITIAL' ;
  120. temps zero ;
  121.  
  122. *-- avec CMCT
  123.  
  124. rigc2 = 'CMCT' cl3 chpm1 ;
  125. mac0 sgac0 = temps sgac ;
  126. depsl2 = '*'
  127. chpm1 (MOTS 'UX' 'UY' 'RZ')
  128. chpfi (MOTS 'FX' 'FY' 'MZ') (MOTS 'UX' 'UY' 'RZ') ;
  129. mac1 sgac1 = temps sgac ;
  130. difsg = sgac1 '-' sgac0 ;
  131. difmsg = mac1 '-' mac0 ;
  132.  
  133. * second membre du systeme sur les multiplicateurs
  134. fcl2 = cl3 '*' depsl2 '-' chpd0 ;
  135. depfl2 = 'RESO' rigc2 fcl2 ;
  136.  
  137. chpf2 = ( 'REAC' cl3 depfl2 ) ;
  138.  
  139. depcor2 = '*'
  140. chpm1 (MOTS 'UX' 'UY' 'RZ')
  141. chpf2 (MOTS 'FX' 'FY' 'MZ')
  142. (MOTS 'UX' 'UY' 'RZ')
  143. 'NATURE' 'DIFFUS' ;
  144.  
  145. dep2 = depsl2 + depcor2 ;
  146.  
  147.  
  148. TABTPS = TEMP 'NOEC';
  149. t01 = TABTPS.'TEMPS_CPU'.'INITIAL' ;
  150. t01 = t01 + 1;
  151. temps zero ;
  152. 'MESS' 'Nombre de segment Active' difsg ;
  153. 'MESS' 'Memoire correspondante' difmsg ;
  154. errsg = >eg difsg 17;
  155. 'SI' errsg ;
  156. 'MESS' 'Erreur dans la desactivation des segments' ;
  157. 'ERREUR' 5 ;
  158. 'FINSI' ;
  159.  
  160. *============== resolution directe =======================
  161. *
  162. dep3 = 'RESO' rigt1 ( chpd0 'ET' chpfi ) ;
  163. chpf3 = 'REAC' cl3 dep3 ;
  164. dep3 = 'ENLE' dep3 'LX' ;
  165.  
  166. TABTPS = TEMP 'NOEC';
  167. t1 = TABTPS.'TEMPS_CPU'.'INITIAL' ;
  168. t1 = t1 + 1;
  169.  
  170. 'OPTI' 'ECHO' 0 ;
  171. 'SAUTER' 1 'LIGNE' ;
  172. 'MESS' 'Par rapport à la méthode directe : ' ;
  173. 'MESS' 'SUPER' ( 100. * t0 / t1 - 100.) '% de temps' ;
  174. 'MESS' 'CMCT ' ( 100. * t01 / t1 - 100.)'% de temps - facteur'
  175. (1. * t1 / t01 ) ;
  176. 'SAUTER' 1 'LIGNE' ;
  177. 'OPTI' 'ECHO' 1 ;
  178.  
  179.  
  180.  
  181. *=============== Comparaison ============================
  182. *
  183.  
  184. 'OPTI' 'ECHO' 0 ;
  185. * Erreur sur le deplacement
  186. deperr = dep3 * 2. - dep2 - dep1 / 2. ;
  187. erx = 'EXCO' deperr 'UX' 'SCAL' / ( 'MAXI' 'ABS' ( 'EXCO' dep3 'UX'));
  188. errx = ( 'MAXI' 'ABS' erx) ;
  189. MESS 'Erreur sur X' ( errx * 100.) '%' ;
  190. ery = 'EXCO' deperr 'UY' 'SCAL' / ( 'MAXI' 'ABS' ( 'EXCO' dep3 'UY'));
  191. erry = ( 'MAXI' 'ABS' ery) ;
  192. MESS 'Erreur sur Y' ( erry * 100.) '%' ;
  193. erz = 'EXCO' deperr 'RZ' 'SCAL' ;
  194. errz = ( 'MAXI' 'ABS' erz) ;
  195. MESS 'Erreur sur RZ' (errz * 1.) 'difference absolue' ;
  196. errtot = ( erx ** 2 ) + ( ery ** 2 ) ;
  197. 'OPTI' 'ECHO' 1 ;
  198.  
  199. 'SI' graph ;
  200. 'TITR' 'Niveau d erreur' ;
  201. 'TRAC' errtot su1 ;
  202. 'FINSI' ;
  203.  
  204. * erreur sur les forces
  205. ferr = chpf3 * 2. '-' chpf2 - chpf1 ;
  206. errf = 'MAXI' 'ABS' ferr / ( 'MAXI' 'ABS' chpf3 ) ;
  207. 'MESS' 'Erreur sur les forces' ( errf '*' 100.) '%' ;
  208.  
  209.  
  210. *=====================================
  211. * Code de fonctionnement
  212.  
  213. 'SI' ((errx > 1.D-9) 'OU' (erry > 1.D-9) 'OU' (errf > 1.D-9 ));
  214. 'ERREUR' 5 ;
  215. 'SINON'
  216. 'ERREUR' 0 ;
  217. 'FINSI' ;
  218.  
  219. 'FIN' ;
  220.  
  221.  
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236.  
  237.  
  238.  

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