* fichier : mazars_triaxial.dgibi
************************************************************************
*      Cas test elementaire pour la loi d'endommagement de Mazars      *
*           Essai de compression sous confinement (triaxial)           *
*                                                                      *
* Les donnees experimentales servant de reference sont issues de :     *
* "Study of the Behavior of Concrete under Triaxial Compression"       *
* Domingo Sfer, Ignacio Carol, Ravindra Gettu, Guillermo Etse          *
* Journal of Engineering Mechanics, 2002, 128(2): 156-163              *
* doi : 10.1061/(ASCE)0733-9399(2002)128:2(156)                        *
*                                                                      *
* Le test est realise avec plusieurs modes de calcul :                 *
*   1) 3D volumique                                                    *
*   2) 2D axisymetrique                                                *
*                                                                      *
* Selon deux modes d'integration du comportement                       *
* realises successivement dans la boucle BCALC :                       *
*   passage 1 : Cast3M                                                 *
*   passage 2 : MFront (pas actif a date)                              *
*                                                                      *
* ITRAC1 : VRAI => affichage resultats                                 *
* IPSC1  : VRAI => affichage resultats dans un PostScript Couleur      *
************************************************************************

DIVERS = 'VENV' 'CASTEM_DIVERS' ;
FIC1   = 'CHAI' DIVERS '/mazars_triaxial_0MPa_courbe0.csv';
FIC2   = 'CHAI' DIVERS '/mazars_triaxial_1.5MPa_courbe0.csv';
FIC3   = 'CHAI' DIVERS '/mazars_triaxial_4.5MPa_courbe0.csv';
FIC4   = 'CHAI' DIVERS '/mazars_triaxial_9MPa_courbe0.csv';

** Options
** - pour les traces
ITRAC1 = VRAI ;
IPSC1  = VRAI ;
SI IPSC1 ;
  OPTI 'TRAC' 'PSC' ;
FINSI ;
** - pour la sauvegarde des résultats de calcul Cast3M et/ou MFront
ISAUV = FAUX ;
SI ISAUV ;
  str1 = 'mazars_triaxial_' ;
FINSI ;

** Liste des modes de calcul a traiter
lcas = LECT 1 2 ;
*lcas = LECT 1 ; COMM '3D volumique' ;
*lcas = LECT 2 ; COMM '2D axisymetrique' ;

** Parametres geometriques
rayo = 75.E-3 ;
long = rayo * (pi ** 0.5) ;
haut = 300.E-3 ;

** Parametres materiau
you   = 27.3E9 ;
nu    = 0.2 ;
epd0  = 1.E-5 ;
at    = 1. ;
ac    = 1.1 ;
bt    = 10000. ;
bc    = 1200. ;
beta  = 1. ;
xrho  = 0. ;
xalph = 0. ;
xtalp = 0. ;
xtref = 0. ;	   

** Parametres chargement
lpres = PROG 0. 1.5E6 4.5E6 9.E6 ;
umax  = 0.015 * haut ;

OPTI ECHO 0 ;
SAUT 1 lign ;
mess 'Parametres' ;
mess '==========';
SAUT 1 lign ;
mess 'Parametre de taille de l element' ;
mess '- long = côté de la section : long (m)' ' ' long ;
mess '- rayo = rayon de la section : rayo (m)' ' ' rayo ;
mess '- haut = hauteur de l element : haut (m)' ' ' haut ;
SAUT 1 lign ;
mess 'Coef. elastiques' ;
mess '- Module d Young         : you (Pa)' you ;
mess '- Coefficient de Poisson : nu' ' ' ' ' ' ' ' ' ' ' ' ' nu ;
SAUT 1 lign ;
mess 'Coef. Mazars' ;
mess '- KTR0 : epd0' ' ' epd0 ;
mess '- ATRA : at' ' ' ' ' ' ' at ;
mess '- BTRA : bt' ' ' ' ' ' ' bt ;
mess '- ACOM : ac' ' ' ' ' ' ' ac ;
mess '- BCOM : bc' ' ' ' ' ' ' bc ;
mess '- BETA : beta' ' ' beta ;
SAUT 1 lign ;
mess 'Chargement' ;
mess '- Pression hydrostatique (confinement)';
REPE phydr (DIME lpres) ;
mess '  -' ' ' &phydr ' : P (Pa)' (EXTR &phydr lpres) ;
FIN phydr ;
mess '- Depl.imp.max : -1*umax (m)' ' ' (-1.*umax) ;
SAUT 1 lign ;
OPTI ECHO 1 ;
*opti donn 5;

** Legende courbes
tpre = TABL ;
tpre . 1 = MOT '0MPa' ;
tpre . 2 = MOT '1.5MPa' ;
tpre . 3 = MOT '4.5MPa' ;
tpre . 4 = MOT '9MPa' ;
tleg = TABL ;
tleg.'TITRE' = TABL ;

*nh145313 .csv recuperes => resext = vrai
resexp = vrai ;

** Oeil en 3D_V (trace a conserver pour conserver l'oeil)
  OPTI 'DIME' 3 'MODE' 'TRID' 'ELEM' 'CUB8' ;
  rep1 = @REPERE (-0.55 -0.55 0.) 'TURQ' ;
  TRAC rep1 'TITR' 'Repere 3D';
  o1 = VALE OEIL ;
  xo1 yo1 zo1 = COOR o1 ;


REPE BCALC 2 ;
*============= A ENLEVER pour realiser le calcul MFront =============
SI (NEG &BCALC 1) ;
  QUIT BCALC;
FINSI ;
*========================== fin A ENLEVER ===========================
** Distinction des résultats de calcul Cast3M ou MFront (PostScript, sauvegarde)
SI (EGA &BCALC 1) ;
  str3 = 'Cast3M' ;
  OPTI ECHO 0 ;
  SAUT 1 LIGN ;
  MESS ' ===================  Comportement materiau : Cast3M  ====================';
  SAUT 1 LIGN ;
SINO ;
  str3 = 'MFront' ;
  OPTI ECHO 0 ;
  SAUT 1 LIGN ;
  MESS ' ===================  Comportement materiau : MFront  ====================';
  SAUT 1 LIGN ;
FINSI ;
OPTI ECHO 1 ;
str4 = CHAI str3 ' :' ;


**************** C A S  #1  :  3 D   V O L U M I Q U E ****************
SI (EXIS lcas 1) ;

** Options generales
  OPTI 'DIME' 3 'MODE' 'TRID' 'ELEM' 'CUB8' ;
  o1 = xo1 yo1 zo1 ;

** Maillage
  p1      = 0. 0. 0. ;
  p2      = 0. long 0. ;
  l12     = DROI 1 p1 p2 ;
  sgauche = l12 TRAN 1 (0. 0. haut) ;
  mail    = sgauche VOLU 'TRAN' 1 (long 0. 0.) ;
  sdroite = mail FACE 2 ;
  env1    = ENVE mail ;
  sbas    = env1 ELEM 5 ;
  shaut   = env1 ELEM 2 ;
  sdevant = env1 ELEM 4 ;
  sderrie = env1 ELEM 3 ;
  p3      = mail POIN 'PROC' (0. 0. haut) ;
  mes1    = MESU mail ;
cadr1 = COUL ((0. (long + 5.E-3) (haut + 5.E-3)) ET (long (long + 5.E-3) (haut + 5.E-3))) 'BLAN' ;
rep1 = @REPERE (-0.16 -0.16 0.) (PROG 0.3 0.3 0.3) 'TURQ' ;
SI ITRAC1 ;
  TRAC o1 'QUAL' (cadr1 ET rep1 ET mail ET (sgauche COUL 'VERT') ET (sdroite COUL 'ROUG')) 'TITR' (CHAI str4 ' '  '[3D V] Maillage - sgauche, sdroite') ;
  TRAC o1 'QUAL' (cadr1 ET rep1 ET mail ET (sbas COUL 'VERT') ET (shaut COUL 'ROUG')) 'TITR' (CHAI str4 ' '  '[3D V] Maillage - sbas, shaut') ;
  TRAC o1 'QUAL' (cadr1 ET rep1 ET mail ET (sdevant COUL 'VERT') ET (sderrie COUL 'ROUG')) 'TITR' (CHAI str4 ' '  '[3D V] Maillage - sdevant, sderrie') ;
  TRAC o1 'QUAL' (cadr1 ET rep1 ET mail ET ((sgauche ET sbas ET sdevant) COUL 'VERT') ET ((sdroite ET shaut ET sderrie) COUL 'ROUG')) 'TITR' (CHAI str4 ' '  '[3D V] Maillage') ;
FINSI ;

** Modele et caracteristiques materiau (&BCALC : 1->Cast3M, 2->MFront)
SI (EGA &BCALC 1) ;
  mo = MODE mail 'MECANIQUE' 'ELASTIQUE' 'ENDOMMAGEMENT' 'MAZARS' ;
  ma = MATE mo 'YOUN' you 'NU' nu 'KTR0' epd0 'ACOM' ac 'BCOM' bc 'ATRA' at 'BTRA' bt 'BETA' beta ;
SINO ;
  coel = MOTS 'YOUN' 'NU' 'RHO' 'ALPH' 'TALP' 'TREF' 'ACOM' 'BCOM' 'ATRA' 'BTRA' 'BETA' 'EPSD';
  statev = MOTS 'DCOM' 'DTRA' 'ALFT' 'ALFC' 'D';
  mo = MODE mail 'MECANIQUE' 'ELASTIQUE' 'NON_LINEAIRE' 'UTILISATEUR' 'LIB_LOI' 'libUmatBehaviour.so' 'FCT_LOI' 'umatmazars' 'C_MATERIAU' coel 'C_VARINTER' statev;
  ma = MATE mo 'YOUN' you 'NU' nu 'RHO' xrho 'ALPH' xalph 'TALP' xtalp 'TREF' xtref 'ACOM' ac 'BCOM' bc 'ATRA' at 'BTRA' bt 'BETA' beta 'EPSD' epd0;
FINSI ;

** Blocages
  blhaut = BLOQ 'UZ' shaut ;
  blrig  = (BLOQ 'UX' sgauche) ET (BLOQ 'UY' sdevant) ET (BLOQ 'UZ' sbas) ;
  bl     = blhaut ET blrig ;

** Chargements (plusieurs pressions imposees, deplacement impose)
  evc1   = EVOL 'MANU' (PROG 0. 0.5 1.) (PROG 0. 1. 1.) ;
  lcoul  = MOTS 'BLEU' 'ORAN' 'VERT' 'ROUG' ;
  evtot0 = VIDE 'EVOLUTIO' ;
  evtot1 = VIDE 'EVOLUTIO' ;
  evtot2 = VIDE 'EVOLUTIO' ;
  npres = DIME lpres ;
  REPE b0 npres ;
* On fait un premier calcul avec seulement la pression imposee
    pre1 = EXTR lpres &b0 ;
*    fpr  = PRES 'MASS' mo (sdroite ET sderrie ET shaut) pre1 ;
    fprx = PRES 'MASS' mo sdroite pre1 ;
    fpry = PRES 'MASS' mo sderrie pre1 ;
    fprz = PRES 'MASS' mo shaut pre1 ;
    fpr = fprx ET fpry ET fprz ;
    chap = CHAR 'MECA' fpr evc1 ;
* Resolution 1 : mise en pression seulement
    t                               = TABL ;
    t . 'MODELE'                    = mo ;
    t . 'CARACTERISTIQUES'          = ma ;
    t . 'BLOCAGES_MECANIQUES'       = blrig ;
    t . 'CHARGEMENT'                = chap ;
    t . 'TEMPS_CALCULES'            = PROG 0. 'PAS' 0.1 0.5 ;
    t . 'MOVA'                      = MOT 'D' ;
    t . 'MES_SAUVEGARDES'           = TABL ;
    t . 'MES_SAUVEGARDES' . 'DEFTO' = VRAI ;
    PASAPAS t ;
* Resolution 2 : pression + ajout du deplacement impose
    dihaut = DEPI blhaut -1. ;
    uz0    = ABS (EXTR (t . 'DEPLACEMENTS' . 5) 'UZ' p3) ;
    evc2   = EVOL 'MANU' (PROG 0.5 1.) (PROG uz0 (umax - uz0)) ;
    chadi  = CHAR 'DIMP' dihaut evc2 ;
** Trace des blocages et du chargement
    re0 = @REPERE (-0.16 -0.16 0.) (PROG 0.3 0.3 0.3) ;
    mor = MODE re0 'MECANIQUE' 'BARR' ;
    mar = MATE mor 'YOUN' 1. 'NU' 1. 'SECT' 1. ;
    rir = RIGI mor mar ;
    vprx = VECT fprx 'FORC' 'VERT' 1.2E-6 ;
    vpry = VECT fpry 'FORC' 'VERT' 1.2E-6 ;
    vprz = VECT fprz 'FORC' 'VERT' 1.2E-6 ;
    vdiz = VECT dihaut ' ' ' ' 'FLX' 'ROUG' 3.E-2 ;
    tit0 = CHAI str4 ' '  '[3D V] Blocages et chargements Uimp et P=' tpre . &b0 ;
    SI ITRAC1 ;
      TRAC o1 ((RIGI mo ma) ET bl ET rir) (vprx ET vpry ET vprz ET vdiz) 'NOLE' 'TITR' tit0 ;
    FINSI ;
    t . 'BLOCAGES_MECANIQUES'       = bl ;
    t . 'CHARGEMENT'                = chap ET chadi ;
    t . 'TEMPS_CALCULES'            = PROG 0.505 'PAS' 0.005 1. ;
    PASAPAS t ;
*    SI ITRAC1 ;
*      tit1 = CHAI '[3D V] UZ vs Temps, P=' tpre . &b0 ;
*      DESS (EVOL 'TEMPE' t 'DEPLACEMENTS' 'UZ' p3) 'TITR' tit1 ;
*      tit1 = CHAI '[3D V] SMZZ vs Temps, P=' tpre . &b0 ;
*      DESS (EVOL 'TEMP' t 'CONTRAINTES' 'SMZZ' 1 1 1) 'TITR' tit1 ;
*    FINSI ;

** Post traitement
* courbes : endommagement moyen vs temps
*           contrainte moyenne vs deformation moyenne
*           force de reaction vs deplacement impose
    tt = t . 'TEMPS' ;
    tc = t . 'CONTRAINTES' ;
    te = t . 'DEFORMATIONS' ;
    tv = t . 'VARIABLES_INTERNES' ;
    tu = t . 'DEPLACEMENTS' ;
    tr = t . 'REACTIONS' ;
    def0 = DEFO (cadr1 ET rep1 ET mail) (tu . 0) 0. ;
    ltps = PROG ;
    lend = PROG ;
    lsig = PROG ;
    leps = PROG ;
    ldep = PROG ;
    lrea = PROG ;
    REPE b1 (DIME tt) ;
      tps1 = tt . (&b1 - 1) ;
      end1 = (INTG mo (tv . (&b1 - 1)) 'D')    / mes1 ;
      sig1 = (INTG mo (tc . (&b1 - 1)) 'SMZZ') / mes1 ;
      eps1 = (INTG mo (te . (&b1 - 1)) 'EPZZ') / mes1 ;
      u1   = EXTR (tu . (&b1 - 1)) 'UZ' p3 ;
      rea1 = 0. ;
      SI ((NBEL (EXTR (tr . (&b1 - 1)) 'MAIL')) > 0) ;
        rea1 = @TOTAL (tr . (&b1 - 1)) sbas 'FZ' ;
      FINSI ;
      ltps = ltps ET tps1 ;
      lend = lend ET end1 ;
      lsig = lsig ET sig1 ;
      leps = leps ET eps1 ;
      ldep = ldep ET u1 ;
      lrea = lrea ET rea1 ;
    FIN b1 ;
    evtot0 = evtot0 ET (EVOL (EXTR lcoul &b0) 'MANU' 'Temps' ltps 'D'    lend) ;
    evtot1 = evtot1 ET (EVOL (EXTR lcoul &b0) 'MANU' 'EPZZ'  leps 'SMZZ' lsig) ;
    evtot2 = evtot2 ET (EVOL (EXTR lcoul &b0) 'MANU' 'UZ'    ldep 'FZ'   lrea) ;
    tleg . 'TITRE' . &b0 = CHAI 'Cal_mat_' str3 '_P=' tpre . &b0 ;
*   deformee
    REPE b2 (DIME t . 'TEMPS') ;
      tps1 = t . 'TEMPS' . (&b2 - 1) ;
*     trace de la deformee (pour controle visuel)
      def1 = DEFO mail (tu . (&b2 - 1)) 10. 'ROUG' ;
      tit1 = CHAI str4 ' '  '[3D V] Deformee au temps' ' ' tps1 's, P=' tpre . &b0 ;
      SI ITRAC1 ;
*        TRAC (def0 ET def1) 'TITR' tit1 ;
      FINSI ;
    FIN b2 ;
    SI ITRAC1 ;
      TRAC (def0 ET def1) 'TITR' tit1 ;
    FINSI ;
  FIN b0 ;

** Donnees experimentales de reference (FdP)
si (resexp);
  evxp0 = VIDE 'EVOLUTIO' ;
  tleg . (npres + 1) = MOT 'MARQ CROI NOLI' ;
  tleg . 'TITRE' . (npres + 1) = CHAI 'Ref_Exp_P=' tpre . 1 ;
  tleg . (npres + 2) = MOT 'MARQ CROI NOLI' ;
  tleg . 'TITRE' . (npres + 2) = CHAI 'Ref_Exp_P=' tpre . 2 ;
  tleg . (npres + 3) = MOT 'MARQ CROI NOLI' ;
  tleg . 'TITRE' . (npres + 3) = CHAI 'Ref_Exp_P=' tpre . 3 ;
  tleg . (npres + 4) = MOT 'MARQ CROI NOLI' ;
  tleg . 'TITRE' . (npres + 4) = CHAI 'Ref_Exp_P=' tpre . 4 ;
* Courbes Contrainte vs Deformation  
  txp  = LIRE 'CSV' FIC1 'DEBU' 2 'SEPA' ',' ;
  evxp0 = evxp0 ET (EVOL 'BLEU' 'MANU' 'EPZZ' (txp . 1) 'SMZZ' (txp . 2)) ;
  txp  = LIRE 'CSV' FIC2 'DEBU' 2 'SEPA' ',' ;
  evxp0 = evxp0 ET (EVOL 'ORAN' 'MANU' 'EPZZ' (txp . 1) 'SMZZ' (txp . 2)) ;
  txp  = LIRE 'CSV' FIC3 'DEBU' 2 'SEPA' ',' ;
  evxp0 = evxp0 ET (EVOL 'VERT' 'MANU' 'EPZZ' (txp . 1) 'SMZZ' (txp . 2)) ;
  txp  = LIRE 'CSV' FIC4 'DEBU' 2 'SEPA' ',' ;
  evxp0 = evxp0 ET (EVOL 'ROUG' 'MANU' 'EPZZ' (txp . 1) 'SMZZ' (txp . 2)) ;
fins;
** Trace des courbes pour les differents niveaux de pression
  SI ITRAC1 ;
    DESS evtot0 'TITR' (CHAI str4 ' '  '[3D V] Endommagement vs Temps') 'LEGE' 'NO' tleg ;
    SI RESEXP ;
      DESS (evtot1 ET evxp0) 'TITR' (CHAI str4 ' '  '[3D V] Contrainte vs Deformation') 'LEGE' 'NO' tleg ;
    SINON ;
      DESS evtot1 'TITR' (CHAI str4 ' '  '[3D V] Contrainte vs Deformation') 'LEGE' 'NO' tleg ;
    FINSI ;
    DESS evtot2 'TITR' (CHAI str4 ' '  '[3D V] Reaction vs Deplacement')   'LEGE' 'NO' tleg ;
  FINSI ;

** Sauvegarde des resultats de calcul Cast3M (reference pour MFront) au format .csv
  SI ISAUV ;
    str2 = '3D_V_' ;
    fic5 = CHAI str1 str2 str3 '.csv' ;
    OPTI SORT fic5 ;
    SORT 'EXCE' evtot1 'DIGI' 16 ;
  FINSI ;
  
** Comparaison Calcul / Reference sur l'aire sous la courbe Valeur absolue de l'Ecart relatif en Contrainte vs Deformation
si (resexp);
  tbid = TABL ;
  tbid . 'TITRE' = TABL ;
  tbid . 1 = MOT 'MARQ PLUS' ;
  tbid . 'TITRE' . 1 = CHAI 'Cal_mat_' str3 '_init.' ;
  tbid . 2 = MOT 'MARQ ETOI' ;
  tbid . 'TITRE' . 2 = MOT 'Ref.exp.' ;
  tbid . 3 = MOT 'MARQ CROI' ;
  tbid . 'TITRE' . 3 =  CHAI 'Cal_mat_' str3 '_ipol.' ;
  OPTI ECHO 0 ;
  SAUT 1 LIGN ;
  MESS ' **************** C A S  #1  :  3 D   V O L U M I Q U E ****************';
  SAUT 1 LIGN ;
  REPE b2 (DIME evxp0) ;
    evexp = EXTR evxp0  'COUR' &b2 ;
    evcal = EXTR evtot1 'COUR' &b2 ;
    evcou = EXTR evtot1 'COUL' &b2 ;
    xexp = ENLE (EXTR evexp 'ABSC') 1 ;
    yexp = ENLE (EXTR evexp 'ORDO') 1 ;
    xcal = ENLE (EXTR evcal 'ABSC') 1 ;
    ycal = ENLE (EXTR evcal 'ORDO') 1 ;
    xcal2 = xexp ;
    ycal2 = IPOL xcal2 evcal ;
    evcal = EVOL evcou MANU (CHAI 'EPZZ_' tpre . &b2) xcal  (CHAI 'SMZZ_' tpre . &b2) ycal ;
    evcal2 = EVOL MANU 'EPZZ' xcal2 'SMZZ' ycal2 ;
    ler1 = ABS ((ycal2 - yexp) / yexp) ;
    eer1 = EVOL 'ORAN' 'MANU' (CHAI '|EPZZ|_' tpre . &b2) (ABS xexp) (CHAI '|Er|SMZZ_' tpre . &b2) ler1 ;
    deps = ABS ((EXTR xexp (DIME xexp)) - (EXTR xexp 1)) ;
    mer1 = (INTG eer1) / deps ;
    tit1 = CHAI str4 ' '  '[3D V]' ' P=' (tpre . &b2) ' Val.abs.Ec.rel.|(Calc-Ref)/Ref|, Moy.int.|Ecart|' ' ' mer1 ;
    SI ITRAC1 ;
      DESS eer1 'TITR' tit1 ;
    FINSI ;
    SI ITRAC1 ;
      DESS (evcal et evexp et evcal2) 'TITR' tit1 'LEGE' 'NO' tbid ;
    FINSI ;
    MESS ' Pres.hydrostat.: P=' (tpre . &b2) ', Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ;
    SAUT 1 LIGN ;
  FIN b2 ;
  OPTI ECHO 1 ;
fins;
FINSI ;





************ C A S  #2  :  2 D   A X I S Y M E T R I Q U E *************
SI (EXIS lcas 2) ;
** Options generales
  OPTI 'DIME' 2 'MODE' 'AXIS' 'ELEM' 'QUA4' ;

** Maillage
  p1      = 0. 0. ;
  p2      = rayo 0. ;
  sbas    = DROI 1 p1 p2 ;
  mail    = sbas TRAN 1 (0. haut) ;
  shaut   = mail COTE 3 ;
  p3      = shaut POIN 'PROC' (0. haut) ;
  sdroite = mail COTE 2 ;
cadr1 = COUL ((0. (haut + 5.E-3)) ET (long (haut + 5.E-3))) 'BLAN' ;
rep1 = @REPERE (-0.1 -0.06) (PROG 0.15 0.15) 'TURQ' ;
SI ITRAC1 ;
  TRAC 'QUAL' (cadr1 ET rep1 ET mail ET (sbas COUL 'VERT') ET ((sdroite ET shaut) COUL 'ROUG')) 'TITR' (CHAI str4 ' '  '[2D Axi.] Maillage') ;
FINSI ;

** Modele et caracteristiques materiau (&BCALC : 1->Cast3M, 2->MFront)
SI (EGA &BCALC 1) ;
  mo = MODE mail 'MECANIQUE' 'ELASTIQUE' 'ENDOMMAGEMENT' 'MAZARS' ;
  ma = MATE mo 'YOUN' you 'NU' nu 'KTR0' epd0 'ACOM' ac 'BCOM' bc 'ATRA' at 'BTRA' bt 'BETA' beta ;
SINO ;
  coel = MOTS 'YOUN' 'NU' 'RHO' 'ALPH' 'TALP' 'TREF' 'ACOM' 'BCOM' 'ATRA' 'BTRA' 'BETA' 'EPSD';
  statev = MOTS 'DCOM' 'DTRA' 'ALFT' 'ALFC' 'D';
  mo = MODE mail 'MECANIQUE' 'ELASTIQUE' 'NON_LINEAIRE' 'UTILISATEUR' 'LIB_LOI' 'libUmatBehaviour.so' 'FCT_LOI' 'umatmazars' 'C_MATERIAU' coel 'C_VARINTER' statev;
  ma = MATE mo 'YOUN' you 'NU' nu 'RHO' xrho 'ALPH' xalph 'TALP' xtalp 'TREF' xtref 'ACOM' ac 'BCOM' bc 'ATRA' at 'BTRA' bt 'BETA' beta 'EPSD' epd0;
FINSI ;
  zer = ZERO mo 'NOEUD' ;
  un  = zer + 1. ;
  mes1 = INTG mo un 'SCAL' ;

** Blocages
  blhaut = BLOQ 'UZ' shaut ;
  blrig  = BLOQ 'UZ' sbas ;
  bl     = blhaut ET blrig ;

** Chargements (plusieurs pressions imposees, deplacement impose)
  evc1   = EVOL 'MANU' (PROG 0. 0.5 1.) (PROG 0. 1. 1.) ;
  lcoul  = MOTS 'BLEU' 'ORAN' 'VERT' 'ROUG' ;
  evtot0 = VIDE 'EVOLUTIO' ;
  evtot1 = VIDE 'EVOLUTIO' ;
  evtot2 = VIDE 'EVOLUTIO' ;
  npres = DIME lpres ;
  REPE b0 npres ;
* On fait un premier calcul avec seulement la pression imposee
    pre1 = EXTR lpres &b0 ;
*    fpr  = PRES 'MASS' mo (sdroite ET shaut) pre1 ;
    fprr = PRES 'MASS' mo sdroite pre1 ;
    fprz = PRES 'MASS' mo shaut pre1 ;
    fpr = fprr ET fprz ;
    chap = CHAR 'MECA' fpr evc1 ;
* Resolution 1 : mise en pression seulement
    t                               = TABL ;
    t . 'MODELE'                    = mo ;
    t . 'CARACTERISTIQUES'          = ma ;
    t . 'BLOCAGES_MECANIQUES'       = blrig ;
    t . 'CHARGEMENT'                = chap ;
    t . 'TEMPS_CALCULES'            = PROG 0. 'PAS' 0.1 0.5 ;
    t . 'MOVA'                      = MOT 'D' ;
    t . 'MES_SAUVEGARDES'           = TABL ;
    t . 'MES_SAUVEGARDES' . 'DEFTO' = VRAI ;
    PASAPAS t ;
* Resolution 2 : pression + ajout du deplacement impose
    dihaut = DEPI blhaut -1. ;
    uz0    = ABS (EXTR (t . 'DEPLACEMENTS' . 5) 'UZ' p3) ;
    evc2   = EVOL 'MANU' (PROG 0.5 1.) (PROG uz0 (umax - uz0)) ;
    chadi  = CHAR 'DIMP' dihaut evc2 ;
** Trace des blocages et du chargement
    re0 = @REPERE (-0.1 -0.06) (PROG 0.15 0.15) ;

    mor = MODE re0 'MECANIQUE' 'BARR' ;
    mar = MATE mor 'YOUN' 1. 'NU' 1. 'SECT' 1. ;
    rir = RIGI mor mar ;
    vprr = VECT fprr 'FORC' 'VERT' 1.E-7 ;
    vprz = VECT fprz 'FORC' 'VERT' 1.E-7 ;
    vdiz = VECT dihaut ' ' 'FLX' 'ROUG' 3.E-2 ;
    tit0 = CHAI str4 ' '  '[2D Axi.] Blocages et chargements Uimp et P=' tpre . &b0 ;
    SI ITRAC1 ;
      TRAC ((RIGI mo ma) ET bl ET rir) (vprr ET vprz ET vdiz) 'NOLE' 'TITR' tit0 ;
    FINSI ;
    t . 'BLOCAGES_MECANIQUES'       = bl ;
    t . 'CHARGEMENT'                = chap ET chadi ;
    t . 'TEMPS_CALCULES'            = PROG 0.505 'PAS' 0.005 1. ;
    PASAPAS t ;
*    SI ITRAC1 ;
*      tit1 = CHAI str4 ' '  '[2D Axi.] UZ vs Temps, P=' tpre . &b0 ;
*      DESS (EVOL 'TEMPE' t 'DEPLACEMENTS' 'UZ' p3) 'TITR' tit1 ;
*      tit1 = CHAI str4 ' '  '[2D Axi.] SMZZ vs Temps, P=' tpre . &b0 ;
*      DESS (EVOL 'TEMP' t 'CONTRAINTES' 'SMZZ' 1 1 1) 'TITR' tit1 ;
*    FINSI ;

** Post traitement
* courbes : endommagement moyen vs temps
*           contrainte moyenne vs deformation moyenne
*           force de reaction vs deplacement impose
    tt = t . 'TEMPS' ;
    tc = t . 'CONTRAINTES' ;
    te = t . 'DEFORMATIONS' ;
    tv = t . 'VARIABLES_INTERNES' ;
    tu = t . 'DEPLACEMENTS' ;
    tr = t . 'REACTIONS' ;
    def0 = DEFO (cadr1 ET rep1 ET mail) (tu . 0) 0. ;
    ltps = PROG ;
    lend = PROG ;
    lsig = PROG ;
    leps = PROG ;
    ldep = PROG ;
    lrea = PROG ;
    REPE b1 (DIME tt) ;
      tps1 = tt . (&b1 - 1) ;
      end1 = (INTG mo (tv . (&b1 - 1)) 'D')    / mes1 ;
      sig1 = (INTG mo (tc . (&b1 - 1)) 'SMZZ') / mes1 ;
      eps1 = (INTG mo (te . (&b1 - 1)) 'EPZZ') / mes1 ;
      u1   = EXTR (tu . (&b1 - 1)) 'UZ' p3 ;
      rea1 = 0. ;
      SI ((NBEL (EXTR (tr . (&b1 - 1)) 'MAIL')) > 0) ;
        rea1 = @TOTAL (tr . (&b1 - 1)) sbas 'FZ' ;
      FINSI ;
      ltps = ltps ET tps1 ;
      lend = lend ET end1 ;
      lsig = lsig ET sig1 ;
      leps = leps ET eps1 ;
      ldep = ldep ET u1 ;
      lrea = lrea ET rea1 ;
    FIN b1 ;
    evtot0 = evtot0 ET (EVOL (EXTR lcoul &b0) 'MANU' 'Temps' ltps 'D'    lend) ;
    evtot1 = evtot1 ET (EVOL (EXTR lcoul &b0) 'MANU' 'EPZZ'  leps 'SMZZ' lsig) ;
    evtot2 = evtot2 ET (EVOL (EXTR lcoul &b0) 'MANU' 'UZ'    ldep 'FZ'   lrea) ;
    tleg . 'TITRE' . &b0 = CHAI 'Cal_mat_' str3 '_P=' tpre . &b0 ;
*   deformee
    REPE b2 (DIME t . 'TEMPS') ;
      tps1 = t . 'TEMPS' . (&b2 - 1) ;
*     trace de la deformee (pour controle visuel)
      def1 = DEFO mail (tu . (&b2 - 1)) 10. 'ROUG' ;
      tit1 = CHAI str4 ' '  '[2D Axi.] Deformee au temps' ' ' tps1 's, P=' tpre . &b0 ;
      SI ITRAC1 ;
*        TRAC (def0 ET def1) 'TITR' tit1 ;
      FINSI ;
    FIN b2 ;
    SI ITRAC1 ;
      TRAC (def0 ET def1) 'TITR' tit1 ;
    FINSI ;
  FIN b0 ;

** Donnees experimentales de reference (FdP)
si (resexp);
  evxp0 = VIDE 'EVOLUTIO' ;
  tleg . (npres + 1) = MOT 'MARQ CROI NOLI' ;
  tleg . 'TITRE' . (npres + 1) = CHAI 'Ref_Exp_P=' tpre . 1 ;
  tleg . (npres + 2) = MOT 'MARQ CROI NOLI' ;
  tleg . 'TITRE' . (npres + 2) = CHAI 'Ref_Exp_P=' tpre . 2 ;
  tleg . (npres + 3) = MOT 'MARQ CROI NOLI' ;
  tleg . 'TITRE' . (npres + 3) = CHAI 'Ref_Exp_P=' tpre . 3 ;
  tleg . (npres + 4) = MOT 'MARQ CROI NOLI' ;
  tleg . 'TITRE' . (npres + 4) = CHAI 'Ref_Exp_P=' tpre . 4 ;
* Courbes Contrainte vs Deformation  
  txp  = LIRE 'CSV' FIC1 'DEBU' 2 'SEPA' ',' ;
  evxp0 = evxp0 ET (EVOL 'BLEU' 'MANU' 'EPZZ' (txp . 1) 'SMZZ' (txp . 2)) ;
  txp  = LIRE 'CSV' FIC2 'DEBU' 2 'SEPA' ',' ;
  evxp0 = evxp0 ET (EVOL 'ORAN' 'MANU' 'EPZZ' (txp . 1) 'SMZZ' (txp . 2)) ;
  txp  = LIRE 'CSV' FIC3 'DEBU' 2 'SEPA' ',' ;
  evxp0 = evxp0 ET (EVOL 'VERT' 'MANU' 'EPZZ' (txp . 1) 'SMZZ' (txp . 2)) ;
  txp  = LIRE 'CSV' FIC4 'DEBU' 2 'SEPA' ',' ;
  evxp0 = evxp0 ET (EVOL 'ROUG' 'MANU' 'EPZZ' (txp . 1) 'SMZZ' (txp . 2)) ;
fins;
** Trace des courbes pour les differents niveaux de pression
  SI ITRAC1 ;
    DESS evtot0 'TITR' (CHAI str4 ' ' '[2D Axi.] Endommagement vs Temps') 'LEGE' 'NO' tleg ;
    SI RESEXP ;
      DESS (evtot1 ET evxp0) 'TITR' (CHAI str4 ' '  '[2D Axi.] Contrainte vs Deformation') 'LEGE' 'NO' tleg ;
    SINON ;
      DESS evtot1 'TITR' (CHAI str4 ' '  '[2D Axi.] Contrainte vs Deformation') 'LEGE' 'NO' tleg ;
    FINSI ;
    DESS evtot2 'TITR' (CHAI str4 ' '  '[2D Axi.] Reaction vs Deplacement')   'LEGE' 'NO' tleg ;
  FINSI ;

** Sauvegarde des resultats de calcul Cast3M (reference pour MFront) au format .csv
  SI ISAUV ;
    str2 = '2D_axi_' ;
    fic5 = CHAI str1 str2 str3 '.csv' ;
    OPTI SORT fic5 ;
    SORT 'EXCE' evtot1 'DIGI' 16 ;
  FINSI ;
  
** Comparaison Calcul / Reference sur l'aire sous la courbe Valeur absolue de l'Ecart relatif en Contrainte vs Deformation
si (resexp);
  tbid = TABL ;
  tbid . 'TITRE' = TABL ;
  tbid . 1 = MOT 'MARQ PLUS' ;
  tbid . 'TITRE' . 1 = CHAI 'Cal_mat_' str3 '_init.' ;
  tbid . 2 = MOT 'MARQ ETOI' ;
  tbid . 'TITRE' . 2 = MOT 'Ref.exp.' ;
  tbid . 3 = MOT 'MARQ CROI' ;
  tbid . 'TITRE' . 3 =  CHAI 'Cal_mat_' str3 '_ipol.' ;
  OPTI ECHO 0 ;
  SAUT 1 LIGN ;
  MESS ' ************ C A S  #2  :  2 D   A X I S Y M E T R I Q U E *************';
  SAUT 1 LIGN ;
  REPE b2 (DIME evxp0) ;
    evexp = EXTR evxp0  'COUR' &b2 ;
    evcal = EXTR evtot1 'COUR' &b2 ;
    evcou = EXTR evtot1 'COUL' &b2 ;
    xexp = ENLE (EXTR evexp 'ABSC') 1 ;
    yexp = ENLE (EXTR evexp 'ORDO') 1 ;
    xcal = ENLE (EXTR evcal 'ABSC') 1 ;
    ycal = ENLE (EXTR evcal 'ORDO') 1 ;
    xcal2 = xexp ;
    ycal2 = IPOL xcal2 evcal ;
    evcal = EVOL evcou MANU (CHAI 'EPZZ_' tpre . &b2) xcal  (CHAI 'SMZZ_' tpre . &b2) ycal ;
    evcal2 = EVOL MANU 'EPZZ' xcal2 'SMZZ' ycal2 ;
    ler1 = ABS ((ycal2 - yexp) / yexp) ;
    eer1 = EVOL 'ORAN' 'MANU' (CHAI '|EPZZ|_' tpre . &b2) (ABS xexp) (CHAI '|Er|SMZZ_' tpre . &b2) ler1 ;
    deps = ABS ((EXTR xexp (DIME xexp)) - (EXTR xexp 1)) ;
    mer1 = (INTG eer1) / deps ;
    tit1 = CHAI str4 ' '  '[2D Axi.]' ' P=' (tpre . &b2) ' Val.abs.Ec.rel.|(Calc-Ref)/Ref|, Moy.int.|Ecart|' ' ' mer1 ;
    SI ITRAC1 ;
      DESS eer1 'TITR' tit1 ;
    FINSI ;
    SI ITRAC1 ;
      DESS (evcal et evexp et evcal2) 'TITR' tit1 'LEGE' 'NO' tbid ;
    FINSI ;
    MESS ' Pres.hydrostat.: P=' (tpre . &b2) ', Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ;
    SAUT 1 LIGN ;
  FIN b2 ;
  OPTI ECHO 1 ;
fins;
FINSI ;


FIN BCALC ;


FIN ;