* fichier : dyne02.dgibi
************************************************************************
* Section : Mecanique Dynamique
************************************************************************
*
************************************************************************
*
*  Impact d'une barre elastique sur un plan rigide :
*  p1             p2   p2b
*   +--------------+   +
*      ->           <->
*      P=M*g        JeuLiai=0.01
*
*  Calcul sur base modale via l'operateur DYNE
* 
*  Creation, modifs :
*  D. Combescure, Aout 2006
*  BP, 2018-12-19 : Test des differents algorithme d'integration 
*                   et de l'option K et M PLEINE
*
************************************************************************
*
GRAPH_boucle = FAUX;
GRAPH = VRAI; opti trac PSC EPTR 5 POTR HELVETICA_16;
*GRAPH = FAUX;
*
opti dime 3 elem seg2;
*
************************************************************************
*
* MAILLAGE, MODELE, MATRICES ET BASE MODALE
*
************************************************************************

L1 = 1.;
n1 = 10;
L =  2.;
n2 = 10;
*
p1 = 0. 0. 0.;
p2 = L1 0. 0.;
p2b = L1 0. 0.;
p3 = L 0. 0.; 
lig1 = d n1 p1 p2;
lig2 = d n2 p2b p3;
mesh = lig1;
*
S1 = 0.1;
*
E1 = 2000.;
nu1 =0.2;
ro1 = 7.8D-3;
*
mod1 = MODE mesh mecanique elastique BARR;
mat1 = MATE mod1 YOUN E1 NU nu1 RHO ro1
                 SECT S1;
* 
rig1 = RIGI mod1 mat1;
mas1 = MASS mod1 mat1;
bl1 = (BLOQ lig1 UY) et (BLOQ lig1 UZ);
bl2 = (BLOQ lig2 UY) et (BLOQ lig2 UZ);
bl2b = (BLOQ p2b DEPL);
*
For1 = MAS1 * (MANU CHPO mesh 1 UX (-9.81));
*
nmode = 2;
TMOD = VIBR 'PROCHE' (PROG nmode) (LECT nmode) 
      (RIG1 et bl1 et bl2b) MAS1 ;
*
SI GRAPH_boucle;
oeil = 0. 0. 100.;
Repeter bou1 nmode;
  FIFI = TMOD . MODES . &Bou1 . 'FREQUENCE'      ;
  TT = TMOD . MODES . &Bou1 .  'DEFORMEE_MODALE'  ;
  DEf1 = DEFO mesh TT;
  titre (CHAIN 'Frequence=' Fifi);
  trac oeil def1;
fin bou1;
FINSI;


************************************************************************
*
* CALCUL NON-LINEAIRE : Appel à DYNE
*
************************************************************************

*************************
* TEMPS
*************************
*
* DT  = 1.E-4 ; NTT = 5000 ; NSORT = 1; 
*--> resultats non converges !
* DT  = 2.E-5 ; NTT = 25000; NSORT = 5; 
*--> convergence sur 1er rebond mais pas sur le 2eme
DT  = 1.E-5 ; NTT = 50000; NSORT = 10;
*--> convergence jusqu'au 3eme rebond : suffisant pour le TNR
* DT  = 1.E-6 ; NTT = 500000; NSORT = 50;
* *--> convergence jusqu'au 5eme rebond
*
*************************
* TABLE AMORTISSEMENT
*************************
*
TAMOR = TABLE 'AMORTISSEMENT'       ;
AMO = PROG NMODE*0.0                ;
TAMOR . 'AMORTISSEMENT' = AMOR TMOD AMO ;
*
***************************************************
* TABLE CHARGEMENT (CHARG1 CHARG2 CHARG3)
***************************************************
TCHAR = TABLE 'CHARGEMENT'          ;
*
Ampl = 1.0;
Evotime = evol manu (prog 0. 0.01 100.)
                    (Ampl*(prog 0. 1. 1.));
CHA1 = CHAR For1 Evotime; 
CHAP1 = PJBA TMOD CHA1;     
TCHAR.'BASE_A' = CHAP1;

***********************************************************
* LIAISONS
***********************************************************
*
TILIA = TABLE 'LIAISON'  ;
TTLB = TABLE 'LIAISON_B' ;
TILIA.'LIAISON_B' = TTLB ;

* LIAISON POINT POINT
Kliai = 30.0*E1*S1/L1;
mess Kliai;
Jeuliai = .01;
TLA12 = TABLE 'LIAISON_ELEMENTAIRE' ;
TTLB .  1 = TLA12 ;
TLA12 . 'TYPE_LIAISON' = MOT 'POINT_POINT';
TLA12 . 'POINT_A' = P2 ;
TLA12 . 'POINT_B' = P2b ;
TLA12. 'NORMALE' = ( -1. 0. 0.);
TLA12 . 'RAIDEUR' = Kliai ;
TLA12 . 'JEU' = JeuLiai      ;
TLA12 . 'AMORTISSEMENT' = 0. ;

*************************
* TABLE SORTIE
*************************
*
TABSOR = TABLE 'SORTIE';
*
TABSOR2 = TABLE 'LIAISON_B' ;
TABSOR.'LIAISON_B' = TABSOR2  ;
*
TVARPP = TABLE 'VARIABLE'    ;
TVARPP.'FORCE_DE_CHOC_POINT_A' = VRAI ;
TVARPP.'FORCE_DE_CHOC_POINT_B' = VRAI ;
TVARPP.'UY_POINT_A' = VRAI ;
TVARPP.'UY_POINT_B' = VRAI ;
TVARPP.'VITESSE' = VRAI ;
*
TVARGD = TABLE 'VARIABLE'    ;
TVARGD.'FORCE_DE_CHOC' = VRAI ;
TVARGD.'VITESSE' = VRAI ;
TVARGD.'UY' = VRAI ;
*
TABSOR2. TLA12 = TVARPP;

liste_coul = MOTS 'MARI' 'BLEU' 'BOUT' 'VERT' 'BRIQ' 'ROUG';
Tdess1 = tabl;
Tdess1 . 'TITRE' = tabl;
*=======================================================================
*======>  BOUCLE SUR DIVERSES SYNTAXES ET ALGO POUR LES TESTER

*on a le droit a des nom de variables jusqu'a 24 caracteres
prog_ERR1 = PROG ;
prog_ERR2 = PROG ;
prog_ERR3 = PROG ;
UXP2_TOT    = VIDE 'EVOLUTIO';
EVFchoA_tot = VIDE 'EVOLUTIO'; 

REPE BTEST 6;

* choix algo
si ((EGA &BTEST 1) ou (EGA &BTEST 2));
  mot_algo = MOT 'DE_VOGELAERE'; mot_algo_court=mot 'DV';
finsi;
si ((EGA &BTEST 3) ou (EGA &BTEST 4));
  mot_algo = MOT 'DIFFERENCES_CENTREES'; mot_algo_court=mot 'CD';
finsi;
si ((EGA &BTEST 5) ou (EGA &BTEST 6));
  mot_algo = MOT 'ACCELERATION_MOYENNE';mot_algo_court=mot 'AA';
finsi;
Tdess1 . 'TITRE' . &BTEST = mot_algo_court;

* choix syntaxe (base modale / table raideur et masse)
BPAIR = MULT &BTEST 2;

************************************************
* CALCUL TEMPOREL PROPREMENT DIT
***********************************************

si BPAIR;
* SYNTAXE BASE_MODALE :
 titre1 = chai 'DYNE ' ' ' mot_algo ' (syntaxe base_modale)'  ;
 TRESU = DYNE mot_algo TMOD TAMOR TCHAR TILIA
                      NTT DT NSORT TABSOR  ;
sinon;
* SYNTAXE TABLE K ET M
 titre1 = chai 'DYNE ' ' ' mot_algo ' (syntaxe table K et M)'  ;
  si (&BTEST <eg 2);
    TKETM = TABL 'RAIDEUR_ET_MASSE';
    TKETM . 'RAIDEUR' = PJBA RIG1 TMOD;
    TKETM . 'MASSE'   = PJBA MAS1 TMOD;
    TKETM . 'NATURE_RAIDEUR' = MOT 'PLEINE';
    TKETM . 'NATURE_MASSE'   = MOT 'PLEINE';
    TKETM . 'BASE_MODALE'    = TMOD;
  finsi;
  TRESU = DYNE mot_algo TKETM TAMOR TCHAR TILIA
                        NTT DT NSORT TABSOR  ;
finsi;


************************************************
* TRACE DES DEPLACEMENTS EN FONCTION DU TEMPS
**********************************************
COUL1 = EXTR liste_coul &BTEST ;

UXP2 = EVOL COUL1 RECO TRESU TMOD 'DEPL' P2  UX  ;
UXP3 = EVOL COUL1 RECO TRESU TMOD 'DEPL' P3  UX  ;
UXP3 = evol COUL1 manu (extr UXP3 absc) 
  ((extr UXP3 ordo) + (PROG (dime (extr UXP3 ordo) )*(-1.*JeuLiai)));
VXP2 = EVOL COUL1 RECO TRESU TMOD 'VITE' P2  UX  ;
VXP3 = EVOL COUL1 RECO TRESU TMOD 'VITE' P3  UX  ;
*
TDS = TRESU . 'TEMPS_DE_SORTIE' ;
RESCHOC = TRESU . TLA12;
FchocB = RESCHOC. FORCE_DE_CHOC_POINT_B;
FchocA = RESCHOC. FORCE_DE_CHOC_POINT_A;
*
EVFchoA = 'EVOL' COUL1 'MANU' TDS FChocA;
EVFchoB = 'EVOL' COUL1 'MANU' TDS FChocB;
*
UXP2_TOT    = UXP2_TOT    ET UXP2;
EVFchoA_tot = EVFchoA_tot ET EVFchoA;
SI GRAPH_boucle;
 ta = table;
 ta.2 = 'TIRR';
 DESS (UXP2 ET UXP3) ta XBORD 0. 0.20 TITRE titre1 ;
 DESS (VXP2 ET VXP3) ta XBORD 0. 0.20 TITRE titre1 ;
*  DESS EVFchoA  XBORD 0. 0.20 ;
*  DESS EVFchoB  XBORD 0. 0.20 ;
FINSI;


***********************************************
* TEST DE NON REGRESSION (TNR)
* par rapport a une precedente execution du cas test
***********************************************
* on teste le deplacement à t=0.1s (fin du 1er rebond)
tt = 0.1;

* ref pour DT=1.E-4 :  uref = -6.71056E-03;
* nouvelle ref. (BP, 2018-12-12)
uref = -6.63459E-03;
uevol = IPOL UXP2 tt;
* test du fonctionnement de RECO 
uphys = RECO TRESU TMOD 'DEPL' tt;
ureco = EXTR uphys P2 'UX';

* on compare le resultat de RECO et de EVOL RECO 
MESS 'valeur calculee : ' uref ureco uevol;
ERR1 = ABS ((ureco - uref) / uref);
ERR2 = ABS ((uevol - uref) / uref);
ERR3 = ABS ((ureco - uevol) / uref);
MESS 'erreur relative : ' ERR1 ERR2 ERR3;
MESS '----------------------------------------------------------------';

*on a le droit a des nom de variables jusqu'a 24 caracteres
prog_ERR1 = prog_ERR1 et ERR1 ;
prog_ERR2 = prog_ERR2 et ERR2 ;
prog_ERR3 = prog_ERR3 et ERR3 ;

FIN BTEST;
*<======  FIN DE BOUCLE SUR DIVERSES SYNTAXES ET ALGO POUR LES TESTER
*=======================================================================

ERR1 = maxi prog_ERR1 ;
ERR2 = maxi prog_ERR2 ;
ERR3 = maxi prog_ERR3 ;

SI GRAPH;
 DESS UXP2_TOT XBORD 0. 0.5 TITRE 'dyne02' LEGE  Tdess1;
 DESS UXP2_TOT XBORD 0. 0.20 TITRE 'dyne02' LEGE 'NE' Tdess1; 
 DESS UXP2_TOT XBORD 0.2 0.3 TITRE 'dyne02' LEGE 'NO' Tdess1;;
 DESS UXP2_TOT XBORD 0.44 0.5 XGRA 0.01 TITRE 'dyne02' LEGE 'NO' Tdess1;;
 DESS EVFchoA_tot  XBORD 0. 0.5 TITRE 'dyne02' LEGE Tdess1;;
 DESS EVFchoA_tot  XBORD 0. 0.20 TITRE 'dyne02';;
 DESS EVFchoA_tot  XBORD 0.2 0.30 TITRE 'dyne02';;
FINSI;

SI ((ERR1 > 0.005) OU  (ERR2 > 0.005) OU  (ERR3 > 0.0001));
  ERRE 5;
FINSI;

* OPTI DONN 5 TRAC 'X';
FIN ;