// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL  : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.

#include "UmatAbaqus.h" 

#include "NevezTenseur.h"
#include "NevezTenseurQ.h"
#include "ParaGlob.h"
#include "TenseurQ-3.h"
#include "TenseurQ-2.h"
#include "Tenseur3.h"
#include "Tenseur2.h"
#include "ConstMath.h"

// --- fichier pour le tube ------
#ifdef SYSTEM_MAC_OS_CARBON
  #include <types.h>
  #include <stat.h>
#else  
  #include <sys/types.h>
  #include <sys/stat.h>
#endif

#include <stdio.h>
#include <fcntl.h> /* Pour O_WRONLY, etc */

#include <unistd.h>    // pour les ordres read write close


//------------------------- variables statiques --------------------------- 

const UmatAbaqus::ChangementIndex  UmatAbaqus::cdex;

// init de la variable statique cdex: cas 3D
// fonctionne aussi avec le cas 3D axi, c'est uniquement le terme (1,2) (= (2,1)) qui est
// non nul
UmatAbaqus::ChangementIndex::ChangementIndex() :
  idx_i(6),idx_j(6),odVect(3)
  { idx_i(1)=1;idx_i(2)=2;idx_i(3)=3;idx_i(4)=1;idx_i(5)=1;idx_i(6)=2;
    idx_j(1)=1;idx_j(2)=2;idx_j(3)=3;idx_j(4)=2;idx_j(5)=3;idx_j(6)=3;
    odVect(1,1)=1;odVect(2,1)=4;odVect(3,1)=5;
    odVect(1,2)=4;odVect(2,2)=2;odVect(3,2)=6;
    odVect(1,3)=5;odVect(2,3)=6;odVect(3,3)=3;
  }; 
                           
const UmatAbaqus::ChangementIndexCP  UmatAbaqus::cdexCP;

// init de la variable statique cdexCP : cas 3D CP
UmatAbaqus::ChangementIndexCP::ChangementIndexCP() :
// 11 -> 1  22 -> 2  12 -> 3
  idx_i(3),idx_j(3),odVect(2)
  { idx_i(1)=1;idx_i(2)=2;idx_i(3)=1;
    idx_j(1)=1;idx_j(2)=2;idx_j(3)=2;
    odVect(1,1)=1;odVect(1,2)=3;odVect(2,2)=2;
    odVect(2,1)=3; // stockage symétrique
  };

//---------------------- fin variables statiques --------------------------- 

    // par défaut
    UmatAbaqus::UmatAbaqus(int dimension_espace, int nb_vecteur):
    // si dimension = 2, on ne peut accepter que nb_vecteur = 2
    // donc on prend MiN(dimension_espace,nb_vecteur)
       STRESS(NULL),t_sigma(NevezTenseurHH(MiN(dimension_espace,nb_vecteur)))
       ,STATEV(NULL),v_statev()
      ,DDSDDE(NULL)
      ,d_sigma_deps(NevezTenseurHHHH(MiN(dimension_espace,nb_vecteur)*10+MiN(dimension_espace,nb_vecteur)))
      ,SSE(NULL),SPD(NULL),SCD(NULL)
      ,energie_elastique(0.),dissipation_plastique(0.),dissipation_visqueuse(0.)
      ,RPL(NULL),ener_therm_vol_par_utemps(0.)
      ,DDSDDT(NULL)
      ,d_sigma_d_tempera(NevezTenseurHH(MiN(dimension_espace,nb_vecteur)))
      ,DRPLDE(NULL)
      ,d_ener_therm_vol_deps(NevezTenseurHH(MiN(dimension_espace,nb_vecteur)))
      ,DRPLDT(NULL),d_ener_therm_dtemper(0.)
      ,PNEWDT(NULL),pnewdt(0.)
      ,STRAN(NULL)
      ,eps_meca(NevezTenseurBB(MiN(dimension_espace,nb_vecteur)))
      ,DSTRAN(NULL)
      ,delta_eps_meca(NevezTenseurBB(MiN(dimension_espace,nb_vecteur)))
      ,TIME(NULL),temps_t(0.),temps_tdt(0.)
      ,DTIME(NULL),delta_t(0.)
      ,TemP(NULL),temper_t(0.)
      ,DTEMP(NULL),delta_temper(0.)
      ,PREDEF(NULL),ddlNoeA_ptinteg()
      ,DPRED(NULL),delta_ddlNoeA_ptinteg()
      ,CMNAME(NULL),nom_materiau()
      ,NDI(NULL),dim_tens(0)
      ,NSHR(NULL),nb_tau_ij(0)
      ,NTENS(NULL),nb_ij_tenseur(0)
      ,NSTATV(NULL),nb_variables_etat(0)
      ,PROPS(NULL),proprietes_materiau()
      ,NPROPS(NULL),nb_proprietes_materiau(0)
      ,COORDS(NULL),coor_pt(dimension_espace)
      ,DROT(NULL),mat_corota(3,3) // je crois que l'on est obligé d'être toujours en 3x3
                                   // pour abaqus
      ,CELENT(NULL),longueur_characteristique(-1)
      ,DFGRD0(NULL)
      ,giB_t(dimension_espace,MiN(dimension_espace,nb_vecteur))
      ,DFGRD1(NULL)
      ,giB_tdt(dimension_espace,MiN(dimension_espace,nb_vecteur))
      ,NOEL(NULL),nb_elem(0)
      ,NPT(NULL),nb_pt_integ(0)
      ,LAYER(NULL),nb_du_plis(0)
      ,KSPT(NULL),nb_pt_dans_le_plis(0)
      ,KSTEP(NULL),nb_step(0)
      ,KINC(NULL),nb_increment(0)
      ,passage(0),envoi("Umat_envoi_Hz"),reception("Umat_reception_Hz")
      ,N_t(dimension_espace),N_tdt(dimension_espace)
      	{
       
       
       
         // on définit les grandeurs pointées au bon endroit dans la chaine de 
      	  // caractères
          // on considère que les int sont codés sur 4 octets, les doubles sur 8 
          // On se met dans le cas en 3D (taille maxi) d'où l'ensemble des infos
          // STRESS 6*8=48,  STATEV rien on n'en tiend pas compte
          // DDSDDE 36*8= 288, SSE+SPD+SCD 3*8=24
          // RPL 8,  DDSDDT 6*8=48,  DRPLDE 6*8=48,  DRPLDT 8,  PNEWDT 8,
          // STRAN 6*8=48, DSTRAN 6*8=48
          // TIME 2*8=16, DTIME 8, TemP 8,  DTEMP 8,
          // PREDEF rien on n'en tiens pas compte
          // DPRED rien on n'en tiens pas compte
          // CMNAME 24 (par défaut sur 24 caractères: dans abaqus c'est sur 80)
          // NDI+NSHR+NTENS+NSTATV 4*4=16,
          // PROPS rien on n'en tiens pas compte
          // NPROPS rien on n'en tiens pas compte
          // COORDS 3*8=24, DROT 9*8=72,
          // CELENT 8,
          // DFGRD0 9*8=72,  DFGRD1 9*8=72,
          // NOEL+NPT+LAYER+KSPT+KSTEP+KINC 6*4=24
          // total 952 octets
          // --- en sortie uniquement
          DDSDDT  =  &t_car_x_n.x[0]; // 6
          DRPLDE  =  &t_car_x_n.x[6]; // 6
          DDSDDE =   &t_car_x_n.x[12]; // 36
          // --- en entrée - sortie
          RPL  =     &t_car_x_n.x[48]; // 1
          STRESS =   &t_car_x_n.x[49]; // 6
          SSE =      &t_car_x_n.x[55]; // 1 
          SPD =      &t_car_x_n.x[56]; // 1
          SCD =      &t_car_x_n.x[57]; // 1
          DRPLDT  =  &t_car_x_n.x[58]; // 1
          PNEWDT  =  &t_car_x_n.x[59]; // 1
          // --- en entrée seulement 
          STRAN  =   &t_car_x_n.x[60]; // 6
          DSTRAN  =  &t_car_x_n.x[66]; // 6
          TIME  =    &t_car_x_n.x[72]; // 2
          DTIME =    &t_car_x_n.x[74]; // 1
          TemP  =    &t_car_x_n.x[75]; // 1
          DTEMP  =   &t_car_x_n.x[76]; // 1
          COORDS  =  &t_car_x_n.x[77]; // 3
          DROT  =    &t_car_x_n.x[80]; // 9
          CELENT  =  &t_car_x_n.x[89]; // 1
          DFGRD0  =  &t_car_x_n.x[90]; // 9
          DFGRD1  =  &t_car_x_n.x[99]; // 9

          NDI  =     &t_car_x_n.n[216]; // 1
          NSHR  =    &t_car_x_n.n[217]; // 1
          NTENS =    &t_car_x_n.n[218]; // 1
          NSTATV  =  &t_car_x_n.n[219]; // 1
          NOEL  =    &t_car_x_n.n[220]; // 1
		        NPT =      &t_car_x_n.n[221];  // 1
          LAYER =    &t_car_x_n.n[222]; // 1
          KSPT =     &t_car_x_n.n[223]; // 1
          KSTEP =    &t_car_x_n.n[224]; // 1
          KINC =     &t_car_x_n.n[225]; // 1

          CMNAME  =  &t_car_x_n.tampon[904];
          // le nom
          nom_materiau=CMNAME;
      	};

    // constructeur de copie
    UmatAbaqus::UmatAbaqus(const UmatAbaqus& a) :       
      STRESS(NULL),t_sigma(NULL)
      ,STATEV(NULL),v_statev(a.v_statev)
      ,DDSDDE(NULL),d_sigma_deps(NULL)
      ,SSE(NULL),SPD(NULL),SCD(NULL)
      ,energie_elastique(a.energie_elastique)
      ,dissipation_plastique(a.dissipation_plastique)
      ,dissipation_visqueuse(a.dissipation_visqueuse)
      ,RPL(NULL),ener_therm_vol_par_utemps(a.ener_therm_vol_par_utemps)
      ,DDSDDT(NULL),d_sigma_d_tempera(NULL)
      ,DRPLDE(NULL),d_ener_therm_vol_deps(NULL)
      ,DRPLDT(NULL),d_ener_therm_dtemper(a.d_ener_therm_dtemper)
      ,PNEWDT(NULL),pnewdt(a.pnewdt)
      ,STRAN(NULL),eps_meca(NULL)
      ,DSTRAN(NULL),delta_eps_meca(NULL)
      ,TIME(NULL),temps_t(a.temps_t),temps_tdt(a.temps_tdt)
      ,DTIME(NULL),delta_t(a.delta_t)
      ,TemP(NULL),temper_t(a.temper_t)
      ,DTEMP(NULL),delta_temper(a.delta_temper)
      ,PREDEF(NULL),ddlNoeA_ptinteg(a.ddlNoeA_ptinteg)
      ,DPRED(NULL),delta_ddlNoeA_ptinteg(a.delta_ddlNoeA_ptinteg)
      ,CMNAME(NULL),nom_materiau(a.nom_materiau)
      ,NDI(NULL),dim_tens(a.dim_tens)
      ,NSHR(NULL),nb_tau_ij(a.nb_tau_ij)
      ,NTENS(NULL),nb_ij_tenseur(a.nb_ij_tenseur)
      ,NSTATV(NULL),nb_variables_etat(a.nb_variables_etat)
      ,PROPS(NULL),proprietes_materiau(a.proprietes_materiau)
      ,NPROPS(NULL),nb_proprietes_materiau(a.nb_proprietes_materiau)
      ,COORDS(NULL),coor_pt(a.coor_pt)
      ,DROT(NULL),mat_corota(a.mat_corota)
      ,CELENT(NULL),longueur_characteristique(a.longueur_characteristique)
      ,DFGRD0(NULL),giB_t(a.giB_t)
      ,DFGRD1(NULL),giB_tdt(a.giB_tdt)
      ,NOEL(NULL),nb_elem(a.nb_elem)
      ,NPT(NULL),nb_pt_integ(a.nb_pt_integ)
      ,LAYER(NULL),nb_du_plis(a.nb_du_plis)
      ,KSPT(NULL),nb_pt_dans_le_plis(a.nb_pt_dans_le_plis)
      ,KSTEP(NULL),nb_step(a.nb_step)
      ,KINC(NULL),nb_increment(a.nb_increment)
      ,passage(a.passage),envoi(a.envoi),reception(a.reception)
      ,N_t(a.N_t),N_tdt(a.N_tdt)
      	{ // on définit les grandeurs pointées au bon endroit dans la chaine de
      	  // caractères
          // on considère que les int sont codés sur 4 octets, les doubles sur 8 
          // On se met dans le cas en 3D (taille maxi) d'où l'ensemble des infos
          // STRESS 6*8=48,  STATEV rien on n'en tiend pas compte
          // DDSDDE 36*8= 288, SSE+SPD+SCD 3*8=24
          // RPL 8,  DDSDDT 6*8=48,  DRPLDE 6*8=48,  DRPLDT 8,  PNEWDT 8,
          // STRAN 6*8=48, DSTRAN 6*8=48
          // TIME 2*8=16, DTIME 8, TemP 8,  DTEMP 8,
          // PREDEF rien on n'en tiens pas compte
          // DPRED rien on n'en tiens pas compte
          // CMNAME 20 (par défaut sur 20 caractères: dans abaqus c'est sur 80)
          // NDI+NSHR+NTENS+NSTATV 4*4=16,
          // PROPS rien on n'en tiens pas compte
          // NPROPS rien on n'en tiens pas compte
          // COORDS 3*8=24, DROT 9*8=72,
          // CELENT 8,
          // DFGRD0 9*8=72,  DFGRD1 9*8=72,
          // NOEL+NPT+LAYER+KSPT+KSTEP+KINC 6*4=24
          // total 924 octets
		        int* pt1 = &(t_car_x_n.n[0]); const int *pt2= &(a.t_car_x_n.n[0]);
		        for (int i=0;i<232;i++,pt1++,pt2++) *pt1=*pt2;
		    
          // --- en sortie uniquement
          DDSDDT  =  &t_car_x_n.x[0]; // 6
          DRPLDE  =  &t_car_x_n.x[6]; // 6
          DDSDDE =   &t_car_x_n.x[12]; // 36
          // --- en entrée - sortie
          RPL  =     &t_car_x_n.x[48]; // 1
          STRESS =   &t_car_x_n.x[49]; // 6
          SSE =      &t_car_x_n.x[55]; // 1 
          SPD =      &t_car_x_n.x[56]; // 1
          SCD =      &t_car_x_n.x[57]; // 1
          DRPLDT  =  &t_car_x_n.x[58]; // 1
          PNEWDT  =  &t_car_x_n.x[59]; // 1
          // --- en entrée seulement 
          STRAN  =   &t_car_x_n.x[60]; // 6
          DSTRAN  =  &t_car_x_n.x[66]; // 6
          TIME  =    &t_car_x_n.x[72]; // 2
          DTIME =    &t_car_x_n.x[74]; // 1
          TemP  =    &t_car_x_n.x[75]; // 1
          DTEMP  =   &t_car_x_n.x[76]; // 1
          COORDS  =  &t_car_x_n.x[77]; // 3
          DROT  =    &t_car_x_n.x[80]; // 9
          CELENT  =  &t_car_x_n.x[89]; // 1
          DFGRD0  =  &t_car_x_n.x[90]; // 9
          DFGRD1  =  &t_car_x_n.x[99]; // 9

          NDI  =     &t_car_x_n.n[216]; // 1
          NSHR  =    &t_car_x_n.n[217]; // 1
          NTENS =    &t_car_x_n.n[218]; // 1
          NSTATV  =  &t_car_x_n.n[219]; // 1
          NOEL  =    &t_car_x_n.n[220]; // 1
		        NPT =      &t_car_x_n.n[221];  // 1
          LAYER =    &t_car_x_n.n[222]; // 1
          KSPT =     &t_car_x_n.n[223]; // 1
          KSTEP =    &t_car_x_n.n[224]; // 1
          KINC =     &t_car_x_n.n[225]; // 1

          CMNAME  =  &t_car_x_n.tampon[904];
          // le nom
          nom_materiau=CMNAME;
      	  // cas des grandeurs évoluées
      	  if (a.t_sigma != NULL) 
      	     t_sigma = NevezTenseurHH(*(a.t_sigma));
      	  if (a.d_sigma_deps != NULL )
      	     d_sigma_deps = NevezTenseurHHHH(*(a.d_sigma_deps));
          if (a.d_sigma_d_tempera != NULL) 
             d_sigma_d_tempera = NevezTenseurHH(*(a.d_sigma_d_tempera));
          if (a.d_ener_therm_vol_deps != NULL) 
             d_ener_therm_vol_deps =NevezTenseurHH(*(a.d_ener_therm_vol_deps));
          if (a.eps_meca != NULL) 
             eps_meca = NevezTenseurBB(*(a.eps_meca)); 
          if (a.delta_eps_meca != NULL) 
             delta_eps_meca = NevezTenseurBB(*(a.delta_eps_meca));
      	};

    // DESTRUCTEUR :
    UmatAbaqus::~UmatAbaqus()
    	{ // cas des grandeurs évoluées
      	  if (t_sigma != NULL) 
      	       delete t_sigma;
      	  if (d_sigma_deps != NULL )
      	       delete d_sigma_deps;
          if (d_sigma_d_tempera != NULL) 
               delete d_sigma_d_tempera;
          if (d_ener_therm_vol_deps != NULL) 
               delete d_ener_therm_vol_deps;
          if (eps_meca != NULL) 
               delete eps_meca; 
          if (delta_eps_meca != NULL) 
              delete delta_eps_meca;
    	};

// Surcharge de l'operateur = 
UmatAbaqus& UmatAbaqus::operator= (UmatAbaqus& a)
{   t_sigma = a.t_sigma;
    v_statev = a.v_statev;
    d_sigma_deps = a.d_sigma_deps;
    energie_elastique = a.energie_elastique;
    dissipation_plastique = a.dissipation_plastique;
    dissipation_visqueuse = a.dissipation_visqueuse;
    ener_therm_vol_par_utemps = a.ener_therm_vol_par_utemps;
    d_sigma_d_tempera = a.d_sigma_d_tempera;
    d_ener_therm_vol_deps = a.d_ener_therm_vol_deps;
    d_ener_therm_dtemper = a.d_ener_therm_dtemper;
    pnewdt = a.pnewdt;
    eps_meca = a.eps_meca;
    delta_eps_meca = a.delta_eps_meca;
    temps_t = a.temps_t; temps_tdt = a.temps_tdt;
    delta_t = a.delta_t;
    temper_t = a.temper_t;
    delta_temper = a.delta_temper;
    ddlNoeA_ptinteg = a.ddlNoeA_ptinteg;
    delta_ddlNoeA_ptinteg = a.delta_ddlNoeA_ptinteg;
    nom_materiau = a.nom_materiau;
    dim_tens = a.dim_tens;
    nb_tau_ij = a.nb_tau_ij;
    nb_ij_tenseur = a.nb_ij_tenseur;
    nb_variables_etat = a.nb_variables_etat;
    proprietes_materiau = a.proprietes_materiau;
    nb_proprietes_materiau = a.nb_proprietes_materiau;
    coor_pt = a.coor_pt;
    mat_corota = a.mat_corota;
    longueur_characteristique = a.longueur_characteristique;
    giB_t = a.giB_t;
    giB_tdt = a.giB_tdt;
    nb_elem = a.nb_elem;
    nb_pt_integ = a.nb_pt_integ;
    nb_du_plis = a.nb_du_plis;
    nb_pt_dans_le_plis = a.nb_pt_dans_le_plis;
    nb_step = a.nb_step;
    nb_increment = a.nb_increment;
    N_t = a.N_t; N_tdt = a.N_tdt;
    return *this;		
	};    
    
// lecture des données Umat
void UmatAbaqus::LectureDonneesUmat(bool utilisation_umat_interne,const int niveau)
{ //  Creation d'un processus de reception des données
  // en fait l'objectif est de réceptionner que les variables d'entrées ou 
  // d'entrée sortie 
  if (niveau > 5)
       cout << "\n hzpp:  ouverture du tube en lecture, lecture sur le tube" << endl; 
  if (!utilisation_umat_interne)
    {// dans le cas où l'état actuel de l'incrément est 0,1 ou 2, on initialise le champ de caractères
     if (nb_increment<=2)
         for (int i=0;i<24;i++) CMNAME[i]='\0';
     char* tampon_envoi =  &(t_car_x_n.tampon[384]);  
     int tub = open(reception.c_str(),O_RDONLY);  
     read (tub,tampon_envoi,544);
     close (tub); // fermeture du tampon
     // transformation de l'information 
     Copie_base_vers_evolue(niveau);
     if (niveau > 5)
  	   { cout << "\n *********** grandeurs lues ************ ";
        cout << "\n sigma_t "; for (int i=0;i<6;i++) cout << "("<<i<<")= " << STRESS[i];
        cout << "\n def "; for (int i=0;i<6;i++) cout << "("<<i<<")= " << STRAN[i];
        cout << "\n deltatdef "; for (int i=0;i<6;i++) cout << "("<<i<<")= " << DSTRAN[i];
        cout << "\n rotation "; for (int i=0;i<9;i++) cout << "("<<i<<")= " << DROT[i];
        cout << "\n gradiant_t "; for (int i=0;i<9;i++) cout << "("<<i<<")= " << DFGRD0[i];
        cout << "\n gradiant_tdt "; for (int i=0;i<9;i++) cout << "("<<i<<")= " << DFGRD1[i];
        cout << "\n temps_t= " << temps_t << " temps_tdt= " <<  temps_tdt 
             << " delta_t= " << delta_t;
        cout << "\n temper_t= " <<  temper_t << " delta_temper= " <<  delta_temper;
        cout << "\n nom_materiau= " << nom_materiau; // nom du matériau
        cout << "\n dim_tens= " << dim_tens << " nb_tau_ij= " << nb_tau_ij
             << " nb_ij_tenseur= " <<  nb_ij_tenseur; 
        cout << "\n nb_variables_etat= " << nb_variables_etat; 
        cout << "\n Coordonnee= " <<  coor_pt; // coordonnee du point
        cout << "\n mat_corota= " << mat_corota; // matrice de rotation du corotationnelle
        cout << "\n longueur_characteristique= " << longueur_characteristique;
        cout << "\n BaseB giB_t= " << giB_t; // base naturelle à t
        cout << "\n BaseB giB_tdt = " << giB_tdt; // base natuelle à tdt
        if (giB_t.NbVecteur()==2)
          cout <<"\n N_t= "<<N_t <<", N_tdt= "<<N_tdt;
        cout << "\n nb_elem= " << nb_elem ; // numéro de l'élément
        cout << " nb_pt_integ= " << nb_pt_integ; // numéro du point d'intégration
        cout << " nb_du_plis= " << nb_du_plis; // numéro du plis pour les multicouches 
        cout << " nb_pt_dans_le_plis= " << nb_pt_dans_le_plis; // numéro du po     cout << "dans le plis
        cout << " nb_step= " << nb_step; // numéro du step
        cout << " nb_increment= " << nb_increment; // numéro de l'incrément
        cout << "\n  ******** fin grandeurs lues ************** " << endl;
     	};

     #ifdef SYSTEM_MAC_OS_CARBON
      // def de valeur par défaut pour tester
      STRESS[0]=0.;STRESS[1]=0.;STRESS[2]=0.;STRESS[3]=0.;STRESS[4]=0.;STRESS[5]=0.;
      *SSE = 0.;*SPD = 0.;*SCD = 0.;*RPL = 0.;
      STRAN[0] = 0.01; STRAN[1] = 0.; STRAN[2] = 0.; STRAN[3] = 0.; STRAN[4] = 0.; STRAN[5] = 0.;
      DSTRAN[0] = 0.01; DSTRAN[1] = 0.; DSTRAN[2] = 0.; DSTRAN[3] = 0.; DSTRAN[4] = 0.; DSTRAN[5] = 0.;
      TIME[0] = 0.; TIME[1] = 1.;*DTIME = 1.;
      CMNAME="acier";//CMNAME[5]='\0';
      *NDI=3;*NSHR = 3;*NTENS=6;*NSTATV=0;
      COORDS[0]=1.;COORDS[1]=0.;COORDS[2]=0.;
      DROT[0]=1.;DROT[1]=0.;DROT[2]=0.;DROT[3]=0.;DROT[4]=1.;DROT[5]=0.;DROT[6]=0.;DROT[7]=0.;DROT[8]=1.;
      *CELENT=0.;
      DFGRD0[0]=1.;DFGRD0[1]=0.;DFGRD0[2]=0.;DFGRD0[3]=0.;DFGRD0[4]=1.;
      DFGRD0[5]=0.;DFGRD0[6]=0.;DFGRD0[7]=0.;DFGRD0[8]=1.;
      DFGRD1[0]=1.;DFGRD1[1]=0.;DFGRD1[2]=0.;DFGRD1[3]=0.;DFGRD1[4]=1.;
      DFGRD1[5]=0.;DFGRD1[6]=0.;DFGRD1[7]=0.;DFGRD1[8]=1.;
      *NOEL=1;*NPT=1;*KSTEP=1;*KINC=1;
     #endif
    };
  if (niveau > 5)
       cout << "\n HZpp: fermeture du tube en lecture" << endl;
//  // transformation de l'information 
//  Copie_base_vers_evolue();
  // à chaque lecture on incrémente un compteur
  passage++; // pour des comptes éventuelles
};

// écriture des données Umat
void UmatAbaqus::EcritureDonneesUmat(bool utilisation_umat_interne,const int niveau)
{ //   Creation d'un processus d'ecriture des variables 
  // en fait l'objectif est d'ecrire que les variables qui évoluent 
  if (niveau > 5)
       cout << "\n hzpp: ouverture du tube en ecriture, ecriture  sur le tube " << endl; 
  if (!utilisation_umat_interne)
    { // transformation de l'information 
      Copie_evolue_vers_base(niveau);
      // ouverture du tube nomme en ecriture 
      int tub = open(envoi.c_str(),O_WRONLY);  
      write (tub,t_car_x_n.tampon,480);    
      close (tub); // fermeture du tampon
      #ifdef SYSTEM_MAC_OS_CARBON
      if (niveau > 5)
       {// on ecrit les résultats à l'écran
        if ((*NDI == 3) && (*NSHR == 3))
         { // cas classique 3D
           Tenseur3HH&  sigma =  *((Tenseur3HH*) t_sigma);
           cout << "\n ************* contraintes calculees: *********\n" << sigma;
           Tenseur3HHHH& dd_sigma_deps =  *((Tenseur3HHHH*) d_sigma_deps);
           cout << "\n ************* matrice tangente ***************\n" << dd_sigma_deps;
  //      Tenseur3HH&  d_sigma_d_tempe =  *((Tenseur3HH*) d_sigma_d_tempera);
  //      Tenseur3HH&  d_ener_therm_vol_dep =  *((Tenseur3HH*) d_ener_therm_vol_deps);
         }
        else if ((*NDI == 2) && (*NSHR == 1))
         { // cas contraintes planes
           Tenseur2HH&  sigma =  *((Tenseur2HH*) t_sigma);
           cout << "\n ************* contraintes calculees: *********\n" << sigma;
           Tenseur2HHHH& dd_sigma_deps =  *((Tenseur2HHHH*) d_sigma_deps);
           cout << "\n ************* matrice tangente ***************\n" << dd_sigma_deps;
  //      Tenseur3HH&  d_sigma_d_tempe =  *((Tenseur3HH*) d_sigma_d_tempera);
  //      Tenseur3HH&  d_ener_therm_vol_dep =  *((Tenseur3HH*) d_ener_therm_vol_deps);
           cout <<"\n N_t= "<<N_t <<", N_tdt= "<<N_tdt;
         }
        else
         { cout << "\n erreur , pour l'instant seule la dimension 3 est en place"
                << "\n il est vivement conseille de ce plaindre pour remedier a ce probleme "
                << "\n void UmatAbaqus::EcritureDonneesUmat()" << endl ;
           Sortie(1);     
          };
        cout << "\n ************* les energies *********" 
             << "\n energie_elastique: " << energie_elastique 
             << "  dissipation_plastique: " << dissipation_plastique
             << "  dissipation_visqueuse: " << dissipation_visqueuse;
        cout << "\n rentrer un caractere pour continuer ? ";
        string toto; toto= lect_chaine();       
//   *PNEWDT = pnewdt;
//   *RPL =  ener_therm_vol_par_utemps; 
//   *DRPLDT = d_ener_therm_dtemper; 
       };
     #endif
    };
  if (ParaGlob::NiveauImpression() > 7)
       cout << "\n HZpp: fermeture du tube en ecriture" << endl;
 };

// <<<<<< 2) <<<<<< contexte: utilisation par herezh d'une umat extérieure
// écriture des données pour l'umat 
void UmatAbaqus::EcritureDonneesPourUmat(bool utilisation_umat_interne,const int niveau)
{ //   Creation d'un processus d'ecriture des variables 
  //   en fait l'objectif est d'ecrire que les variables nécessaires pour le calcul 
  if (niveau > 5)
       cout << "\n hzpp-utilisation d'umat: ouverture du tube en ecriture, ecriture  sur le tube " 
            << endl; 
  if (!utilisation_umat_interne)
    {//transformation de l'information 
     Copie_evolue_vers_base_PourUmat(niveau);
     // ouverture du tube nomme en ecriture 
     char* tampon_envoi =  &(t_car_x_n.tampon[384]);  
     int tub = open(reception.c_str(),O_WRONLY);
     if (tub < 0)
      {cout << "\n *** erreur en ouverture du pipe " << reception.c_str()
            << " peut-etre n'existe t'il pas ? en tout cas il ne veut pas s'ouvrir !! ";
       Sortie(1);
      };
     
     write (tub,tampon_envoi,544);    
     if (niveau > 5)
  	   { cout << "\n *********** donnees envoyees ************ ";
        cout << "\n sigma_t "; for (int i=0;i<6;i++) cout << "("<<i<<")= " << STRESS[i];
        cout << "\n def "; for (int i=0;i<6;i++) cout << "("<<i<<")= " << STRAN[i];
        cout << "\n deltatdef "; for (int i=0;i<6;i++) cout << "("<<i<<")= " << DSTRAN[i];
        cout << "\n rotation "; for (int i=0;i<9;i++) cout << "("<<i<<")= " << DROT[i];
        cout << "\n gradiant_t "; for (int i=0;i<9;i++) cout << "("<<i<<")= " << DFGRD0[i];
        cout << "\n gradiant_tdt "; for (int i=0;i<9;i++) cout << "("<<i<<")= " << DFGRD1[i];
        cout << "\n temps_t= " << temps_t << " temps_tdt= " <<  temps_tdt 
             << " delta_t= " << delta_t;
        cout << "\n temper_t= " <<  temper_t << " delta_temper= " <<  delta_temper;
        cout << "\n nom_materiau= " << nom_materiau; // nom du matériau
        cout << "\n dim_tens= " << dim_tens << " nb_tau_ij= " << nb_tau_ij
             << " nb_ij_tenseur= " <<  nb_ij_tenseur; 
        cout << "\n nb_variables_etat= " << nb_variables_etat; 
        cout << "\n Coordonnee= " <<  coor_pt; // coordonnee du point
        cout << "\n mat_corota= " << mat_corota; // matrice de rotation du corotationnelle
        cout << "\n longueur_characteristique= " << longueur_characteristique;
        cout << "\n BaseB giB_t= " << giB_t; // base naturelle à t
        cout << "\n BaseB giB_tdt = " << giB_tdt; // base natuelle à tdt
        if (giB_t.NbVecteur()==2)
          cout <<"\n N_t= "<<N_t <<", N_tdt= "<<N_tdt;
        cout << "\n nb_elem= " << nb_elem ; // numéro de l'élément
        cout << " nb_pt_integ= " << nb_pt_integ; // numéro du point d'intégration
        cout << " nb_du_plis= " << nb_du_plis; // numéro du plis pour les multicouches 
        cout << " nb_pt_dans_le_plis= " << nb_pt_dans_le_plis; // numéro du po     cout << "dans le plis
        cout << " nb_step= " << nb_step; // numéro du step
        cout << " nb_increment= " << nb_increment; // numéro de l'incrément
        cout << "\n  ******** fin grandeurs envoyees ************** " << endl;
     	};
     close (tub); // fermeture du tampon
     };
  if (niveau > 5)
       cout << "\n HZpp-utilisation d'umat: fermeture du tube en ecriture" << endl;
 };
 
// Puis lecture des résultats de l'umat
void UmatAbaqus::LectureResultatUmat(bool utilisation_umat_interne,const int niveau)
{ //     Creation d'un processus de reception des données
  //     en fait l'objectif est de réceptionner que les variables calculées
  //     ouverture du tube nomme en lecture 
  //     lecture dans le t_car_x_n.tampon
  if (niveau > 5)
       cout << "\n hzpp-utilisation d'umat:  ouverture du tube en lecture, lecture sur le tube" 
       << endl; 
  if (!utilisation_umat_interne)
    { int tub = open(envoi.c_str(),O_RDONLY);  
      read (tub,t_car_x_n.tampon,480);
      close (tub); // fermeture du tampon
      // transformation de l'information 
      Copie_base_vers_evolue_PourUmat(niveau);
     };
  if (niveau > 5)
       cout << "\n HZpp-utilisation d'umat: fermeture du tube en lecture" << endl;
 };
 
// initialisation des données d'entrée tous à zéro
// sauf: nom_materiau: qui reste inchangé
//       giB_t et BaseB: qui sont initialisée à une base absolue identité
//       pnewdt: très grand, matrice de rotation identité
// sinon on a:
//       *eps_meca, *delta_eps_meca : nulles, temps_t=0, temps_tdt et delta_t=0;
//       temper_t et delta_temper: nulles,
//       pas de variables aux noeuds, pas de proprietes_materiau
//       coordonnées à 0, 
//       longueur caractéristique=0,
// et tous les indices=0: 
//       nb_elem nb_pt_integ nb_du_plis nb_pt_dans_le_plis nb_step nb_increment
void UmatAbaqus::Init_zero()
	{t_sigma->Inita(0.);v_statev.Change_taille(0);
	 energie_elastique=0.;dissipation_plastique=0.;dissipation_visqueuse=0.;
	 pnewdt = ConstMath::tresgrand;
	 eps_meca->Inita(0);delta_eps_meca->Inita(0.);
	 temps_t=0.;temps_tdt=0.;delta_t=0.;
	 temper_t=0.;delta_temper=0.;
	 ddlNoeA_ptinteg.Change_taille(0);
	 delta_ddlNoeA_ptinteg.Change_taille(0);
  
	 //dim_tens=3;nb_tau_ij=3;nb_ij_tenseur=6;
  dim_tens=Abs(eps_meca->Dimension());
  // les valeurs de nb_tau_ij et nb_ij_tenseur
  // dépendent du cas de travail, donc on ne les changes pas ici
  // éventuellement en lecture il seront modifié
  
	 nb_variables_etat=0.;
	 proprietes_materiau.Change_taille(0);nb_proprietes_materiau=0;
	 coor_pt.Zero();
	 mat_corota.Initialise(0.);mat_corota(1,1)=1.;mat_corota(2,2)=1.;mat_corota(3,3)=1.;
	 longueur_characteristique=0.;
	 giB_tdt = giB_t=BaseB(giB_t.Dimension(),giB_t.NbVecteur());
  if (giB_t.NbVecteur()==2)
    N_tdt.Zero(); N_t=N_tdt;
//  giB_tdt=BaseB(3);
  nb_elem=1;
  nb_pt_integ=1;
  nb_du_plis=1;
  nb_pt_dans_le_plis=1;
  nb_step=1;
  nb_increment=1;
	};
	
// initialisation des données d'entrée à un démarrage correcte
// *t_sigma=0, v_statev de taille nulle, les énergies nulles, pnewdt très grand
//       *eps_meca, *delta_eps_meca : nulles, temps_t=0, temps_tdt et delta_t=1;
//       temper_t et delta_temper: nulles,
//       pas de variables aux noeuds, pas de proprietes_materiau
//       coordonnées à 0, matrice de rotation identité
//       longueur caractéristique=1, giB_t et giB_tdt: base identité en absolue
// et tous les indices=1: 
//       nb_elem nb_pt_integ nb_du_plis nb_pt_dans_le_plis nb_step nb_increment
// par défaut la dimension des tenseurs est 3, mais pour les CP par exemple
// la dimension réelle utilisée est 2
// de même pour nb_reel_tau_ij qui indique le nombre de contrainte de cisaillement
// réellement utilisé: qui en CP = 1
void UmatAbaqus::Init_un(int dim_reel_tens ,int nb_reel_tau_ij)
	{t_sigma->Inita(0.);v_statev.Change_taille(0); 
	 energie_elastique=0.;dissipation_plastique=0.;dissipation_visqueuse=0.;
	 pnewdt = ConstMath::tresgrand;
	 eps_meca->Inita(0);delta_eps_meca->Inita(0.);
	 temps_t=0.;temps_tdt=1.;delta_t=1.;
	 temper_t=0.;delta_temper=0.;
	 ddlNoeA_ptinteg.Change_taille(0);
	 delta_ddlNoeA_ptinteg.Change_taille(0);
	 dim_tens=dim_reel_tens;
  nb_tau_ij=nb_reel_tau_ij;
  nb_ij_tenseur=dim_tens+nb_tau_ij; // correspond à NTENS,
  // en CP on n'utilise normalement que 4, mais de toute manière
  // on passe systématiquement des tenseurs 3D en paramètre
	 nb_variables_etat=0.;
	 proprietes_materiau.Change_taille(0);nb_proprietes_materiau=0;
	 coor_pt.Zero();
	 mat_corota.Initialise(0.);mat_corota(1,1)=1.;mat_corota(2,2)=1.;mat_corota(3,3)=1.;
	 longueur_characteristique=1.;
//	 giB_t=BaseB(3);giB_tdt=BaseB(3);
  giB_tdt = giB_t=BaseB(giB_t.Dimension(),giB_t.NbVecteur());
  if (giB_t.NbVecteur()==2)
    N_tdt.Zero(); N_t=N_tdt;

//     nb_elem=1;
//     nb_pt_integ=1;
  nb_du_plis=1;
  nb_pt_dans_le_plis=1;
//     nb_step=1;
//     nb_increment=1;
	};
		 
// transfert de l'Umat vers le point d'intégration à t
void UmatAbaqus::Transfert_Umat_ptInteg_t(PtIntegMecaInterne& ptIntegMeca)
{ *(ptIntegMeca.EpsBB()) = *eps_meca;
  *(ptIntegMeca.DeltaEpsBB()) = *delta_eps_meca;
  *(ptIntegMeca.SigHH_t()) = *t_sigma;
};	 
// transfert de l'Umat vers le point d'intégration à tdt
void UmatAbaqus::Transfert_Umat_ptInteg_tdt(PtIntegMecaInterne& ptIntegMeca)
{ *(ptIntegMeca.SigHH()) = *t_sigma;
};	 
	

// affichage en fonction d'un niveau
// niveau > 0: données
void UmatAbaqus::Affiche(ofstream& sort,const int niveau)
 { sort << "\n *********** info umat ************ ";
   sort << "\n sigma_t "; for (int i=0;i<6;i++) sort << "("<<i<<")= " << STRESS[i];
   sort << "\n def "; for (int i=0;i<6;i++) sort << "("<<i<<")= " << STRAN[i];
   sort << "\n deltatdef "; for (int i=0;i<6;i++) sort << "("<<i<<")= " << DSTRAN[i];
   sort << "\n rotation "; for (int i=0;i<9;i++) sort << "("<<i<<")= " << DROT[i];
   sort << "\n gradiant_t "; for (int i=0;i<9;i++) sort << "("<<i<<")= " << DFGRD0[i];
   sort << "\n gradiant_tdt "; for (int i=0;i<9;i++) sort << "("<<i<<")= " << DFGRD1[i];
   sort << "\n temps_t= " << temps_t << " temps_tdt= " <<  temps_tdt
        << " delta_t= " << delta_t;
   sort << "\n temper_t= " <<  temper_t << " delta_temper= " <<  delta_temper;
   sort << "\n nom_materiau= " << nom_materiau; // nom du matériau
   sort << "\n dim_tens= " << dim_tens << " nb_tau_ij= " << nb_tau_ij
        << " nb_ij_tenseur= " <<  nb_ij_tenseur;
   sort << "\n nb_variables_etat= " << nb_variables_etat;
   sort << "\n Coordonnee= " <<  coor_pt; // coordonnee du point
   sort << "\n mat_corota= " << mat_corota; // matrice de rotation du corotationnelle
   sort << "\n longueur_characteristique= " << longueur_characteristique;
   sort << "\n BaseB giB_t= " << giB_t; // base naturelle à t
   sort << "\n BaseB giB_tdt = " << giB_tdt; // base natuelle à tdt
   if (giB_t.NbVecteur()==2)
     sort <<"\n N_t= "<<N_t <<", N_tdt= "<<N_tdt;
   sort << "\n nb_elem= " << nb_elem ; // numéro de l'élément
   sort << " nb_pt_integ= " << nb_pt_integ; // numéro du point d'intégration
   sort << " nb_du_plis= " << nb_du_plis; // numéro du plis pour les multicouches
   sort << " nb_pt_dans_le_plis= " << nb_pt_dans_le_plis; // numéro du po     sort << "dans le plis
   sort << " nb_step= " << nb_step; // numéro du step
   sort << " nb_increment= " << nb_increment; // numéro de l'incrément
   sort << "\n  ******** fin grandeurs  ************** " << endl;
 };

	//----- lecture écriture de restart -----
	// cas donne le niveau de la récupération
       // = 1 : on récupère tout
       // = 2 : on récupère uniquement les données variables (supposées comme telles)
void UmatAbaqus::Lecture_base_info(ifstream& ent,const int cas)
	{ string toto,nom; 
	  if (cas == 1) 
	   { ent >> toto;
	     if (toto != "StructureUmat_Abaqus")
        { cout << "\n erreur en lecture des parametres d'une structure Umat abaqus !! "
               << "\n UmatAbaqus::Lecture_base_info(... ";
          Sortie(1);
        };
      ent >> toto >> envoi >> toto >> reception;
    };
 };

       // cas donne le niveau de sauvegarde
       // = 1 : on sauvegarde tout
       // = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void UmatAbaqus::Ecriture_base_info(ofstream& sort,const int cas)
 { if (cas == 1)
       { sort << "\n StructureUmat_Abaqus: " 
              << " nom_tube_envoi " << envoi <<" nom_tube_reception= " << reception << " ";  
       }; 
	};

// ---------------------------------------------------------------------------------------------
//                        méthodes protégées 
// ---------------------------------------------------------------------------------------------

    // 1)-- utilisation en tant qu'umat
// copie des grandeurs évoluées dans les grandeurs de base
// uniquement celles qui varient
void UmatAbaqus::Copie_evolue_vers_base(const int niveau)
{  if ((dim_tens == 3) && (nb_tau_ij == 3))
	   { // cas classique 3D
      Tenseur3HH&  sigma =  *((Tenseur3HH*) t_sigma);
	     Tenseur3HHHH& dd_sigma_deps =  *((Tenseur3HHHH*) d_sigma_deps);
      Tenseur3HH&  d_sigma_d_tempe =  *((Tenseur3HH*) d_sigma_d_tempera);
      Tenseur3HH&  d_ener_therm_vol_dep =  *((Tenseur3HH*) d_ener_therm_vol_deps);
      for (int ij=1;ij<=6;ij++)
       {int i=cdex.idx_i(ij);int j=cdex.idx_j(ij);
        STRESS[ij-1] = sigma(i,j); 
        DDSDDT[ij-1] = d_sigma_d_tempe(i,j);
        DRPLDE[ij-1] = d_ener_therm_vol_dep(i,j);
        for (int kl=1;kl<=6;kl++)
         {int k=cdex.idx_i(kl);int l=cdex.idx_j(kl);
          int r=(ij-1)*6+kl-1;		   
          //		   DDSDDE[r] = dd_sigma_deps(i,j,k,l);
          // en considérant que la matrice tangente est relative à gamma_ij = eps_ij + eps_ji
          // et non la déformation !!, de plus en considérant que dd_sigma_deps (i,j,k,l)
          // devrait être symétrique par rapport aux deux derniers indices
          // enfin dans l'utilisation on fait directement sig(i,j) = dd(i,j,k,l) * gamma(k,l)
          // et on ne considère pas le cas gamma(l,k), donc il ya seulement 6 gamma
          if (k == l)
           {DDSDDE[r] = dd_sigma_deps(i,j,k,l);}
          else 
     //		    {DDSDDE[r] = 2. * dd_sigma_deps(i,j,k,l);};
           {DDSDDE[r] = 0.5 * (dd_sigma_deps(i,j,k,l)+dd_sigma_deps(i,j,l,k));};
         };
       };
      if (niveau > 5)
  	     {  cout << "\n ************* contraintes calculees: *********\n" << sigma;
            cout << "\n ************* matrice tangente ***************\n" << dd_sigma_deps;
  	     };
	   }
   else if ((dim_tens == 2) && (nb_tau_ij == 1))
    { // cas des contraintes planes
      Tenseur2HH&  sigma =  *((Tenseur2HH*) t_sigma);
      Tenseur2HHHH& dd_sigma_deps =  *((Tenseur2HHHH*) d_sigma_deps);
      Tenseur2HH&  d_sigma_d_tempe =  *((Tenseur2HH*) d_sigma_d_tempera);
      Tenseur2HH&  d_ener_therm_vol_dep =  *((Tenseur2HH*) d_ener_therm_vol_deps);
     
      for (int ij=1;ij<=3;ij++)
       {int i=cdexCP.idx_i(ij);int j=cdexCP.idx_j(ij);
        STRESS[ij-1] = sigma(i,j);
        DDSDDT[ij-1] = d_sigma_d_tempe(i,j);
        DRPLDE[ij-1] = d_ener_therm_vol_dep(i,j);
        for (int kl=1;kl<=3;kl++)
         {int k=cdexCP.idx_i(kl);int l=cdexCP.idx_j(kl);
          int r=(ij-1)*3+kl-1;
          //     DDSDDE[r] = dd_sigma_deps(i,j,k,l);
          // en considérant que la matrice tangente est relative à gamma_ij = eps_ij + eps_ji
          // et non la déformation !!, de plus en considérant que dd_sigma_deps (i,j,k,l)
          // devrait être symétrique par rapport aux deux derniers indices
          // enfin dans l'utilisation on fait directement sig(i,j) = dd(i,j,k,l) * gamma(k,l)
          // et on ne considère pas le cas gamma(l,k), donc il ya seulement 3 gamma
          if (k == l)
           {DDSDDE[r] = dd_sigma_deps(i,j,k,l);}
          else
           {DDSDDE[r] = 0.5 * (dd_sigma_deps(i,j,k,l)+dd_sigma_deps(i,j,l,k));};
         };
       };
      // on s'occupe de mettre à jour la normale
      int i=3;
      for (int j=1;j<=3;j++)
       {int r=(i-1)*3+j-1;
        DFGRD1[r] = N_tdt(j);
       };
     
      if (niveau > 5)
        {  cout << "\n ************* contraintes calculees: *********\n" << sigma;
           cout << "\n ************* matrice tangente ***************\n" << dd_sigma_deps;
           cout <<"\n N_t= "<<N_t <<", N_tdt= "<<N_tdt;
        };
    }
   else if ((dim_tens == 3) && (nb_tau_ij == 1))
	   { // cas d'éléments axisymétrique 3D
      Tenseur3HH&  sigma =  *((Tenseur3HH*) t_sigma);
	     Tenseur3HHHH& dd_sigma_deps =  *((Tenseur3HHHH*) d_sigma_deps);
      Tenseur3HH&  d_sigma_d_tempe =  *((Tenseur3HH*) d_sigma_d_tempera);
      Tenseur3HH&  d_ener_therm_vol_dep =  *((Tenseur3HH*) d_ener_therm_vol_deps);
      for (int ij=1;ij<=4;ij++)
       {int i=cdex.idx_i(ij);int j=cdex.idx_j(ij);
        STRESS[ij-1] = sigma(i,j); 
        DDSDDT[ij-1] = d_sigma_d_tempe(i,j);
        DRPLDE[ij-1] = d_ener_therm_vol_dep(i,j);
        for (int kl=1;kl<=4;kl++)
         {int k=cdex.idx_i(kl);int l=cdex.idx_j(kl);
          int r=(ij-1)*4+kl-1;
          //		   DDSDDE[r] = dd_sigma_deps(i,j,k,l);
          // en considérant que la matrice tangente est relative à gamma_ij = eps_ij + eps_ji
          // et non la déformation !!, de plus en considérant que dd_sigma_deps (i,j,k,l)
          // devrait être symétrique par rapport aux deux derniers indices
          // enfin dans l'utilisation on fait directement sig(i,j) = dd(i,j,k,l) * gamma(k,l)
          // et on ne considère pas le cas gamma(l,k), donc il ya seulement 6 gamma
          if (k == l)
           {DDSDDE[r] = dd_sigma_deps(i,j,k,l);}
          else 
     //		    {DDSDDE[r] = 2. * dd_sigma_deps(i,j,k,l);};
           {DDSDDE[r] = 0.5 * (dd_sigma_deps(i,j,k,l)+dd_sigma_deps(i,j,l,k));};
         };
       };
      if (niveau > 5)
  	     {  cout << "\n ************* contraintes calculees: *********\n" << sigma;
           cout << "\n ************* matrice tangente ***************\n" << dd_sigma_deps;
  	     };
	   }
   else
    { cout << "\n erreur , pour l'instant seules la dimension 3 et le cas axi sont en place"
           << "\n il est vivement conseille de ce plaindre pour remedier a ce probleme "
           << "\n void UmatAbaqus::Copie_evolue_vers_base()" << endl ;
      Sortie(1);     
    };
   *SSE = energie_elastique; *SPD = dissipation_plastique; *SCD = dissipation_visqueuse; 
   *PNEWDT = pnewdt;
   *RPL =  ener_therm_vol_par_utemps; 
   *DRPLDT = d_ener_therm_dtemper; 
//cout << "\n sigma_t "; for (int i=0;i<6;i++) cout << " ("<<i<<")=" << STRESS[i];
//cout << endl;
};

// copie des grandeurs de base vers les grandeurs évoluées
// uniquement les grandeurs d'entrée (supposées à lire)
void UmatAbaqus::Copie_base_vers_evolue(const int niveau)
{  if ((*NDI == 3) && (*NSHR == 3))
    { // cas classique 3D
      Tenseur3HH&  sigma =  *((Tenseur3HH*) t_sigma);
      Tenseur3BB&  eps_m =  *((Tenseur3BB*) eps_meca);
      Tenseur3BB&  delta_eps_m =  *((Tenseur3BB*) delta_eps_meca);
      for (int ij=1;ij<=6;ij++)
      {int i=cdex.idx_i(ij);int j=cdex.idx_j(ij);
       sigma.Coor(i,j) = STRESS[ij-1];
       // transformation des gamma en epsilon pour les termes ij, i diff de j
       if (i==j){eps_m.Coor(i,j) = STRAN[ij-1] + DSTRAN[ij-1];
                 delta_eps_m.Coor(i,j) = DSTRAN[ij-1];}
       else     {eps_m.Coor(i,j) = 0.5 * (STRAN[ij-1] + DSTRAN[ij-1]);
                 delta_eps_m.Coor(i,j) = 0.5 * DSTRAN[ij-1];}
      }
    }
   else if ((*NDI == 3) && (*NSHR == 1))
    { // cas d'éléments axisymétrique 3D
      Tenseur3HH&  sigma =  *((Tenseur3HH*) t_sigma);
      Tenseur3BB&  eps_m =  *((Tenseur3BB*) eps_meca);
      Tenseur3BB&  delta_eps_m =  *((Tenseur3BB*) delta_eps_meca);
      for (int ij=1;ij<=4;ij++)
        {int i=cdex.idx_i(ij);int j=cdex.idx_j(ij);
         sigma.Coor(i,j) = STRESS[ij-1];
         // transformation des gamma en epsilon pour les termes ij, i diff de j
         if (i==j){eps_m.Coor(i,j) = STRAN[ij-1] + DSTRAN[ij-1];
                   delta_eps_m.Coor(i,j) = DSTRAN[ij-1];}
         else     {eps_m.Coor(i,j) = 0.5 * (STRAN[ij-1] + DSTRAN[ij-1]);
                   delta_eps_m.Coor(i,j) = 0.5 * DSTRAN[ij-1];}
        };
      // maintenant je m'occupe des autres valeurs qui sont nulles
      sigma.Coor(1,3)=sigma.Coor(2,3)=0.;
      eps_m.Coor(1,3)=eps_m.Coor(2,3)=0.;
      delta_eps_m.Coor(1,3)=delta_eps_m.Coor(2,3)=0.;
    }
   else if ((*NDI == 2) && (*NSHR == 1))
    { // cas d'éléments en contraintes plane
      Tenseur2HH&  sigma =  *((Tenseur2HH*) t_sigma);
      Tenseur2BB&  eps_m =  *((Tenseur2BB*) eps_meca);
      Tenseur2BB&  delta_eps_m =  *((Tenseur2BB*) delta_eps_meca);
      for (int ij=1;ij<=3;ij++)
        {int i=cdexCP.idx_i(ij);int j=cdexCP.idx_j(ij);
         sigma.Coor(i,j) = STRESS[ij-1];
         // transformation des gamma en epsilon pour les termes ij, i diff de j
         if (i==j){eps_m.Coor(i,j) = STRAN[ij-1] + DSTRAN[ij-1];
                   delta_eps_m.Coor(i,j) = DSTRAN[ij-1];}
         else     {eps_m.Coor(i,j) = 0.5 * (STRAN[ij-1] + DSTRAN[ij-1]);
                   delta_eps_m.Coor(i,j) = 0.5 * DSTRAN[ij-1];}
        };
      // maintenant je m'occupe des autres valeurs qui sont nulles
//      sigma(1,3)=sigma(2,3)=sigma(3,3)=0.;
//      eps_m(1,3)=eps_m(2,3)=0.;
//      delta_eps_m(1,3)=delta_eps_m(2,3)=0.;
    }
   else
    { cout << "\n erreur , pour l'instant seules la dimension 3, le cas contrainte plane et le cas axi sont en place"
           << "\n il est vivement conseille de ce plaindre pour remedier a ce probleme "
           << "\n void UmatAbaqus::Copie_base_vers_evolue()" << endl ;
      Sortie(1);     
    };
   energie_elastique = dissipation_plastique = dissipation_visqueuse = 0.;  // init
   pnewdt=  *PNEWDT;
   temps_t= TIME[1];  delta_t = *DTIME; temps_tdt= temps_t+delta_t;
   temper_t =  *TemP;  delta_temper = *DTEMP;
   dim_tens =  *NDI; nb_tau_ij =  *NSHR; nb_ij_tenseur =  *NTENS; nb_variables_etat =  *NSTATV;
   for (int i=1;i<=3;i++) 
      {coor_pt(i)= COORDS[i-1];};
   // au niveau des bases il faut intégrer le cas des contraintes planes
   if ((*NDI == 2) && (*NSHR == 1)) // cas contraintes planes
    {for (int i=1;i<=2;i++) // uniquement 2 vecteurs pour les gi
      {for (int j=1;j<=3;j++)
        {int r=(i-1)*3+j-1;
         mat_corota(i,j)= DROT[r]; giB_t.CoordoB(i)(j)=DFGRD0[r]; giB_tdt.CoordoB(i)(j)=DFGRD1[r];
        };
      };
     // on stocke les normales
     int i=3;
     for (int j=1;j<=3;j++)
      {int r=(i-1)*3+j-1;
       N_t(j) = DFGRD0[r];
       N_tdt(j) = DFGRD1[r];
      };
    }
   else // autres cas on est en 3D avec 3 vecteurs
    {for (int i=1;i<=3;i++)
      {for (int j=1;j<=3;j++)
        {int r=(i-1)*3+j-1;
         mat_corota(i,j)= DROT[r]; giB_t.CoordoB(i)(j)=DFGRD0[r]; giB_tdt.CoordoB(i)(j)=DFGRD1[r];
        };
      };
    };
 
   longueur_characteristique = *CELENT; 
   nb_elem = *NOEL; nb_pt_integ = *NPT; nb_du_plis = *LAYER; nb_pt_dans_le_plis = *KSPT;
   nb_step = *KSTEP; nb_increment = *KINC;
   if (nb_increment <= 1)
      nom_materiau=CMNAME; // le nom du matériau
};

    // 2) -- utilisation en tant qu'utilisation d 'une umat extérieure
// copie des grandeurs de base vers les grandeurs évoluées
// uniquement celles qui sont des résultats
void UmatAbaqus::Copie_base_vers_evolue_PourUmat(const int niveau)
{  if ((*NDI == 3) && (*NSHR == 3))
	   { // cas classique 3D
      Tenseur3HH&  sigma =  *((Tenseur3HH*) t_sigma);
	     Tenseur3HHHH& dd_sigma_deps =  *((Tenseur3HHHH*) d_sigma_deps);
      Tenseur3HH&  d_sigma_d_tempe =  *((Tenseur3HH*) d_sigma_d_tempera);
      Tenseur3HH&  d_ener_therm_vol_dep =  *((Tenseur3HH*) d_ener_therm_vol_deps);
      for (int ij=1;ij<=6;ij++)
       {int i=cdex.idx_i(ij);int j=cdex.idx_j(ij);
        sigma.Coor(i,j) = STRESS[ij-1];
        d_sigma_d_tempe.Coor(i,j) = DDSDDT[ij-1];
        d_ener_therm_vol_dep.Coor(i,j) = DRPLDE[ij-1];
        for (int kl=1;kl<=6;kl++)
         {int k=cdex.idx_i(kl);int l=cdex.idx_j(kl);
          int r=(ij-1)*6+kl-1;
          // en considérant que la matrice tangente est relative à gamma_ij = eps_ij + eps_ji
          // et non la déformation !!, de plus en considérant que dd_sigma_deps (i,j,k,l)
          // devrait être symétrique par rapport aux deux derniers indices
          // enfin dans l'utilisation on fait directement sig(i,j) = dd(i,j,k,l) * gamma(k,l)
          // et on ne considère pas le cas gamma(l,k), donc il ya seulement 6 gamma
           dd_sigma_deps.Change(i,j,k,l,DDSDDE[r]);
//          if (k == l)
//           {dd_sigma_deps.Change(i,j,k,l,DDSDDE[r]);}
//          else
//           {dd_sigma_deps.Change(i,j,k,l,0.5*DDSDDE[r]);};
         };
       };
      if (niveau > 5)
  	     {  cout << "\n ************* contraintes calculees: *********\n" << sigma;
           cout << "\n ************* matrice tangente ***************\n" << dd_sigma_deps;
  	     };
	   }
   else if ((*NDI == 2) && (*NSHR == 1))
    { // cas des contraintes planes
      Tenseur2HH&  sigma =  *((Tenseur2HH*) t_sigma);
      Tenseur2HHHH& dd_sigma_deps =  *((Tenseur2HHHH*) d_sigma_deps);
      Tenseur2HH&  d_sigma_d_tempe =  *((Tenseur2HH*) d_sigma_d_tempera);
      Tenseur2HH&  d_ener_therm_vol_dep =  *((Tenseur2HH*) d_ener_therm_vol_deps);
      for (int ij=1;ij<=3;ij++)
       {int i=cdexCP.idx_i(ij);int j=cdexCP.idx_j(ij);
        sigma.Coor(i,j) = STRESS[ij-1];
        d_sigma_d_tempe.Coor(i,j) = DDSDDT[ij-1];
        d_ener_therm_vol_dep.Coor(i,j) = DRPLDE[ij-1];
        for (int kl=1;kl<=3;kl++)
         {int k=cdexCP.idx_i(kl);int l=cdexCP.idx_j(kl);
          int r=(ij-1)*3+kl-1;
          dd_sigma_deps.Change(i,j,k,l,DDSDDE[r]);
         };
       };
      // on met également à jour la normale finale
      // ce qui permet de récupérer ensuite la déformation d'épaisseur
      int i=3;
      for (int j=1;j<=3;j++)
       {int r=(i-1)*3+j-1;
        N_tdt(j) = DFGRD1[r];
       };
      if (niveau > 5)
        {  cout << "\n ************* contraintes calculees: *********\n" << sigma;
           cout << "\n ************* matrice tangente ***************\n" << dd_sigma_deps;
           cout <<"\n N_t= "<<N_t <<", N_tdt= "<<N_tdt;
        };
    }
   else if ((*NDI == 3) && (*NSHR == 1))
	   { // cas d'éléments axisymétrique 3D
      Tenseur3HH&  sigma =  *((Tenseur3HH*) t_sigma);
	     Tenseur3HHHH& dd_sigma_deps =  *((Tenseur3HHHH*) d_sigma_deps);
      Tenseur3HH&  d_sigma_d_tempe =  *((Tenseur3HH*) d_sigma_d_tempera);
      Tenseur3HH&  d_ener_therm_vol_dep =  *((Tenseur3HH*) d_ener_therm_vol_deps);
      for (int ij=1;ij<=4;ij++)
       {int i=cdex.idx_i(ij);int j=cdex.idx_j(ij);
        sigma.Coor(i,j) = STRESS[ij-1];
        d_sigma_d_tempe.Coor(i,j) = DDSDDT[ij-1];
        d_ener_therm_vol_dep.Coor(i,j) = DRPLDE[ij-1];
        for (int kl=1;kl<=4;kl++)
         {int k=cdex.idx_i(kl);int l=cdex.idx_j(kl);
          int r=(ij-1)*4+kl-1;
          dd_sigma_deps.Change(i,j,k,l,DDSDDE[r]);
         };
       };
      // maintenant je m'occupe des autres valeurs qui sont nulles
      sigma.Coor(1,3)=sigma.Coor(2,3)=0.;
      d_sigma_d_tempe.Coor(1,3)=d_sigma_d_tempe.Coor(2,3)=0.;
      d_ener_therm_vol_dep.Coor(1,3)=d_ener_therm_vol_dep.Coor(2,3)=0.;
      if (niveau > 5)
  	     {  cout << "\n ************* contraintes calculees: *********\n" << sigma;
           cout << "\n ************* matrice tangente ***************\n" << dd_sigma_deps;
  	     };
	   }
   else
    { cout << "\n erreur , pour l'instant seules la dimension 3 et le cas axi sont en place"
           << "\n il est vivement conseille de ce plaindre pour remedier a ce probleme "
           << "\n void UmatAbaqus::Copie_evolue_vers_base()" << endl ;
      Sortie(1);     
    };
   energie_elastique = *SSE; dissipation_plastique = *SPD; dissipation_visqueuse = *SCD; 
   pnewdt = *PNEWDT ;
   ener_therm_vol_par_utemps = *RPL ; 
   d_ener_therm_dtemper = *DRPLDT; 
   if (niveau > 5)
    { cout << "\n sigma_t ";
      for (int i=0;i<(*NDI+*NSHR);i++) cout << " ("<<i<<")=" << STRESS[i];
      cout << endl;
    };
};

// copie des grandeurs évoluées dans les grandeurs de base
// uniquement celles qui sont des données
void UmatAbaqus::Copie_evolue_vers_base_PourUmat(const int niveau)
{  if ((dim_tens == 3) && (nb_tau_ij == 3))
	   { // cas classique 3D
      Tenseur3HH&  sigma =  *((Tenseur3HH*) t_sigma);
	     Tenseur3BB&  eps_m =  *((Tenseur3BB*) eps_meca);
	     Tenseur3BB&  delta_eps_m =  *((Tenseur3BB*) delta_eps_meca);
      for (int ij=1;ij<=6;ij++)
       {int i=cdex.idx_i(ij);int j=cdex.idx_j(ij);
        STRESS[ij-1] = sigma(i,j);
        // passage des déformations en gamma pour les termes ij, i diff  de j
        if (i == j)
            {STRAN[ij-1] = eps_m(i,j) - delta_eps_m(i,j);
             DSTRAN[ij-1] = delta_eps_m(i,j);
            }
        else    
            {STRAN[ij-1] = 2. * (eps_m(i,j) - delta_eps_m(i,j));
             DSTRAN[ij-1] = 2. * delta_eps_m(i,j);
           };
       };
    }
   else if ((dim_tens == 2) && (nb_tau_ij == 1))
    { // cas des contraintes planes
      Tenseur2HH&  sigma =  *((Tenseur2HH*) t_sigma);
      Tenseur2BB&  eps_m =  *((Tenseur2BB*) eps_meca);
      Tenseur2BB&  delta_eps_m =  *((Tenseur2BB*) delta_eps_meca);
      for (int ij=1;ij<=3;ij++)
       {int i=cdexCP.idx_i(ij);int j=cdexCP.idx_j(ij);
        STRESS[ij-1] = sigma(i,j);
        // passage des déformations en gamma pour les termes ij, i diff  de j
        if (i == j)
            {STRAN[ij-1] = eps_m(i,j) - delta_eps_m(i,j);
             DSTRAN[ij-1] = delta_eps_m(i,j);
            }
        else
            {STRAN[ij-1] = 2. * (eps_m(i,j) - delta_eps_m(i,j));
             DSTRAN[ij-1] = 2. * delta_eps_m(i,j);
           };
       };
    }
   else if ((dim_tens == 3) && (nb_tau_ij == 1))
	   { // cas d'éléments axisymétrique 3D
      Tenseur3HH&  sigma =  *((Tenseur3HH*) t_sigma);
	     Tenseur3BB&  eps_m =  *((Tenseur3BB*) eps_meca);
	     Tenseur3BB&  delta_eps_m =  *((Tenseur3BB*) delta_eps_meca);
      for (int ij=1;ij<=4;ij++)
       {int i=cdex.idx_i(ij);int j=cdex.idx_j(ij);
        STRESS[ij-1] = sigma(i,j);
        // passage des déformations en gamma pour les termes ij, i diff  de j
        if (i == j)
            {STRAN[ij-1] = eps_m(i,j) - delta_eps_m(i,j);
             DSTRAN[ij-1] = delta_eps_m(i,j);
            }
        else    
            {STRAN[ij-1] = 2. * (eps_m(i,j) - delta_eps_m(i,j));
             DSTRAN[ij-1] = 2. * delta_eps_m(i,j);
           };
       };
    }
   else
    { cout << "\n erreur , pour l'instant seules la dimension 3 et le cas axi sont en place"
           << "\n il est vivement conseille de ce plaindre pour remedier a ce probleme "
           << "\n void UmatAbaqus::Copie_base_vers_evolue()" << endl ;
      Sortie(1);     
    };
   *PNEWDT = pnewdt;
   TIME[0] = temps_t; // par défaut on met un step = la durée totale
   TIME[1] = temps_t; *DTIME = delta_t;
   *TemP = temper_t ;  *DTEMP = delta_temper ;
   *NDI = dim_tens ; *NSHR = nb_tau_ij ; *NTENS = nb_ij_tenseur ; *NSTATV = nb_variables_etat;
   for (int i=1;i<=3;i++)
      {COORDS[i-1] = coor_pt(i);};
   // au niveau des bases il faut intégrer le cas des contraintes planes
   if ((*NDI == 2) && (*NSHR == 1)) // cas contraintes planes
    { for (int i=1;i<=2;i++) // uniquement 2 vecteurs pour les gi
       {for (int j=1;j<=3;j++)
           {int r=(i-1)*3+j-1;
            DROT[r] = mat_corota(i,j); DFGRD0[r] = giB_t (i)(j); DFGRD1[r] =giB_tdt(i)(j);
           };
       };
      // pour la troisième direction, mat_corota normalement est bien initialisé
      // et on init les autres à un vecteur unitaire normale
      int i=3; // le cas de la troisième direction
       {for (int j=1;j<=3;j++)
           {int r=(i-1)*3+j-1;
            DROT[r] = mat_corota(i,j);
            DFGRD0[r] = N_t(j);
            DFGRD1[r] = N_tdt(j);
           };
       };
    }
   else // autres cas on est en 3D avec 3 vecteurs
    {for (int i=1;i<=3;i++)
         {COORDS[i-1] = coor_pt(i);
          for (int j=1;j<=3;j++)
           {int r=(i-1)*3+j-1;
            DROT[r] = mat_corota(i,j); DFGRD0[r] = giB_t (i)(j); DFGRD1[r] =giB_tdt(i)(j);
           };
         };
    };
 
   *CELENT = longueur_characteristique ; 
   *NOEL = nb_elem ; *NPT = nb_pt_integ ; *LAYER = nb_du_plis ;  *KSPT = nb_pt_dans_le_plis ;
   *KSTEP = nb_step; *KINC = nb_increment;
   if (nb_increment <= 1)
     { int longueur = nom_materiau.length();
//	   longueur = 5;
       if (longueur > 24)
       	{cout << "\n erreur: le nom du materiau a une longueur superieur a 24, pas possible actuellement"
       	      << "\n UmatAbaqus::Copie_evolue_vers_base_PourUmat()";
       	 Sortie(1);     
       	};
        for (int i=0;i<longueur;i++)
          	CMNAME[i]=(nom_materiau.c_str())[i];
        CMNAME[longueur]='\0'; 	
     };
};