* fichier : cmct1.dgibi ************************************************************************ * Section : Mathematiques Fonctions * Section : Mecanique Elastique ************************************************************************ * * CAST TEST PORTANT SUR L'OPÉRATEUR CMCT et l'opérateur * * * Principe : * On considère une matrice formée de l'assemblage * de matrice de blocage nuls cl1 et cl2 * d'une matrie de masse diagonale massc1 (LUMP) * et un chargement formé * de forces reparties : chpfi * de second membre sur les relations * * On resoud indirectement par condensation sur les * multiplicateurs de Lagrange de cl3 puis par * une redescente sur tous les noeuds. * * On résoud directement * * On compare * * * degay 22 04 97 * *========== taille du maillage ====================== m = 30 ; n = m ; E1 = 1.d15 ; r0 = 7.d35 ; graph = faux ; * coefficient qui sert à abimer le conditionnement du système alpha = 1.d6 ; *========== maillage ================================ 'OPTI' 'DIME' 2 'MODE' 'PLAN' 'DEFO'; 'OPTI' 'ELEM' 'SEG2'; p1 = 0. 0. ; p2 = 10. 0. ; li1 = d p1 n p2 ; 'OPTI' 'ELEM' 'QUA4' ; su1 = 'TRAN' li1 m ( 0. 1. ) ; ls1 = 'COTE' 3 su1 ; P3 = 'POINT' su1 'PROC' ( 10. 1. ) ; P4 = 'POINT' su1 'PROC' ( 0. 1. ) ; lr0 = 'COTE' 4 su1 ; lr10 = 'COTE' 2 su1 ; * rotation du maillage de 45 pour avoir des conditions * aux limites composées 'DEPL' su1 'TOUR' 45. ( 0. 5. ) ; *============== modele et matériau =================== mod1 = 'MODE' su1 mecanique elastique ; mod2 = 'MODE' ( li1 et ls1 ) mecanique elastique coq2 ; mat1 = 'MATE' mod1 'YOUN' E1 'NU' 0.3 'RHO' r0 ; mat2 = 'MATE' mod2 'YOUN' E1 'NU' 0.3 'RHO' r0 'EPAI' 0.05 ; *=============== matrices classique ================== mas1 = 'LUMP' mod1 mat1 ; mas2 = 'LUMP' mod2 mat2 ; massc1 = mas1 'ET' mas2 ; * stockage sous la forme d'un champ point chpu1 = 'MANU' 'CHPO' su1 4 'UX' 1. 'UY' 1. 'UZ' 1. 'RZ' 1. ; chpm1 = massc1 '*' chpu1 ; chpm1 = 'NOMC' chpm1 ( 'MOTS' 'FX' 'FY' 'MZ' ) ( 'MOTS' 'UX' 'UY' 'RZ' ) 'NATU' 'DISCRET' ; * inversion du champ point chpm1 = 'INVE' chpm1 ; *================= conditions aux limites ============= cl0 = 'SYMT' lr10 'DEPL' P2 P3 0.01 ; * on multiplie la matrice par alpha pour tester le conditionnement cl1 = ( 'BLOQ' p1 'RZ' ) * alpha ; cl2 = 'SYMT' lr0 'DEPL' P1 P4 0.01 ; cl3 = cl0 'ET' cl1 'ET' cl2 ; * matrice globale du systeme rigt1 = cl3 'ET' massc1 ; * maillage support des multiplicateurs mailmult = 'EXTR' cl3 'MAIL' 'MULT' ; *================= Chargement =========================== * on impose un deplacement unité de lr10 selon l'axe de la plaque chpd0 = ( 'DEPI' cl0 1.d0 ) 'ET' ( 'DEPI' cl1 ALPHA ) ; * on applique des forces sur le maillage chpfi = 'CHPOINT' 'ALEATOIRE' ( mas1 et mas2 ) ; chpfi = 'NOMC' chpfi ( 'MOTS' 'UX ' 'UY ' 'RZ ' ) ( 'MOTS' 'FX ' 'FY ' 'MZ ' ) 'NATU' 'DISCRET' ; *============ résolution indirecte ==================== temps zero ; *-- avec SUPER sup1 = 'SUPE' 'RIGI' rigt1 mailmult ; rigc1 = 'EXTR' sup1 'RIGI' ; **chpfci1 = 'SUPE' sup1 'CHAR' chpfi ; **chpd01 = 'SUPE' sup1 'CHAR' chpd0; **depfl1 = 'RESO' rigc1 ( chpd01 'ET' chpfci1 ) ; * on redescend **chpf1 = 'REAC' cl3 depfl1 ; **dep1 = 'RESO' massc1 ( chpf1 'ET' chpfi ); chpfci1 = 'SUPE' sup1 'CHAR' (chpd0 'ET' chpfi); depfl1 = 'RESO' rigc1 chpfci1 ; dep1 = super depl sup1 depfl1 (chpd0 'ET' chpfi); chpf1 = 'REAC' cl3 dep1 ; dep1 = 'ENLE' dep1 'LX'; TABTPS = TEMP 'NOEC'; t0 = TABTPS.'TEMPS_CPU'.'INITIAL' ; temps zero ; *-- avec CMCT rigc2 = 'CMCT' cl3 chpm1 ; mac0 sgac0 = temps sgac ; depsl2 = '*' chpm1 (MOTS 'UX' 'UY' 'RZ') chpfi (MOTS 'FX' 'FY' 'MZ') (MOTS 'UX' 'UY' 'RZ') ; mac1 sgac1 = temps sgac ; difsg = sgac1 '-' sgac0 ; difmsg = mac1 '-' mac0 ; * second membre du systeme sur les multiplicateurs fcl2 = cl3 '*' depsl2 '-' chpd0 ; depfl2 = 'RESO' rigc2 fcl2 ; chpf2 = ( 'REAC' cl3 depfl2 ) ; depcor2 = '*' chpm1 (MOTS 'UX' 'UY' 'RZ') chpf2 (MOTS 'FX' 'FY' 'MZ') (MOTS 'UX' 'UY' 'RZ') 'NATURE' 'DIFFUS' ; dep2 = depsl2 + depcor2 ; TABTPS = TEMP 'NOEC'; t01 = TABTPS.'TEMPS_CPU'.'INITIAL' ; t01 = t01 + 1; temps zero ; 'MESS' 'Nombre de segment Active' difsg ; 'MESS' 'Memoire correspondante' difmsg ; errsg = >eg difsg 17; * 2019/12/11 SG Le nouveau paradigme (fiche anomalie 10273) étant de ne * pas désactiver les segments, on commente ce test *new-paradigm 'SI' errsg ; *new-paradigm 'MESS' 'Erreur dans la desactivation des segments' ; *new-paradigm 'ERREUR' 5 ; *new-paradigm 'FINSI' ; *============== resolution directe ======================= * dep3 = 'RESO' rigt1 ( chpd0 'ET' chpfi ) ; chpf3 = 'REAC' cl3 dep3 ; dep3 = 'ENLE' dep3 'LX' ; TABTPS = TEMP 'NOEC'; t1 = TABTPS.'TEMPS_CPU'.'INITIAL' ; t1 = t1 + 1; 'OPTI' 'ECHO' 0 ; 'SAUTER' 1 'LIGNE' ; 'MESS' 'Par rapport à la méthode directe : ' ; 'MESS' 'SUPER' ( 100. * t0 / t1 - 100.) '% de temps' ; 'MESS' 'CMCT ' ( 100. * t01 / t1 - 100.)'% de temps - facteur' (1. * t1 / t01 ) ; 'SAUTER' 1 'LIGNE' ; 'OPTI' 'ECHO' 1 ; *=============== Comparaison ============================ * 'OPTI' 'ECHO' 0 ; * Erreur sur le deplacement deperr = dep3 * 2. - dep2 - dep1 / 2. ; erx = 'EXCO' deperr 'UX' 'SCAL' / ( 'MAXI' 'ABS' ( 'EXCO' dep3 'UX')); errx = ( 'MAXI' 'ABS' erx) ; MESS 'Erreur sur X' ( errx * 100.) '%' ; ery = 'EXCO' deperr 'UY' 'SCAL' / ( 'MAXI' 'ABS' ( 'EXCO' dep3 'UY')); erry = ( 'MAXI' 'ABS' ery) ; MESS 'Erreur sur Y' ( erry * 100.) '%' ; erz = 'EXCO' deperr 'RZ' 'SCAL' ; errz = ( 'MAXI' 'ABS' erz) ; MESS 'Erreur sur RZ' (errz * 1.) 'difference absolue' ; errtot = ( erx ** 2 ) + ( ery ** 2 ) ; 'OPTI' 'ECHO' 1 ; 'SI' graph ; 'TITR' 'Niveau d erreur' ; 'TRAC' errtot su1 ; 'FINSI' ; * erreur sur les forces ferr = chpf3 * 2. '-' chpf2 - chpf1 ; errf = 'MAXI' 'ABS' ferr / ( 'MAXI' 'ABS' chpf3 ) ; 'MESS' 'Erreur sur les forces' ( errf '*' 100.) '%' ; *===================================== * Code de fonctionnement 'SI' ((errx > 1.D-9) 'OU' (erry > 1.D-9) 'OU' (errf > 1.D-9 )); 'ERREUR' 5 ; 'SINON' 'ERREUR' 0 ; 'FINSI' ; 'FIN' ;