* fichier : mazars_compression.dgibi
************************************************************************
*      Cas test elementaire pour la loi d'endommagement de Mazars      *
*                 Essai de compression simple monotone                 *
*                                                                      *
* Le test est realise avec plusieurs modes de calcul :                 *
*   1) 3D volumique                                                    *
*   2) 3D poutre a fibre                                               *
*   3) 2D contraintes planes                                           *
*   4) 2D axisymetrique                                                *
*   5) 2D poutre a fibre                                               *
*                                                                      *
* 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      *
************************************************************************

** 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_compression_' ;
FINSI ;

** Initialisation sortie sur Erreur 1152 "Solution calculee non conforme par rapport a la solution de reference"
IERR = FAUX ;

** Liste des modes de calcul a traiter
lcas = LECT 1 2 3 4 5 ;
*lcas = LECT 1 ; COMM '3D volumique' ;
*lcas = LECT 2 ; COMM '3D poutre a fibre' ;
*lcas = LECT 3 ; COMM '2D contraintes planes' ;
*lcas = LECT 4 ; COMM '2D axisymetrique' ;
*lcas = LECT 5 ; COMM '2D poutre a fibre' ;

** Parametres geometriques
long = 1. ;
epai = 0.1 ;

** Parametres materiau
you   = 20.E9 ;
nu    = 0.2 ;
epd0  = 8.E-5 ;
at    = 1. ;
ac    = 1.5 ;
bt    = 8000. ;
bc    = 1550. ;
beta  = 1. ;
xrho  = 0. ;
xalph = 0. ;
xtalp = 0. ;
xtref = 0. ;

** Parametres chargement
umax  = -5.E-3 ;

OPTI ECHO 0 ;
SAUT 1 lign ;
mess 'Parametres' ;
mess '==========';
SAUT 1 lign ;
mess 'Parametre de taille de l element' ;
mess '- H = arete du cube : long (m)' ' ' long ;
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 '- Depl.imp.max : umax' ' ' umax ;
SAUT 1 lign ;
OPTI ECHO 1 ;
*opti donn 5;

** Legende courbes
  tdess = 'TABLE' ;
  tdess.'TITRE' = 'TABLE' ;
  tdess. 1 = MOT 'MARQ LOSA' ;
  tdess. 2 = MOT 'MARQ CROI' ;
  tdess.'TITRE'. 2 = MOT 'Reference' ;

** 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 ' :' ;
tdess.'TITRE'. 1 = CHAI 'Cal_mat_' 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. long) ;
  mail    = sgauche VOLU 'TRAN' 1 (long 0. 0.) ;
  sdroite = mail FACE 2 ;
  p3      = sdroite POIN 'PROC' (long 0. 0.) ;
  mes1    = MESU mail ;
cadr1 = COUL ((0. 1.05 1.05) ET (1. 1.05 1.05)) 'BLAN' ;
rep1 = @REPERE (-0.55 -0.55 0.) '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') ;
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
  blgauche = BLOQ 'UX' sgauche ;
  bldroite = BLOQ 'UX' sdroite ;
  blrig    = (BLOQ 'UY' 'UZ' p1) ET (BLOQ 'UZ' p2) ;
  bl       = blgauche ET bldroite ET blrig ;

** Chargements
  didroite = DEPI bldroite umax ;
  ev1      = EVOL 'MANU' (PROG 0. 1.) (PROG 0. 1.) ;
  cha      = CHAR 'DIMP' didroite ev1 ;

** Tracé des blocages et du chargement
  re0 = @REPERE (-0.55 -0.55 0.) ;
  mor = MODE re0 'MECANIQUE' 'BARR' ;
  mar = MATE mor 'YOUN' 1. 'NU' 1. 'SECT' 1. ;
  rir = RIGI mor mar ;
  vf  = VECT didroite 'FLX' ' ' ' ' 'ROUGE' ;
  SI ITRAC1 ;
    TRAC o1 ((RIGI mo ma) ET blgauche ET blrig ET rir) vf 'NOLE' 'TITR' (CHAI str4 ' ' '[3D V] Blocages et chargement') ;
  FINSI ;

** Resolution
  t                               = TABL ;
  t . 'MODELE'                    = mo ;
  t . 'CARACTERISTIQUES'          = ma ;
  t . 'BLOCAGES_MECANIQUES'       = bl ;
  t . 'CHARGEMENT'                = cha ;
  t . 'TEMPS_CALCULES'            = PROG 0. 'PAS' 0.01 1. ;
  t . 'MOVA'                      = MOT 'D' ;
  t . 'MES_SAUVEGARDES'           = TABL ;
  t . 'MES_SAUVEGARDES' . 'DEFTO' = VRAI ;
  PASAPAS t ;

** 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) ;
    ltps = ltps ET tps1 ;
    end1 = (INTG mo (tv . (&b1 - 1)) 'D')    / mes1 ;
    sig1 = (INTG mo (tc . (&b1 - 1)) 'SMXX') / mes1 ;
    eps1 = (INTG mo (te . (&b1 - 1)) 'EPXX') / mes1 ;
    u1   = EXTR (tu . (&b1 - 1)) 'UX' p3 ;
    rea1 = 0. ;
    SI (NEG &b1 1) ;
      rea1 = @TOTAL (tr . (&b1 - 1)) sgauche 'FX' ;
    FINSI ;
    lend = lend ET end1 ;
    lsig = lsig ET sig1 ;
    leps = leps ET eps1 ;
    ldep = ldep ET u1 ;
    lrea = lrea ET rea1 ;
  FIN b1 ;
  ev0 = EVOL 'ROUG' 'MANU' 'Temps' ltps 'D'    lend ;
  ev1 = EVOL 'ROUG' 'MANU' 'EPXX'  leps 'SMXX' lsig ;
  ev2 = EVOL 'ROUG' 'MANU' 'UX'    ldep 'FX'   lrea ;
  ltps1 =  ENLE ltps 1 ;
  nb1 = DIME ltps1 ;
  un1 = PROG nb1 * 1.D0 ;
  lepd0 = PROG nb1 * epd0 ;
  lac1 = PROG nb1 * ac ;
  lepu = ltps1 * umax / long ;
  lepe = (2.**0.5) * nu * (ABS lepu) ;
  dc = un1 - ((epd0*(un1-lac1)/lepe) + (ac*(EXP(bc*(lepd0-lepe))))) ;
*mess 'dc avant filtrage';list dc;
  lfiltr = LECT ;
  REPE b2 ;
    dc_i = EXTR dc &b2 ;
    SI (EGA &b2 1) ;
      dc_p = 1. ;
    FINSI ;
    ddc = dc_i - dc_p ;
*mess '=== Boucle b2 passage' ' ' &b2;
*mess '    dom.courant' ' ' dc_i;
*mess '    inc.dommage' ' ' ddc;
    SI ((ddc < 0.) OU (dc_i < 0.) OU (dc_i > 1.)) ;
*mess '    => a filtrer !!!';
      lfiltr = lfiltr ET (LECT &b2) ;
*sino;
*mess '    a garder...';
    FINSI ;
    dc_p = dc_i ;
    SI (EGA &b2 (DIME dc)) ;
*mess '=== Sortie Boucle b2';
*mess '    liste des indices a eliminer';list lfiltr;
      QUIT b2 ;
    FINSI ;
  FIN b2 ;
  ltps1 = ENLE ltps1 lfiltr ;
  un1   = ENLE un1   lfiltr ;
  lepu  = ENLE lepu  lfiltr ;
  dc    = ENLE dc    lfiltr ;
  lsig1 = (un1 - dc) * you * lepu ;
  an0 = EVOL 'TURQ' 'MANU' 'Temps' ltps1 'D'    dc    ;
  an1 = EVOL 'TURQ' 'MANU' 'EPXX'  lepu  'SMXX' lsig1 ;
  ansmxx = EVOL 'TURQ' 'MANU' 'Temps' ltps1 'SMXX' lsig1 ;
  SI ITRAC1 ;
    DESS (ev0 ET an0) 'TITR' (CHAI str4 ' ' '[3D V] Endommagement vs Temps') 'LEGE' 'SE' tdess ;
    DESS (ev1 ET an1) 'TITR' (CHAI str4 ' ' '[3D V] Contrainte vs Deformation') 'LEGE' 'NE' tdess ;
    DESS ev2 'TITR' (CHAI str4 ' ' '[3D V] Reaction vs Deplacement') 'LEGE' 'SE' tdess ;
  FINSI ;
* evolutions temporelles a tous les points de Gauss
  ng = 8 ;
  tleg = TABL ;
  tleg . 1 = MOT 'MARQ CARR NOLI' ;
  tleg . 2 = MOT 'MARQ LOSA NOLI' ;
  tleg . 3 = MOT 'MARQ ROND NOLI' ;
  tleg . 4 = MOT 'MARQ ETOI NOLI' ;
  tleg . 5 = MOT 'MARQ TRID NOLI' ;
  tleg . 6 = MOT 'MARQ TRIU NOLI' ;
  tleg . 7 = MOT 'MARQ TRIL NOLI' ;
  tleg . 8 = MOT 'MARQ TRIR NOLI' ;
  tleg . 'TITRE' = TABL ;
*
  evend  = VIDE 'EVOLUTIO' ;
  evsmxx = VIDE 'EVOLUTIO' ;
  evepxx = VIDE 'EVOLUTIO' ;
  REPE b ng ;
    tleg . 'TITRE' . &b = CHAI 'Cal_mat_' str3 '_Pt_Gauss' ' ' &b ;
    evend  = evend  ET (EVOL 'ROSE' 'TEMP' t 'VARIABLES_INTERNES' 'D'    1 1 &b) ;
    evsmxx = evsmxx ET (EVOL 'BLEU' 'TEMP' t 'CONTRAINTES'        'SMXX' 1 1 &b) ;
    evepxx = evepxx ET (EVOL 'VERT' 'TEMP' t 'DEFORMATIONS'       'EPXX' 1 1 &b) ;
  FIN b ;
*
  tleg . (ng+1) = MOT 'MARQ CROI' ;
  tleg . 'TITRE' . (ng+1) = MOT 'Reference' ;
*
  SI ITRAC1 ;
    DESS (evend ET an0)     'LEGE' 'SE' tleg 'TITR' (CHAI str4 ' ' '[3D V] D vs Temps') ;
    DESS (evsmxx ET ansmxx) 'LEGE' 'NE' tleg 'TITR' (CHAI str4 ' ' '[3D V] SMXX vs Temps') ;
    DESS evepxx             'LEGE' 'NE' tleg 'TITR' (CHAI str4 ' ' '[3D V] EPXX vs Temps') ;
  FINSI ;
* deformee
  REPE b2 (DIME tt) ;
    tps1 = tt . (&b2 - 1) ;
*   trace de la deformee (pour controle visuel)
    def1 = DEFO mail (tu . (&b2 - 1)) 33. 'ROUG' ;
    tit1 = CHAI str4 ' ' '[3D V] Deformee au temps' ' ' tps1 ;
    SI ITRAC1 ;
*      TRAC o1 (def0 ET def1) 'TITR' tit1 ;
    FINSI ;
  FIN b2 ;
  SI ITRAC1 ;
    TRAC o1 (def0 ET def1) 'TITR' tit1 ;
  FINSI ;
* comparaison Calcul / Reference sur la contrainte moyenne
  ler1 = ABS (((ENLE (ENLE lsig 1) lfiltr) - lsig1) / lsig1) ;
  eer1 = EVOL 'ORAN' 'MANU' 'Temps'  ltps1  '|Ec.rel|SMXX' ler1 ;
  dtps = (EXTR ltps1 (DIME ltps1)) - (EXTR ltps1 1) ;
  mer1 = (INTG eer1) / dtps ;
  tit1 = CHAI str4 ' ' '[3D V] Val.abs.Ec.rel.|(Calc-Ref)/Ref|, Moy.int.|Ecart|' ' ' mer1 ;
  SI ITRAC1 ;
    DESS eer1 'TITR' tit1 ;
  FINSI ;

** Sauvegarde des resultats de calcul Cast3M (reference pour MFront) au format .csv
  SI ISAUV ;
    str2 = '3D_V_' ;
    fic1 = CHAI str1 str2 str3 '.csv' ;
    OPTI SORT fic1 ;
    SORT 'EXCE' ev1 'DIGI' 16 ;
  FINSI ;
  
  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 ;
  SI (mer1 < 1.E-6) ;
    MESS ' Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ' < 1.E-06 => OK' ;
    SAUT 1 LIGN ;
  SINO ;
    MESS ' Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ' > 1.E-06 => ECHEC' ;
    IERR = VRAI ;
    SAUT 1 LIGN ;
  FINS ;
  SAUT 1 LIGN ;
  OPTI ECHO 1 ;
FINSI ;





*********** C A S  #2  :  3 D   P O U T R E   A   F I B R E ************
SI ((EXIS lcas 2) ET (EGA &BCALC 1)) ;

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

** Maillage
  p1s     = (-0.5 * epai) (-0.5 * epai) 0. ;
  p2s     = (-0.5 * epai) ( 0.5 * epai) 0. ;
  l12s    = DROI 1 p1s p2s ;
  mails   = l12s TRAN 1 (epai 0. 0.) ;
  p1      = 0. 0. 0. ;
  p2      = long 0. 0. ;
  mail    = DROI 1 p1 p2 ;
  mes1    = MESU mails ;
cadr1 = COUL ((-0.1 0. 0.) ET (1. 0. 0.)) 'BLAN' ;
rep1 = @REPERE (-0.2 -0.2 0.) (PROG 0.2 0.2 0.2) 'TURQ' ;
SI ITRAC1 ;
  TRAC o1 'QUAL' (cadr1 ET rep1 ET (mails COUL 'VERT') ET (mail COUL 'ROUG')) 'TITR' (CHAI str4 ' ' '[3D PaF] Maillages de la Poutre et de sa Section') ;
FINSI ;

** Modele et caracteristiques materiau (&BCALC : 1->Cast3M, 2->MFront)
SI (EGA &BCALC 1) ;
  mos = MODE mails 'MECANIQUE' 'ELASTIQUE' 'PLASTIQUE' 'MAZARS' 'QUAS' ;
  mas = MATE mos 'YOUN' you 'NU' nu 'KTR0' epd0 'ACOM' ac 'BCOM' bc 'ATRA' at 'BTRA' bt 'BETA' beta
                 'ALPY' 1. 'ALPZ' 1. ;
  mo = MODE mail 'MECANIQUE' 'ELASTIQUE' 'SECTION' 'PLASTIQUE' 'SECTION' 'TIMO' ;
  ma = MATE mo 'MODS' mos 'MATS' mas 'VECT' (0. 1. 0.) ;
SINO ;
* Implémentation MFront A FAIRE. En attendant on quitte la boucle...
  QUIT BCALC;
FINSI ;

** Blocages
  blgauche = BLOQ 'DEPL' 'ROTA' p1 ;
  bldroite = BLOQ 'UX' p2 ;
  bl       = blgauche ET bldroite ;
  bldroit2 = BLOQ 'ROTA' p2 ;
  bl       = bl ET bldroit2 ;

** Chargements
  didroite = DEPI bldroite umax ;
  ev1      = EVOL 'MANU' (PROG 0. 1.) (PROG 0. 1.) ;
  cha      = CHAR 'DIMP' didroite ev1 ;

** Tracé des blocages et du chargement
  re0 = @REPERE (-0.2 -0.2 0.) (PROG 0.2 0.2 0.2) ;
  mor = MODE re0 'MECANIQUE' 'BARR' ;
  mar = MATE mor 'YOUN' 1. 'NU' 1. 'SECT' 1. ;
  rir = RIGI mor mar ;
  vf  = VECT didroite 'FLX' ' ' ' ' 'ROUG' 20. ;
  SI ITRAC1 ;
    TRAC o1 ((RIGI mo ma) ET blgauche ET rir) vf 'NOLE' 'TITR' (CHAI str4 ' ' '[3D PaF] Blocages et chargement') ;
  FINSI ;

** Resolution
  t                               = TABL ;
  t . 'MODELE'                    = mo ;
  t . 'CARACTERISTIQUES'          = ma ;
  t . 'BLOCAGES_MECANIQUES'       = bl ;
  t . 'CHARGEMENT'                = cha ;
  t . 'TEMPS_CALCULES'            = PROG 0. 'PAS' 0.01 1. ;
  t . 'MES_SAUVEGARDES'           = TABL ;
  t . 'MES_SAUVEGARDES' . 'DEFTO' = VRAI ;
  PASAPAS t ;

** Post traitement
* courbes : endommagement moyen vs temps
*           contrainte moyenne vs deformation moyenne
*           force de reaction vs deplacement impose
* evolutions temporelles a tous les points de Gauss
  ng = 4 ;
  tleg = TABL ;
  tleg . 1 = MOT 'MARQ CARR NOLI' ;
  tleg . 2 = MOT 'MARQ LOSA NOLI' ;
  tleg . 3 = MOT 'MARQ ROND NOLI' ;
  tleg . 4 = MOT 'MARQ ETOI NOLI' ;
  tleg . 'TITRE' = TABL ;
  REPE b ng ;
    tleg . 'TITRE' . &b = CHAI 'Cal_mat_' str3 '_Pt_Gauss' ' ' &b ;
  FIN b ;
  tleg . (ng+1) = MOT 'MARQ CROI' ;
  tleg . 'TITRE' . (ng+1) = MOT 'Reference' ;
  tt = t . 'TEMPS' ;
  tv = t . 'VARIABLES_INTERNES' ;
  tu = t . 'DEPLACEMENTS' ;
  tr = t . 'REACTIONS' ;
  def0 = DEFO (cadr1 ET rep1 ET mail) (tu . 0) 0. ;
  ltps = PROG ;
  lend1 = PROG ;
  lend2 = PROG ;
  lend3 = PROG ;
  lend4 = PROG ;
  lsig1 = PROG ;
  lsig2 = PROG ;
  lsig3 = PROG ;
  lsig4 = PROG ;
  leps1 = PROG ;
  leps2 = PROG ;
  leps3 = PROG ;
  leps4 = PROG ;
  ldep = PROG ;
  lrea = PROG ;
  REPE b1 (DIME tt) ;
    tps1 = tt . (&b1 - 1) ;
    ltps = ltps ET tps1 ;
    SI (EGA &b1 1) ;
      lend1 = lend1 ET 0. ;
      lend2 = lend2 ET 0. ;
      lend3 = lend3 ET 0. ;
      lend4 = lend4 ET 0. ;
      lsig1 = lsig1 ET 0. ;
      lsig2 = lsig2 ET 0. ;
      lsig3 = lsig3 ET 0. ;
      lsig4 = lsig4 ET 0. ;
      leps1 = leps1 ET 0. ;
      leps2 = leps2 ET 0. ;
      leps3 = leps3 ET 0. ;
      leps4 = leps4 ET 0. ;
    SINON ;
      lend1 = lend1 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VHNI' 1 1 1) ;
      lend2 = lend2 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VHNI' 1 1 2) ;
      lend3 = lend3 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VHNI' 1 1 3) ;
      lend4 = lend4 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VHNI' 1 1 4) ;
      lsig1 = lsig1 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VONS' 1 1 1) 'SMXX' 1 1 1) ;
      lsig2 = lsig2 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VONS' 1 1 1) 'SMXX' 1 1 2) ;
      lsig3 = lsig3 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VONS' 1 1 1) 'SMXX' 1 1 3) ;
      lsig4 = lsig4 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VONS' 1 1 1) 'SMXX' 1 1 4) ;
      leps1 = leps1 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VISX' 1 1 1) ;
      leps2 = leps2 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VISX' 1 1 2) ;
      leps3 = leps3 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VISX' 1 1 3) ;
      leps4 = leps4 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VISX' 1 1 4) ;
    FINSI ;
    u1   = EXTR (tu . (&b1 - 1)) 'UX' p2 ;
    rea1 = EXTR (tr . (&b1 - 1)) 'FX' p1 ;
    ldep = ldep ET u1 ;
    lrea = lrea ET rea1 ;
  FIN b1 ;
  ev0 = (EVOL 'ROUG' 'MANU' 'Temps' ltps 'D' lend1) ET (EVOL 'ROUG' 'MANU' 'Temps' ltps 'D' lend2) ET
        (EVOL 'ROUG' 'MANU' 'Temps' ltps 'D' lend3) ET (EVOL 'ROUG' 'MANU' 'Temps' ltps 'D' lend4) ;
  ev1 = (EVOL 'ROUG' 'MANU' 'EPXX'  leps1 'SMXX' lsig1) ET (EVOL 'ROUG' 'MANU' 'EPXX' leps2 'SMXX' lsig2) ET
        (EVOL 'ROUG' 'MANU' 'EPXX'  leps3 'SMXX' lsig3) ET (EVOL 'ROUG' 'MANU' 'EPXX' leps4 'SMXX' lsig4) ;
  ev2 = EVOL 'ROUG' 'MANU' 'UX' ldep 'FX' lrea ;
  ltps1 =  ENLE ltps 1 ;
  nb1 = DIME ltps1 ;
  un1 = PROG nb1 * 1.D0 ;
  lepd0 = PROG nb1 * epd0 ;
  lac1 = PROG nb1 * ac ;
  lepu = ltps1 * umax / long ;
  lepe = (2.**0.5) * nu * (ABS lepu) ;
  dc = un1 - ((epd0*(un1-lac1)/lepe) + (ac*(EXP(bc*(lepd0-lepe))))) ;
  lfiltr = LECT ;
  REPE b2 ;
    dc_i = EXTR dc &b2 ;
    SI (EGA &b2 1) ;
      dc_p = 1. ;
    FINSI ;
    ddc = dc_i - dc_p ;
    SI ((ddc < 0.) OU (dc_i < 0.) OU (dc_i > 1.)) ;
      lfiltr = lfiltr ET (LECT &b2) ;
    FINSI ;
    dc_p = dc_i ;
    SI (EGA &b2 (DIME dc)) ;
      QUIT b2 ;
    FINSI ;
  FIN b2 ;
  ltps1 = ENLE ltps1 lfiltr ;
  un1   = ENLE un1   lfiltr ;
  lepu  = ENLE lepu  lfiltr ;
  dc    = ENLE dc    lfiltr ;
  lsan1 = (un1 - dc) * you * lepu ;
  an0 = EVOL 'TURQ' 'MANU' 'Temps' ltps1 'D'    dc    ;
  an1 = EVOL 'TURQ' 'MANU' 'EPXX'  lepu  'SMXX' lsan1 ;
  SI ITRAC1 ;
    DESS (ev0 ET an0) 'LEGE' 'SE' tleg 'TITR' (CHAI str4 ' ' '[3D PaF] Endommagement vs Temps') ;
    DESS (ev1 ET an1) 'LEGE' 'NE' tleg 'TITR' (CHAI str4 ' ' '[3D PaF] Contrainte vs Deformation') ;
    DESS ev2 'TITR' (CHAI str4 ' ' '[3D PaF] Reaction vs Deplacement') 'LEGE' 'SE' tdess ;
  FINSI ;
* deformee
  REPE b2 (DIME tt) ;
    tps1 = tt . (&b2 - 1) ;
*   trace de la deformee (pour controle visuel)
    def1 = DEFO mail (tu . (&b2 - 1)) 33. 'ROUG' ;
    tit1 = CHAI str4 ' ' '[3D PaF] Deformee au temps' ' ' tps1 ;
    SI ITRAC1 ;
*      TRAC o1 (def0 ET def1) 'TITR' tit1 ;
    FINSI ;
  FIN b2 ;
  SI ITRAC1 ;
    TRAC o1 (def0 ET def1) 'TITR' tit1 ;
  FINSI ;
* comparaison Calcul / Reference sur la contrainte moyenne
  lsig = (lsig1 + lsig2 + lsig3 + lsig4) / 4. ;
  ler1 = ABS (((ENLE (ENLE lsig 1) lfiltr) - lsan1) / lsan1) ;
  eer1 = EVOL 'ORAN' 'MANU' 'Temps'  ltps1  '|Ec.rel|SMXX' ler1 ;
  dtps = (EXTR ltps1 (DIME ltps1)) - (EXTR ltps1 1) ;
  mer1 = (INTG eer1) / dtps ;
  tit1 = CHAI str4 ' ' '[3D PaF] Val.abs.Ec.rel.|(Calc-Ref)/Ref|, Moy.int.|Ecart|' ' ' mer1 ;
  SI ITRAC1 ;
    DESS eer1 'TITR' tit1 ;
  FINSI ;
  OPTI ECHO 0 ;
  SAUT 1 LIGN ;
  MESS ' *********** C A S  #2  :  3 D   P O U T R E   A   F I B R E ************';
  SAUT 1 LIGN ;
  SI (mer1 < 1.E-6) ;
    MESS ' Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ' < 1.E-06 => OK' ;
    SAUT 1 LIGN ;
  SINO ;
    MESS ' Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ' > 1.E-06 => ECHEC' ;
    IERR = VRAI ;
    SAUT 1 LIGN ;
  FINS ;
  SAUT 1 LIGN ;
  OPTI ECHO 1 ;
FINSI ;





******* C A S  #3  :  2 D   C O N T R A I N T E S   P L A N E S ********
SI (EXIS lcas 3) ;

** Options generales
  OPTI 'DIME' 2 'MODE' 'PLAN' 'CONT' 'ELEM' 'QUA4' ;

** Maillage
  p1      = 0. 0. ;
  p2      = 0. long ;
  sgauche = DROI 1 p1 p2 ;
  mail    = sgauche TRAN 1 (long 0.) ;
  sdroite = mail COTE 3 ;
  p3      = sdroite POIN 'PROC' (long 0.) ;
  mes1    = MESU mail ;
cadr1 = COUL ((0. 1.05) ET (1. 1.05)) 'BLAN' ;
rep1 = @REPERE (-0.2 -0.2) (PROG 0.5 0.5) 'TURQ' ;
SI ITRAC1 ;
  TRAC 'QUAL' (cadr1 ET rep1 ET mail ET (sgauche COUL 'VERT') ET (sdroite COUL 'ROUG')) 'TITR' (CHAI str4 ' ' '[2D Plan cont.] 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 'DIM3' epai ;
SINO ;
  coel = MOTS 'YOUN' 'NU' 'RHO' 'ALPH' 'DIM3' 'TALP' 'TREF' 'ACOM' 'BCOM' 'ATRA' 'BTRA' 'BETA' 'EPSD';
  statev = MOTS 'DCOM' 'DTRA' 'ALFT' 'ALFC' 'D' 'ETOZ';
  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 'DIM3' epai 'TALP' xtalp 'TREF' xtref 'ACOM' ac 'BCOM' bc 'ATRA' at 'BTRA' bt 'BETA' beta 'EPSD' epd0;
FINSI ;

** Blocages
  blgauche = BLOQ 'UX' sgauche ;
  bldroite = BLOQ 'UX' sdroite ;
  blrig    = BLOQ 'UY' p1 ;
  bl       = blgauche ET bldroite ET blrig ;

** Chargements
  didroite = DEPI bldroite umax ;
  ev1      = EVOL 'MANU' (PROG 0. 1.) (PROG 0. 1.) ;
  cha      = CHAR 'DIMP' didroite ev1 ;

** Tracé des blocages et du chargement
  re0 = @REPERE (-0.2 -0.2) (PROG 0.5 0.5) ;
  mor = MODE re0 'MECANIQUE' 'BARR' ;
  mar = MATE mor 'YOUN' 1. 'NU' 1. 'SECT' 1. ;
  rir = RIGI mor mar ;
  vf  = VECT didroite 'FLX' ' ' 'ROUG' ;
  SI ITRAC1 ;
    TRAC ((RIGI mo ma) ET blgauche ET blrig ET rir) vf 'NOLE' 'TITR' (CHAI str4 ' ' '[2D Plan cont.] Blocages et chargement') ;
  FINSI ;

** Resolution
  t                               = TABL ;
  t . 'MODELE'                    = mo ;
  t . 'CARACTERISTIQUES'          = ma ;
  t . 'BLOCAGES_MECANIQUES'       = bl ;
  t . 'CHARGEMENT'                = cha ;
  t . 'TEMPS_CALCULES'            = PROG 0. 'PAS' 0.01 1. ;
  t . 'MOVA'                      = MOT 'D' ;
  t . 'MES_SAUVEGARDES'           = TABL ;
  t . 'MES_SAUVEGARDES' . 'DEFTO' = VRAI ;
  PASAPAS t ;

** 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) ;
    ltps = ltps ET tps1 ;
    end1 = (INTG mo (tv . (&b1 - 1)) 'D')    / mes1 ;
    sig1 = (INTG mo (tc . (&b1 - 1)) 'SMXX') / mes1 ;
    eps1 = (INTG mo (te . (&b1 - 1)) 'EPXX') / mes1 ;
    u1   = EXTR (tu . (&b1 - 1)) 'UX' p3 ;
    rea1 = 0. ;
    SI (NEG &b1 1) ;
      rea1 = @TOTAL (tr . (&b1 - 1)) sgauche 'FX' ;
    FINSI ;
    lend = lend ET end1 ;
    lsig = lsig ET sig1 ;
    leps = leps ET eps1 ;
    ldep = ldep ET u1 ;
    lrea = lrea ET rea1 ;
  FIN b1 ;
  ev0 = EVOL 'ROUG' 'MANU' 'Temps' ltps 'D'    lend ;
  ev1 = EVOL 'ROUG' 'MANU' 'EPXX'  leps 'SMXX' lsig ;
  ev2 = EVOL 'ROUG' 'MANU' 'UX'    ldep 'FX'   lrea ;
  ltps1 =  ENLE ltps 1 ;
  nb1 = DIME ltps1 ;
  un1 = PROG nb1 * 1.D0 ;
  lepd0 = PROG nb1 * epd0 ;
  lac1 = PROG nb1 * ac ;
  lepu = ltps1 * umax / long ;
  lepe = (2.**0.5) * nu * (ABS lepu) ;
  dc = un1 - ((epd0*(un1-lac1)/lepe) + (ac*(EXP(bc*(lepd0-lepe))))) ;
  lfiltr = LECT ;
  REPE b2 ;
    dc_i = EXTR dc &b2 ;
    SI (EGA &b2 1) ;
      dc_p = 1. ;
    FINSI ;
    ddc = dc_i - dc_p ;
    SI ((ddc < 0.) OU (dc_i < 0.) OU (dc_i > 1.)) ;
      lfiltr = lfiltr ET (LECT &b2) ;
    FINSI ;
    dc_p = dc_i ;
    SI (EGA &b2 (DIME dc)) ;
      QUIT b2 ;
    FINSI ;
  FIN b2 ;
  ltps1 = ENLE ltps1 lfiltr ;
  un1   = ENLE un1   lfiltr ;
  lepu  = ENLE lepu  lfiltr ;
  dc    = ENLE dc    lfiltr ;
  lsig1 = (un1 - dc) * you * lepu ;
  an0 = EVOL 'TURQ' 'MANU' 'Temps' ltps1 'D'    dc    ;
  an1 = EVOL 'TURQ' 'MANU' 'EPXX'  lepu  'SMXX' lsig1 ;
  ansmxx = EVOL 'TURQ' 'MANU' 'Temps' ltps1 'SMXX' lsig1 ;
  SI ITRAC1 ;
    DESS (ev0 ET an0) 'TITR' (CHAI str4 ' ' '[2D Plan cont.] Endommagement vs Temps') 'LEGE' 'SE' tdess ;
    DESS (ev1 ET an1) 'TITR' (CHAI str4 ' ' '[2D Plan cont.] Contrainte vs Deformation') 'LEGE' 'NE' tdess ;
    DESS ev2 'TITR' (CHAI str4 ' ' '[2D Plan cont.] Reaction vs Deplacement') 'LEGE' 'SE' tdess ;
  FINSI ;
* evolutions temporelles a tous les points de Gauss
  ng = 4 ;
  tleg = TABL ;
  tleg . 1 = MOT 'MARQ CARR NOLI' ;
  tleg . 2 = MOT 'MARQ LOSA NOLI' ;
  tleg . 3 = MOT 'MARQ ROND NOLI' ;
  tleg . 4 = MOT 'MARQ ETOI NOLI' ;
  tleg . 'TITRE' = TABL ;
*
  evend = VIDE 'EVOLUTIO' ;
  evsmxx = VIDE 'EVOLUTIO' ;
  evepxx = VIDE 'EVOLUTIO' ;
  REPE b ng ;
    tleg . 'TITRE' . &b = CHAI 'Cal_mat_' str3 '_Pt_Gauss' ' ' &b ;
    evend = evend ET (EVOL 'ROSE' 'TEMP' t 'VARIABLES_INTERNES' 'D' 1 1 &b) ;
    evsmxx = evsmxx ET (EVOL 'BLEU' 'TEMP' t 'CONTRAINTES' 'SMXX' 1 1 &b) ;
    evepxx = evepxx ET (EVOL 'VERT' 'TEMP' t 'DEFORMATIONS' 'EPXX' 1 1 &b) ;
  FIN b ;
*
  tleg . (ng+1) = MOT 'MARQ CROI' ;
  tleg . 'TITRE' . (ng+1) = MOT 'Reference' ;
*
  SI ITRAC1 ;
    DESS (evend ET an0) 'LEGE' 'SE' tleg 'TITR' (CHAI str4 ' ' '[2D Plan cont.] D vs Temps') ;
    DESS (evsmxx ET ansmxx) 'LEGE' 'NE' tleg 'TITR' (CHAI str4 ' ' '[2D Plan cont.] SMXX vs Temps') ;
    DESS evepxx 'LEGE' 'NE' tleg 'TITR' (CHAI str4 ' ' '[2D Plan cont.] EPXX vs Temps') ;
  FINSI ;
* deformee
  REPE b2 (DIME tt) ;
    tps1 = tt . (&b2 - 1) ;
*   trace de la deformee (pour controle visuel)
    def1 = DEFO mail (tu . (&b2 - 1)) 33. 'ROUG' ;
    tit1 = CHAI str4 ' ' '[2D Plan cont.] Deformee au temps' ' ' tps1 ;
    SI ITRAC1 ;
*      TRAC (def0 ET def1) 'TITR' tit1 ;
    FINSI ;
  FIN b2 ;
  SI ITRAC1 ;
    TRAC (def0 ET def1) 'TITR' tit1 ;
  FINSI ;
* comparaison Calcul / Reference sur la contrainte moyenne
  ler1 = ABS (((ENLE (ENLE lsig 1) lfiltr) - lsig1) / lsig1) ;
  eer1 = EVOL 'ORAN' 'MANU' 'Temps'  ltps1  '|Ec.rel|SMXX' ler1 ;
  dtps = (EXTR ltps1 (DIME ltps1)) - (EXTR ltps1 1) ;
  mer1 = (INTG eer1) / dtps ;
  tit1 = CHAI str4 ' ' '[2D Plan cont.] Val.abs.Ec.rel.|(Calc-Ref)/Ref|, Moy.int.|Ecart|' ' ' mer1 ;
  SI ITRAC1 ;
    DESS eer1 'TITR' tit1 ;
  FINSI ;

** Sauvegarde des resultats de calcul Cast3M (reference pour MFront) au format .csv
  SI ISAUV ;
    str2 = '2D_CP_' ;
    fic1 = CHAI str1 str2 str3 '.csv' ;
    OPTI SORT fic1 ;
    SORT 'EXCE' ev1 'DIGI' 16 ;
  FINSI ;
  
  OPTI ECHO 0 ;
  SAUT 1 LIGN ;
  MESS ' ******* C A S  #3  :  2 D   C O N T R A I N T E S   P L A N E S ********';
  SAUT 1 LIGN ;
  SI (mer1 < 1.E-6) ;
    MESS ' Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ' < 1.E-06 => OK' ;
    SAUT 1 LIGN ;
  SINO ;
    MESS ' Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ' > 1.E-06 => ECHEC' ;
    IERR = VRAI ;
    SAUT 1 LIGN ;
  FINS ;
  SAUT 1 LIGN ;
  OPTI ECHO 1 ;
FINSI ;





************ C A S  #4  :  2 D   A X I S Y M E T R I Q U E *************
SI (EXIS lcas 4) ;

** Options generales
  OPTI 'DIME' 2 'MODE' 'AXIS' 'ELEM' 'QUA4' ;

** Maillage
  p1      = 0. 0. ;
  p2      = long 0. ;
  sbas    = DROI 1 p1 p2 ;
  mail    = sbas TRAN 1 (0. long) ;
  shaut   = mail COTE 3 ;
  p3      = shaut POIN 'PROC' (0. long) ;
cadr1 = COUL ((0. 1.) ET (1.05 1.)) 'BLAN' ;
rep1 = @REPERE (-0.2 -0.2) (PROG 0.5 0.5) 'TURQ' ;
SI ITRAC1 ;
  TRAC 'QUAL' (cadr1 ET rep1 ET mail ET (sbas COUL 'VERT') 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
  blbas    = BLOQ 'UZ' sbas ;
  blhaut   = BLOQ 'UZ' shaut ;
  bl       = blbas ET blhaut ;

** Chargements
  dihaut   = DEPI blhaut umax ;
  ev1      = EVOL 'MANU' (PROG 0. 1.) (PROG 0. 1.) ;
  cha      = CHAR 'DIMP' dihaut ev1 ;

** Tracé des blocages et du chargement
  re0 = @REPERE (-0.2 -0.2) (PROG 0.5 0.5) ;
  mor = MODE re0 'MECANIQUE' 'BARR' ;
  mar = MATE mor 'YOUN' 1. 'NU' 1. 'SECT' 1. ;
  rir = RIGI mor mar ;
  vf  = VECT dihaut ' ' 'FLX' 'ROUG' ;
  SI ITRAC1 ;
    TRAC ((RIGI mo ma) ET blbas ET rir) vf 'NOLE' 'TITR' (CHAI str4 ' ' '[2D Axi.] Blocages et chargement') ;
  FINSI ;

** Resolution
  t                               = TABL ;
  t . 'MODELE'                    = mo ;
  t . 'CARACTERISTIQUES'          = ma ;
  t . 'BLOCAGES_MECANIQUES'       = bl ;
  t . 'CHARGEMENT'                = cha ;
  t . 'TEMPS_CALCULES'            = PROG 0. 'PAS' 0.01 1. ;
  t . 'MOVA'                      = MOT 'D' ;
  t . 'MES_SAUVEGARDES'           = TABL ;
  t . 'MES_SAUVEGARDES' . 'DEFTO' = VRAI ;
  PASAPAS t ;

** 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) ;
    ltps = ltps ET tps1 ;
    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 (NEG &b1 1) ;
      rea1 = @TOTAL (tr . (&b1 - 1)) sbas 'FZ' ;
    FINSI ;
    lend = lend ET end1 ;
    lsig = lsig ET sig1 ;
    leps = leps ET eps1 ;
    ldep = ldep ET u1 ;
    lrea = lrea ET rea1 ;
  FIN b1 ;
  ev0 = EVOL 'ROUG' 'MANU' 'Temps' ltps 'D'    lend ;
  ev1 = EVOL 'ROUG' 'MANU' 'EPZZ'  leps 'SMZZ' lsig ;
  ev2 = EVOL 'ROUG' 'MANU' 'UZ'    ldep 'FZ'   lrea ;
  ltps1 =  ENLE ltps 1 ;
  nb1 = DIME ltps1 ;
  un1 = PROG nb1 * 1.D0 ;
  lepd0 = PROG nb1 * epd0 ;
  lac1 = PROG nb1 * ac ;
  lepu = ltps1 * umax / long ;
  lepe = (2.**0.5) * nu * (ABS lepu) ;
  dc = un1 - ((epd0*(un1-lac1)/lepe) + (ac*(EXP(bc*(lepd0-lepe))))) ;
  lfiltr = LECT ;
  REPE b2 ;
    dc_i = EXTR dc &b2 ;
    SI (EGA &b2 1) ;
      dc_p = 1. ;
    FINSI ;
    ddc = dc_i - dc_p ;
    SI ((ddc < 0.) OU (dc_i < 0.) OU (dc_i > 1.)) ;
      lfiltr = lfiltr ET (LECT &b2) ;
    FINSI ;
    dc_p = dc_i ;
    SI (EGA &b2 (DIME dc)) ;
      QUIT b2 ;
    FINSI ;
  FIN b2 ;
  ltps1 = ENLE ltps1 lfiltr ;
  un1   = ENLE un1   lfiltr ;
  lepu  = ENLE lepu  lfiltr ;
  dc    = ENLE dc    lfiltr ;
  lsig1 = (un1 - dc) * you * lepu ;
  an0 = EVOL 'TURQ' 'MANU' 'Temps' ltps1 'D'    dc    ;
  an1 = EVOL 'TURQ' 'MANU' 'EPXX'  lepu  'SMXX' lsig1 ;
  ansmzz = EVOL 'TURQ' 'MANU' 'Temps' ltps1 'SMZZ' lsig1 ;
  SI ITRAC1 ;
    DESS (ev0 ET an0) 'TITR' (CHAI str4 ' ' '[2D Axi.] Endommagement vs Temps') 'LEGE' 'SE' tdess ;
    DESS (ev1 ET an1) 'TITR' (CHAI str4 ' ' '[2D Axi.] Contrainte vs Deformation') 'LEGE' 'NE' tdess ;
    DESS ev2 'TITR' (CHAI str4 ' ' '[2D Axi.] Reaction vs Deplacement') 'LEGE' 'SE' tdess ;
  FINSI ;
* evolutions temporelles a tous les points de Gauss
  ng = 4 ;
  tleg = TABL ;
  tleg . 1 = MOT 'MARQ CARR NOLI' ;
  tleg . 2 = MOT 'MARQ LOSA NOLI' ;
  tleg . 3 = MOT 'MARQ ROND NOLI' ;
  tleg . 4 = MOT 'MARQ ETOI NOLI' ;
  tleg . 'TITRE' = TABL ;
*
  evend = VIDE 'EVOLUTIO' ;
  evsmzz = VIDE 'EVOLUTIO' ;
  evepzz = VIDE 'EVOLUTIO' ;
  REPE b ng ;
    tleg . 'TITRE' . &b = CHAI 'Cal_mat_' str3 '_Pt_Gauss' ' ' &b ;
    evend = evend ET (EVOL 'ROSE' 'TEMP' t 'VARIABLES_INTERNES' 'D' 1 1 &b) ;
    evsmzz = evsmzz ET (EVOL 'BLEU' 'TEMP' t 'CONTRAINTES' 'SMZZ' 1 1 &b) ;
    evepzz = evepzz ET (EVOL 'VERT' 'TEMP' t 'DEFORMATIONS' 'EPZZ' 1 1 &b) ;
  FIN b ;
*
  tleg . (ng+1) = MOT 'MARQ CROI' ;
  tleg . 'TITRE' . (ng+1) = MOT 'Reference' ;
*
  SI ITRAC1 ;
    DESS (evend ET an0) 'LEGE' 'SE' tleg 'TITR' (CHAI str4 ' ' '[2D Axi.] D vs Temps') ;
    DESS (evsmzz ET ansmzz) 'LEGE' 'NE' tleg 'TITR' (CHAI str4 ' ' '[2D Axi.] SMZZ vs Temps') ;
    DESS evepzz 'LEGE' 'NE' tleg 'TITR' (CHAI str4 ' ' '[2D Axi.] EPZZ vs Temps') ;
  FINSI ;
* deformee
  REPE b2 (DIME tt) ;
    tps1 = tt . (&b2 - 1) ;
*   trace de la deformee (pour controle visuel)
    def1 = DEFO mail (tu . (&b2 - 1)) 33. 'ROUG' ;
    tit1 = CHAI str4 ' ' '[2D Axi.] Deformee au temps' ' ' tps1 ;
    SI ITRAC1 ;
*      TRAC (def0 ET def1) 'TITR' tit1 ;
    FINSI ;
  FIN b2 ;
  SI ITRAC1 ;
    TRAC (def0 ET def1) 'TITR' tit1 ;
  FINSI ;
* comparaison Calcul / Reference sur la contrainte moyenne
  ler1 = ABS (((ENLE (ENLE lsig 1) lfiltr) - lsig1) / lsig1) ;
  eer1 = EVOL 'ORAN' 'MANU' 'Temps'  ltps1  '|Ec.rel|SMZZ' ler1 ;
  dtps = (EXTR ltps1 (DIME ltps1)) - (EXTR ltps1 1) ;
  mer1 = (INTG eer1) / dtps ;
  tit1 = CHAI str4 ' ' '[2D Axi.] Val.abs.Ec.rel.|(Calc-Ref)/Ref|, Moy.int.|Ecart|' ' ' mer1 ;
  SI ITRAC1 ;
    DESS eer1 'TITR' tit1 ;
  FINSI ;

** Sauvegarde des resultats de calcul Cast3M (reference pour MFront) au format .csv
  SI ISAUV ;
    str2 = '2D_axi_' ;
    fic1 = CHAI str1 str2 str3 '.csv' ;
    OPTI SORT fic1 ;
    SORT 'EXCE' ev1 'DIGI' 16 ;
  FINSI ;
  
  OPTI ECHO 0 ;
  SAUT 1 LIGN ;
  MESS ' ************ C A S  #4  :  2 D   A X I S Y M E T R I Q U E *************';
  SAUT 1 LIGN ;
  SI (mer1 < 1.E-6) ;
    MESS ' Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ' < 1.E-06 => OK' ;
    SAUT 1 LIGN ;
  SINO ;
    MESS ' Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ' > 1.E-06 => ECHEC' ;
    IERR = VRAI ;
    SAUT 1 LIGN ;
  FINS ;
  SAUT 1 LIGN ;
  OPTI ECHO 1 ;
FINSI ;





*********** C A S  #5  :  2 D   P O U T R E   A   F I B R E ************
SI ((EXIS lcas 5) ET (EGA &BCALC 1)) ;

** Options generales
  OPTI 'DIME' 2 'MODE' 'PLAN' 'ELEM' 'QUA4' ;

** Maillage
  p1s     = (-0.5 * epai) (-0.5 * epai) ;
  p2s     = (-0.5 * epai) ( 0.5 * epai) ;
  l12s    = DROI 1 p1s p2s ;
  mails   = l12s TRAN 1 (epai 0.) ;
  p1      = 0. 0. ;
  p2      = long 0. ;
  mail    = DROI 1 p1 p2 ;
  mes1    = MESU mails ;
cadr1 = COUL ((-0.1 0.) ET (1. 0.)) 'BLAN' ;
rep1 = @REPERE (-0.2 -0.2) (PROG 0.2 0.2) 'TURQ' ;
SI ITRAC1 ;
  TRAC 'QUAL' (rep1 ET (mails COUL 'VERT') ET (mail COUL 'ROUG')) 'TITR' (CHAI str4 ' ' '[2D PaF] Maillages de la Poutre et de sa Section') ;
FINSI ;

** Modele et caracteristiques materiau (&BCALC : 1->Cast3M, 2->MFront)
SI (EGA &BCALC 1) ;
  mos = MODE mails 'MECANIQUE' 'ELASTIQUE' 'PLASTIQUE' 'MAZARS' 'QUAS' ;
  mas = MATE mos 'YOUN' you 'NU' nu 'KTR0' epd0 'ACOM' ac 'BCOM' bc 'ATRA' at 'BTRA' bt 'BETA' beta
                 'ALPY' 1. ;
  mo = MODE mail 'MECANIQUE' 'ELASTIQUE' 'SECTION' 'PLASTIQUE' 'SECTION' 'TIMO' ;
  ma = MATE mo 'MODS' mos 'MATS' mas ;
SINO ;
* Implémentation MFront A FAIRE. En attendant on quitte la boucle...
  QUIT BCALC;
FINSI ;

** Blocages
  blgauche = BLOQ 'DEPL' 'ROTA' p1 ;
  bldroite = BLOQ 'UX' p2 ;
  bl       = blgauche ET bldroite ;
  bldroit2 = BLOQ 'ROTA' p2 ;
  bl       = bl ET bldroit2 ;

** Chargements
  didroite = DEPI bldroite umax ;
  ev1      = EVOL 'MANU' (PROG 0. 1.) (PROG 0. 1.) ;
  cha      = CHAR 'DIMP' didroite ev1 ;

** Tracé des blocages et du chargement
  re0 = @REPERE (-0.2 -0.2) (PROG 0.2 0.2) ;
  mor = MODE re0 'MECANIQUE' 'BARR' ;
  mar = MATE mor 'YOUN' 1. 'NU' 1. 'SECT' 1. ;
  rir = RIGI mor mar ;
  vf  = VECT didroite 'FLX' ' ' 'ROUG' 20. ;
  SI ITRAC1 ;
    TRAC ((RIGI mo ma) ET blgauche ET rir) vf 'NOLE' 'TITR' (CHAI str4 ' ' '[2D PaF] Blocages et chargement') ;
  FINSI ;

** Resolution
  t                               = TABL ;
  t . 'MODELE'                    = mo ;
  t . 'CARACTERISTIQUES'          = ma ;
  t . 'BLOCAGES_MECANIQUES'       = bl ;
  t . 'CHARGEMENT'                = cha ;
  t . 'TEMPS_CALCULES'            = PROG 0. 'PAS' 0.01 1. ;
  t . 'MES_SAUVEGARDES'           = TABL ;
  t . 'MES_SAUVEGARDES' . 'DEFTO' = VRAI ;
  PASAPAS t ;

** Post traitement
* courbes : endommagement moyen vs temps
*           contrainte moyenne vs deformation moyenne
*           force de reaction vs deplacement impose
* evolutions temporelles a tous les points de Gauss
  ng = 4 ;
  tleg = TABL ;
  tleg . 1 = MOT 'MARQ CARR NOLI' ;
  tleg . 2 = MOT 'MARQ LOSA NOLI' ;
  tleg . 3 = MOT 'MARQ ROND NOLI' ;
  tleg . 4 = MOT 'MARQ ETOI NOLI' ;
  tleg . 'TITRE' = TABL ;
  REPE b ng ;
    tleg . 'TITRE' . &b = CHAI 'Cal_mat_' str3 '_Pt_Gauss' ' ' &b ;
  FIN b ;
  tleg . (ng+1) = MOT 'MARQ CROI' ;
  tleg . 'TITRE' . (ng+1) = MOT 'Reference' ;
  tt = t . 'TEMPS' ;
  tv = t . 'VARIABLES_INTERNES' ;
  tu = t . 'DEPLACEMENTS' ;
  tr = t . 'REACTIONS' ;
  def0 = DEFO (cadr1 ET rep1 ET mail) (tu . 0) 0. ;
  ltps = PROG ;
  lend1 = PROG ;
  lend2 = PROG ;
  lend3 = PROG ;
  lend4 = PROG ;
  lsig1 = PROG ;
  lsig2 = PROG ;
  lsig3 = PROG ;
  lsig4 = PROG ;
  leps1 = PROG ;
  leps2 = PROG ;
  leps3 = PROG ;
  leps4 = PROG ;
  ldep = PROG ;
  lrea = PROG ;
  REPE b1 (DIME tt) ;
    tps1 = tt . (&b1 - 1) ;
    ltps = ltps ET tps1 ;
    SI (EGA &b1 1) ;
      lend1 = lend1 ET 0. ;
      lend2 = lend2 ET 0. ;
      lend3 = lend3 ET 0. ;
      lend4 = lend4 ET 0. ;
      lsig1 = lsig1 ET 0. ;
      lsig2 = lsig2 ET 0. ;
      lsig3 = lsig3 ET 0. ;
      lsig4 = lsig4 ET 0. ;
      leps1 = leps1 ET 0. ;
      leps2 = leps2 ET 0. ;
      leps3 = leps3 ET 0. ;
      leps4 = leps4 ET 0. ;
    SINON ;
      lend1 = lend1 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VHNI' 1 1 1) ;
      lend2 = lend2 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VHNI' 1 1 2) ;
      lend3 = lend3 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VHNI' 1 1 3) ;
      lend4 = lend4 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VHNI' 1 1 4) ;
      lsig1 = lsig1 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VONS' 1 1 1) 'SMXX' 1 1 1) ;
      lsig2 = lsig2 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VONS' 1 1 1) 'SMXX' 1 1 2) ;
      lsig3 = lsig3 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VONS' 1 1 1) 'SMXX' 1 1 3) ;
      lsig4 = lsig4 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VONS' 1 1 1) 'SMXX' 1 1 4) ;
      leps1 = leps1 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VISX' 1 1 1) ;
      leps2 = leps2 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VISX' 1 1 2) ;
      leps3 = leps3 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VISX' 1 1 3) ;
      leps4 = leps4 ET (EXTR (EXTR (tv . (&b1 - 1)) 'VAIS' 1 1 1) 'VISX' 1 1 4) ;
    FINSI ;
    u1   = EXTR (tu . (&b1 - 1)) 'UX' p2 ;
    rea1 = EXTR (tr . (&b1 - 1)) 'FX' p1 ;
    ldep = ldep ET u1 ;
    lrea = lrea ET rea1 ;
  FIN b1 ;
  ev0 = (EVOL 'ROUG' 'MANU' 'Temps' ltps 'D' lend1) ET (EVOL 'ROUG' 'MANU' 'Temps' ltps 'D' lend2) ET
        (EVOL 'ROUG' 'MANU' 'Temps' ltps 'D' lend3) ET (EVOL 'ROUG' 'MANU' 'Temps' ltps 'D' lend4) ;
  ev1 = (EVOL 'ROUG' 'MANU' 'EPXX'  leps1 'SMXX' lsig1) ET (EVOL 'ROUG' 'MANU' 'EPXX' leps2 'SMXX' lsig2) ET
        (EVOL 'ROUG' 'MANU' 'EPXX'  leps3 'SMXX' lsig3) ET (EVOL 'ROUG' 'MANU' 'EPXX' leps4 'SMXX' lsig4) ;
  ev2 = EVOL 'ROUG' 'MANU' 'UX' ldep 'FX' lrea ;
  ltps1 =  ENLE ltps 1 ;
  nb1 = DIME ltps1 ;
  un1 = PROG nb1 * 1.D0 ;
  lepd0 = PROG nb1 * epd0 ;
  lac1 = PROG nb1 * ac ;
  lepu = ltps1 * umax / long ;
  lepe = (2.**0.5) * nu * (ABS lepu) ;
  dc = un1 - ((epd0*(un1-lac1)/lepe) + (ac*(EXP(bc*(lepd0-lepe))))) ;
  lfiltr = LECT ;
  REPE b2 ;
    dc_i = EXTR dc &b2 ;
    SI (EGA &b2 1) ;
      dc_p = 1. ;
    FINSI ;
    ddc = dc_i - dc_p ;
    SI ((DDC < 0.) OU (dc_i < 0.) OU (dc_i > 1.)) ;
      lfiltr = lfiltr ET (LECT &b2) ;
    FINSI ;
    dc_p = dc_i ;
    SI (EGA &b2 (DIME dc)) ;
      QUIT b2 ;
    FINSI ;
  FIN b2 ;
  ltps1 = ENLE ltps1 lfiltr ;
  un1   = ENLE un1   lfiltr ;
  lepu  = ENLE lepu  lfiltr ;
  dc    = ENLE dc    lfiltr ;
  lsan1 = (un1 - dc) * you * lepu ;
  an0 = EVOL 'TURQ' 'MANU' 'Temps' ltps1 'D'    dc    ;
  an1 = EVOL 'TURQ' 'MANU' 'EPXX'  lepu  'SMXX' lsan1 ;
  SI ITRAC1 ;
    DESS (ev0 ET an0) 'LEGE' 'SE' tleg 'TITR' (CHAI str4 ' ' '[2D PaF] Endommagement vs Temps') ;
    DESS (ev1 ET an1) 'LEGE' 'NE' tleg 'TITR' (CHAI str4 ' ' '[2D PaF] Contrainte vs Deformation') ;
    DESS ev2 'TITR' (CHAI str4 ' ' '[2D PaF] Reaction vs Deplacement') 'LEGE' 'SE' tdess ;
  FINSI ;
* deformee
  REPE b2 (DIME tt) ;
    tps1 = tt . (&b2 - 1) ;
*   trace de la deformee (pour controle visuel)
    def1 = DEFO mail (tu . (&b2 - 1)) 33. 'ROUG' ;
    tit1 = CHAI str4 ' ' '[2D PaF] Deformee au temps' ' ' tps1 ;
    SI ITRAC1 ;
*      TRAC (def0 ET def1) 'TITR' tit1 ;
    FINSI ;
  FIN b2 ;
  SI ITRAC1 ;
    TRAC (def0 ET def1) 'TITR' tit1 ;
  FINSI ;
* comparaison Calcul / Reference sur la contrainte moyenne
  lsig = (lsig1 + lsig2 + lsig3 + lsig4) / 4. ;
  ler1 = ABS (((ENLE (ENLE lsig 1) lfiltr) - lsan1) / lsan1) ;
  eer1 = EVOL 'ORAN' 'MANU' 'Temps'  ltps1  '|Ec.rel|SMXX' ler1 ;
  dtps = (EXTR ltps1 (DIME ltps1)) - (EXTR ltps1 1) ;
  mer1 = (INTG eer1) / dtps ;
  tit1 = CHAI str4 ' ' '[2D PaF] Val.abs.Ec.rel.|(Calc-Ref)/Ref|, Moy.int.|Ecart|' ' ' mer1 ;
  SI ITRAC1 ;
    DESS eer1 'TITR' tit1 ;
  FINSI ;
  OPTI ECHO 0 ;
  SAUT 1 LIGN ;
  MESS ' *********** C A S  #5  :  2 D   P O U T R E   A   F I B R E ************';
  SAUT 1 LIGN ;
  SI (mer1 < 1.E-6) ;
    MESS ' Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ' < 1.E-06 => OK' ;
    SAUT 1 LIGN ;
  SINO ;
    MESS ' Moy.int. Val.abs. Ecart relatif entre Calcul et Reference' ' ' mer1 ' > 1.E-06 => ECHEC' ;
    IERR = VRAI ;
    SAUT 1 LIGN ;
  FINS ;
  SAUT 1 LIGN ;
  OPTI ECHO 1 ;
FINSI ;

** Sortie sur Erreur 1152 en cas de non conformite
SI IERR ;
  ERRE 1152 ;
  ERRE 5 ;
FINSI ;


FIN BCALC ;


FIN ;