* fichier : cmct1.dgibi
*
* 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';
t0 = temps noec ;
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 ;
t01 = temps noec ;
t01 = t01 + 1;
temps zero ;
'MESS' 'Nombre de segment Active' difsg ;
'MESS' 'Memoire correspondante' difmsg ;
errsg = >eg difsg 2 ;
'SI' errsg ;
'MESS' 'Erreur dans la desactivation des segments' ;
'ERREUR' 5 ;
'FINSI' ;
*============== resolution directe =======================
*
dep3 = 'RESO' rigt1 ( chpd0 'ET' chpfi ) ;
chpf3 = 'REAC' cl3 dep3 ;
dep3 = 'ENLE' dep3 'LX' ;
t1 = temps noec ;
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' ;
© Cast3M 2003 - All rights reserved.
Disclaimer