* fichier :  condense1.dgibi
************************************************************************
* Section : Mathematiques Fonctions
************************************************************************
* 
*     CAST TEST PORTANT SUR L'OPÉRATEUR SUPE
*
*      Principe :
*        On considère une matrice formée de l'assemblage
*           d'une matrice de rigidité rig1
*           de matrice de blocage nuls cl1 et cl2
*           d'une matrie de masse massc1
*           d'une  de blocage cl3
*        et un chargement formé de
*          de deplacement non nuls sur cl3: chpd0
*          de forces reparties : chpfi
*
*       On resoud indirectement par condensation sur les
*          multiplicateurs de Lagrange de cl3 puis par
*          une redescente sur tous les neouds
*
*       On résoud directement
*       
*   
*
*     On teste en particulier la normalisation des
*     multiplicateur de Lagrange lorsque'ils sont maitres
*     en choissant un module d'YOUNG une densité  élévés
*
*     degay 20 02 97
*  


*==========  taille du maillage ======================
m = 30 ;
n =  m ;
E1 = 1.d15 ;
r0 = 7.d6 ;
graph = faux ; 
*==========  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. ) ; 

*============== modle 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 = 'MASS' mod1 mat1 ;
mas2 = 'MASS' mod2 mat2 ;
massc1 = mas1 'ET' mas2 ;

rig1 = 'RIGI' ( mod1 'ET' mod2 ) ( mat1 'ET' mat2 ) ;


*============ matrice lumpée "elementaire" ============

masd12 = lump  mod2 mat2 ;

*================= conditions aux limites =============
cl1 = bloq p1 'UX' ;
cl2 = 'SYMT' lr0 'DEPL' P1  P4 0.01 ; 
*
cl3 = 'SYMT' lr10 'DEPL' P2  P3 0.01 ;
* cl3 = 'BLOQ' lr10 'UX' ;

rigt1 = rig1 et cl1 et cl2 et massc1 ;
rigt2 = rigt1 et cl3  ;

*================= Condensation ========================
mailmult = 'EXTR' cl3 'MAIL' 'MULT' ;
sup1 = 'SUPE' 'RIGI'  rigt2  mailmult ;
rigc1 = 'EXTR' sup1 'RIGI' ;

*================= Chargement ===========================
* on impose un deplacement unité de lr10 selon l'axe de la plaque
chpd0 = 'DEPI' cl3 1. ;
*  une alternative à DEPI est la commande suivante
* chp0 = 'MANU' 'CHPO' ( mailmult ) 1 'FLX' 1.  ; 


*  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   ====================
* 
chpfci = 'SUPE'  sup1 'CHAR' (chpfi et chpd0) ;

depfl0 = 'RESO' rigc1  chpfci  ;

* on redescend

dep0 = 'SUPER' 'DEPL' sup1  depfl0 (chpfi et chpd0);
chpf0 = 'REAC' cl3 dep0;

*============== resolution directe =======================
*
dept1 = 'RESO'   rigt2  ( chpd0 'ET' chpfi ) ;
dep1 = 'ENLE'  dept1  'LX' ; 
chpf1 = 'REAC' cl3 dept1 ;


*=============== Comparaison  ============================
* 

'OPTI' 'ECHO' 0 ;
* Erreur sur le deplacement 
deperr =  dep0 - dep1 ;
erx = 'EXCO' deperr 'UX' 'SCAL' / ( 'MAXI' 'ABS' ( 'EXCO'  dep1 'UX'));
errx = ( 'MAXI' 'ABS' erx) ;
MESS 'Erreur sur X' ( errx  * 100.) '%'  ;
ery = 'EXCO' deperr 'UY' 'SCAL' / ( 'MAXI' 'ABS' ( 'EXCO'  dep1 '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 = chpf1 '-' chpf0 ;
errf = 'MAXI' 'ABS' ferr / ( 'MAXI' 'ABS' chpf0 ) ;
'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' ;