// FICHIER : Loi_maxwell1D.cp
// CLASSE : Loi_maxwell1D


// 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-2021 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 "Debug.h"

# include <iostream>
using namespace std;  //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "TypeConsTens.h"
#include "ParaGlob.h"
#include "ConstMath.h"
#include "CharUtil.h"

#include "Loi_maxwell1D.h"


Loi_maxwell1D::Loi_maxwell1D () :  // Constructeur par defaut
  Loi_comp_abstraite(MAXWELL1D,CAT_THERMO_MECANIQUE,1),E(-ConstMath::trespetit),mu(-ConstMath::trespetit)
  ,xn(1.),simple(true)
  ,E_temperature(NULL),mu_temperature(NULL),xn_temperature(NULL)
  ,type_derive(-1)
   { };
	

// Constructeur de copie
Loi_maxwell1D::Loi_maxwell1D (const Loi_maxwell1D& loi) :
	Loi_comp_abstraite(loi),E(loi.E),mu(loi.mu),type_derive(loi.type_derive)
	,xn(loi.xn),simple(loi.xn)
    ,E_temperature(loi.E_temperature),mu_temperature(loi.mu_temperature)
    ,xn_temperature(loi.xn_temperature)
	 { // on regarde s'il s'agit d'une courbe locale ou d'une courbe globale
    if (E_temperature != NULL)
	  	if (E_temperature->NomCourbe() == "_") 
	  		  E_temperature = Courbe1D::New_Courbe1D(*(loi.E_temperature));
    if (mu_temperature != NULL)	
      	if (mu_temperature->NomCourbe() == "_") 
      	  mu_temperature = Courbe1D::New_Courbe1D(*(loi.mu_temperature));
    if (xn_temperature != NULL)	
       if (xn_temperature->NomCourbe() == "_") 
            xn_temperature = Courbe1D::New_Courbe1D(*(loi.xn_temperature));;
	  };

Loi_maxwell1D::~Loi_maxwell1D ()
// Destructeur 
{ if (E_temperature != NULL)
      if (E_temperature->NomCourbe() == "_") delete E_temperature;
  if (mu_temperature != NULL)
      if (mu_temperature->NomCourbe() == "_") delete mu_temperature;
  if (xn_temperature != NULL)
      if (xn_temperature->NomCourbe() == "_") delete xn_temperature;
  };

// Lecture des donnees de la classe sur fichier
void Loi_maxwell1D::LectureDonneesParticulieres (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D
                                             ,LesFonctions_nD& lesFonctionsnD)
  { // module d'young
    string nom;*(entreePrinc->entree) >> nom; 
    if (nom != "E=") 
         { cout << "\n erreur en lecture du module d'young, on aurait du lire le mot E=";
           entreePrinc->MessageBuffer("**erreur1 Loi_maxwell1D::LectureDonneesParticulieres (...**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
         };
    // on regarde si le module d'young est thermo dépendant
    if(strstr(entreePrinc->tablcar,"E_thermo_dependant_")!=0)
     { thermo_dependant=true;
       *(entreePrinc->entree) >> nom;
       if (nom != "E_thermo_dependant_") 
         { cout << "\n erreur en lecture de la thermodependance de E, on aurait du lire le mot cle E_thermo_dependant_"
                << " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
           entreePrinc->MessageBuffer("**erreur2 Loi_maxwell1D::LectureDonneesParticulieres (...**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
         };
       // lecture de la loi d'évolution du module d'young en fonction de la température
       *(entreePrinc->entree) >>  nom;    
       // on regarde si la courbe existe, si oui on récupère la référence 
       if (lesCourbes1D.Existe(nom)) 
         { E_temperature = lesCourbes1D.Trouve(nom);
         }
       else
        { // sinon il faut la lire maintenant
          string non_courbe("_");   
          E_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
          // lecture de la courbe
          E_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
        };
       entreePrinc->NouvelleDonnee(); // prepa du flot de lecture
     }
    else
     { // lecture du module d'young 
       *(entreePrinc->entree) >> E ;
     }; 

    // lecture de la viscosité 
    *(entreePrinc->entree) >> nom; 
    if (nom != "mu=") 
         { cout << "\n erreur en lecture de la viscosite, on aurait du lire le mot mu=";
           entreePrinc->MessageBuffer("**erreur3 Loi_maxwell1D::LectureDonneesParticulieres (...**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
         };
    // on regarde si la viscosité est thermo dépendante
    if(strstr(entreePrinc->tablcar,"mu_thermo_dependant_")!=0)
     { thermo_dependant=true;
       *(entreePrinc->entree) >> nom;
       if (nom != "mu_thermo_dependant_") 
         { cout << "\n erreur en lecture de la thermodependance de mu, on aurait du lire le mot cle mu_thermo_dependant_"
                << " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
           entreePrinc->MessageBuffer("**erreur4 Loi_maxwell1D::LectureDonneesParticulieres (...**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
         };
       // lecture de la loi d'évolution en fonction de la température
       *(entreePrinc->entree) >>  nom;    
       // on regarde si la courbe existe, si oui on récupère la référence 
       if (lesCourbes1D.Existe(nom)) 
         { mu_temperature = lesCourbes1D.Trouve(nom);
         }
       else
        { // sinon il faut la lire maintenant
          string non_courbe("_");   
          mu_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
          // lecture de la courbe
          mu_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
        };
       // prepa du flot de lecture
       if(strstr(entreePrinc->tablcar,"fin_coeff_MAXWELL1D")==0)
          entreePrinc->NouvelleDonnee();
     }
    else
     { // lecture de mu 
       *(entreePrinc->entree) >> mu ;
     }; 
    // on regarde ensuite si le type de dérivée est indiqué
    string toto;
    if (strstr(entreePrinc->tablcar,"type_derivee")!=NULL)
     { // lecture du type
       *(entreePrinc->entree) >> toto >> type_derive;
       if ((type_derive!=0)&&(type_derive!=-1)&&(type_derive!=1))
         { cout << "\n le type de derivee indique pour la loi de maxwell1D:  "<< type_derive
                << " n'est pas acceptable (uniquement -1 ou 0 ou 1), on utilise le type par defaut (-1)"
                << " qui correspon à la derivee de Liedeux fois covariantes";
           type_derive = -1;
         };
     };
    // on regarde s'il y a un coefficient non linéaire
    simple = true;   // par défaut
    if (strstr(entreePrinc->tablcar,"xn=")!=NULL)
      { *(entreePrinc->entree) >> toto ; //>> xn ;
        simple = false;
        #ifdef MISE_AU_POINT
          if (toto != "xn=")
             { cout << "\n erreur en lecture de la loi de maxwell1D, on attendait xn= et un nombre ";
               entreePrinc->MessageBuffer("**erreur5 Loi_maxwell1D::LectureDonneesParticulieres (...**");
               throw (UtilLecture::ErrNouvelleDonnee(-1));
               Sortie(1);
             };
        #endif
        // on regarde si le coefficient non linéaire  est thermo dépendant
        if(strstr(entreePrinc->tablcar,"xn_thermo_dependant_")!=0)
         { thermo_dependant=true; *(entreePrinc->entree) >> nom;
           if (nom != "xn_thermo_dependant_") 
            { cout << "\n erreur en lecture de la thermodependance de xn, on aurait du lire le mot cle xn_thermo_dependant_"
                  << " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
              entreePrinc->MessageBuffer("**erreur6 Loi_maxwell1D::LectureDonneesParticulieres (...**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);     
            };
           // lecture de la loi d'évolution du coefficient non linéaire en fonction de la température
           *(entreePrinc->entree) >>  nom;    
           // on regarde si la courbe existe, si oui on récupère la référence 
           if (lesCourbes1D.Existe(nom))  { xn_temperature = lesCourbes1D.Trouve(nom);}
           else { // sinon il faut la lire maintenant
                  string non_courbe("_");   
                  xn_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
                  // lecture de la courbe
                  xn_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
                };
           // prepa du flot de lecture
           if(strstr(entreePrinc->tablcar,"fin_coeff_MAXWELL1D")==0)
               entreePrinc->NouvelleDonnee();
         }
        else
         { // lecture du coeff 
           *(entreePrinc->entree) >> xn ;
         }; 
      };    
    // prepa du flot de lecture
    if(strstr(entreePrinc->tablcar,"fin_coeff_MAXWELL1D")==0) entreePrinc->NouvelleDonnee(); 
    // appel au niveau de la classe mère
    Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
                 (*entreePrinc,lesFonctionsnD);

  };
// affichage de la loi
void Loi_maxwell1D::Affiche() const 
  { cout << " \n loi de comportement maxwell 1D ";
    if ( E_temperature != NULL) { cout << " module d'young thermo dependant " 
                                       << " courbe E=f(T): " << E_temperature->NomCourbe() <<" ";}
    else                        { cout << " module d'young  E = " << E ;}
    if ( mu_temperature != NULL){ cout << " viscosite thermo dependant " 
                                       << " courbe mu=f(T): " << mu_temperature->NomCourbe() <<" ";}
    else                        { cout << " viscosite mu= " << mu ;} 
    switch (type_derive)
     { case -1: cout << ", et derivee de jauman pour la contrainte" << endl;break;
       case  0: cout << ", et derivee de Liedeux fois covariantes pour la contrainte" << endl;     break;
       case 1:  cout << ", et derivee de Liedeux fois contravariantes pour la contrainte" << endl; break;
     };
    if (!simple) 
     { if (xn_temperature != NULL) { cout << " coef non lineaire thermo dependant " 
                                          << " courbe xn=f(T): " << xn_temperature->NomCourbe() <<" ";}
       else                        { cout << " coef non lin xn= " << xn ;} 
     };
    cout << endl;
    // appel de la classe mère
    Loi_comp_abstraite::Affiche_don_classe_abstraite();
  };
            
// affichage et definition interactive des commandes particulières à chaques lois
void Loi_maxwell1D::Info_commande_LoisDeComp(UtilLecture& entreePrinc)
 {  ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
	   cout << "\n definition standart (rep o) ou exemples exhaustifs (rep n'importe quoi) ? ";
	   string rep = "_";
    // procédure de lecture avec prise en charge d'un retour chariot
    rep = lect_return_defaut(true,"o");
    if (E == -ConstMath::trespetit) 
      { // on initialise à une valeur arbitraire
        E = 110000;}
    if (mu == -ConstMath::trespetit) 
      { // on initialise à une valeur arbitraire
        mu = 0.15; xn = 1.1;}
    sort << "\n# .......  loi de comportement maxwell 1D .................................................................."
         << "\n# |  module d'Young | viscosite      | type de derivee objective utilisee |et eventuellement une puissance |"
         << "\n# |                 |                |  pour le calcul de la contrainte   | pour une evolution non lineaire|"
         << "\n# |       E         |mu (obligatoire)|  type_derivee  (facultatif)        |            xn (facultatif)     |"
         << "\n#..........................................................................................................."
         << "\n    E= "<< setprecision(8) << E << "    mu=       " << setprecision(8) << mu 
         << " type_derivee  " << type_derive << " xnu= " << xn   
         << "\n                fin_coeff_MAXWELL1D         " << endl;
	   if ((rep != "o") && (rep != "O" ) && (rep != "0") )
      { sort << "\n# le type de derivee est optionnel: = -1 -> derivee de jauman, (valeur par defaut) "
         << "\n# = 0 -> derivee deux fois covariantes , "
         << "\n# = 1 -> derivee deux fois contravariantes"
         << "\n# dans le cas ou l'on veut une valeur differente de la valeur par defaut il faut mettre le mot cle"
         << "\n# <type_deriveee> suivi de la valeur -1 ou 0 ou 1" 
         << "\n# \n# chaque parametre peut etre remplace par une fonction dependante de la temperature " 
         << "\n# pour ce faire on utilise un mot cle puis une nom de courbe ou la courbe directement comme avec "
         << "\n# les autre loi de comportement "
         << "\n# exemple pour le module d'young:             E=   E_thermo_dependant_   courbe1  "
         << "\n# exemple pour la viscosite:                  mu=   mu_thermo_dependant_   courbe2  "
         << "\n# exemple pour la coef non lineaire:          xnu=   xnu_thermo_dependant_   courbe3  "
         << "\n# IMPORTANT: a chaque fois qu'il y a une thermodependence, il faut passer une ligne apres la description"
         << "\n# de la grandeur thermodependante, mais pas de passage à la ligne si se n'est pas thermo dependant "
         << "\n#   la derniere ligne doit contenir uniquement le mot cle:     fin_coeff_MAXWELL1D "
         << endl;
       };
    // appel de la classe mère
    Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc);     
  };  		  	  

// test si la loi est complete
int Loi_maxwell1D::TestComplet()
    { int ret = LoiAbstraiteGeneral::TestComplet();
      if (E_temperature == NULL)
       {if (E == -ConstMath::trespetit) 
        { cout << " \n le module d'young n'est pas defini pour  la loi " << Nom_comp(id_comp)
               << '\n';
          ret = 0;
        } 
       }                       
      if ((mu_temperature == NULL) && (mu == -ConstMath::trespetit)) 
       { cout << " \n la viscosite n'est pas defini pour  la loi " << Nom_comp(id_comp)
              << '\n';
         ret = 0;
       }                       
      return ret; 
    }; 
          
// calcul d'un module d'young équivalent à la loi
double Loi_maxwell1D::Module_young_equivalent(Enum_dure temps,const Deformation & def,SaveResul * )
 { double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
   if (thermo_dependant) temperature_tdt = def.DonneeInterpoleeScalaire(TEMP,temps);
   if (E_temperature != NULL) E = E_temperature->Valeur(temperature_tdt);
   if (mu_temperature != NULL) mu = mu_temperature->Valeur(temperature_tdt);
   if (xn_temperature != NULL) xn = xn_temperature->Valeur(temperature_tdt);
   return (E*mu/(mu+E*deltat));
  };
  
//----- 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 Loi_maxwell1D::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D
                                             ,LesFonctions_nD& lesFonctionsnD)
  { string nom; 
	if (cas == 1)
	  { ent >> nom; 
	    if (nom != "maxwell1D")
	     { cout << "\n erreur en lecture de la loi : maxwell1D, on attendait le mot cle : maxwell1D "
	            << "\n Loi_maxwell1D::Lecture_base_info_loi(...";
	       Sortie(1);
	      }
	    // ensuite normalement il n'y a pas de pb de lecture puisque c'est écrit automatiquement (sauf en debug)        
	    ent >> nom >> type_derive;
	    ent >> nom; // module_d_young
	    bool test; ent >> test;
	    if (!test)
	     { ent >> E;
	       if (E_temperature != NULL) {if (E_temperature->NomCourbe() == "_") delete E_temperature; E_temperature = NULL;};
	      }
	    else
	     { ent >> nom; E_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,E_temperature); };
	    // viscosité
	    ent >> nom >> test;
	    if (!test)
	     { ent >> mu;
	       if (mu_temperature != NULL) {if (mu_temperature->NomCourbe() == "_") delete mu_temperature; mu_temperature = NULL;};
	      }
	    else
	     { ent >> nom; mu_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,mu_temperature); };
	    // le coef non linéaire
	    ent >> nom >> nom >> simple >> test;  
	    if (!test)
	     { ent >> xn;
	       if (xn_temperature != NULL) {if (xn_temperature->NomCourbe() == "_") delete xn_temperature; xn_temperature = NULL;};
	      }
	    else
	     { ent >> nom; xn_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,xn_temperature); };
	   } 
	  Loi_comp_abstraite::Lecture_don_base_info(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD); 
  };
	     
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void Loi_maxwell1D::Ecriture_base_info_loi(ofstream& sort,const int cas)
{ if (cas == 1)
   { sort << " maxwell1D " << " type_derivee " << type_derive << " ";
     sort << "\n  module_d_young ";
     if (E_temperature == NULL)
      { sort << false << " " << E << " ";}
     else 
      { sort << true << " fonction_E_temperature ";
        LesCourbes1D::Ecriture_pour_base_info(sort,cas,E_temperature);
       }; 
	    sort << "\n  viscosite ";             
     if (mu_temperature == NULL)
      { sort << false << " " << mu << " ";}
     else 
      { sort << true << " fonction_mu_temperature ";
        LesCourbes1D::Ecriture_pour_base_info(sort,cas,mu_temperature);
       }; 
	    sort << "\n  non_lineaire " << " simple " << simple << " ";             
     if (xn_temperature == NULL)
      { sort << false << " " << xn << " ";}
     else 
      { sort << true << " fonction_xn_temperature ";
        LesCourbes1D::Ecriture_pour_base_info(sort,cas,xn_temperature);
       }; 
   };
  // appel de la classe mère
  Loi_comp_abstraite::Ecriture_don_base_info(sort,cas);
  };

 // ========== codage des METHODES VIRTUELLES  protegees:================
        // calcul des contraintes a t+dt
void Loi_maxwell1D::Calcul_SigmaHH (TenseurHH& sigHH_t,TenseurBB& DepsBB_,DdlElement & ,
             TenseurBB & ,TenseurHH & ,BaseB& ,BaseH& ,TenseurBB&  ,
             TenseurBB& , TenseurBB& gijBB_ ,
		  	 TenseurHH & gijHH_,Tableau <TenseurBB *>& ,double& ,double& ,
		  	  TenseurHH & sigHH_
		  	,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double&  module_cisaillement
		  	,const Met_abstraite::Expli_t_tdt& )
 {   
   #ifdef MISE_AU_POINT	 	 
    if (DepsBB_.Dimension() != 1)
	    { cout << "\nErreur : la dimension devrait etre 1 !\n";
		     cout << " Loi_maxwell1D::Calcul_SigmaHH\n";
		     Sortie(1);
     };
    #endif
    const Tenseur1BB & DepsBB = *((Tenseur1BB*) &DepsBB_); // passage en dim 1
    const Tenseur1HH & gijHH = *((Tenseur1HH*) &gijHH_); //   "      "  "  "
    const Tenseur1BB & gijBB = *((Tenseur1BB*) &gijBB_); //   "      "  "  "
    Tenseur1HH & sigHH = *((Tenseur1HH*) &sigHH_);  //   "      "  "  "
    Tenseur1HH & sigHH_n = *((Tenseur1HH*) &sigHH_t);  //   "      "  "  "

    // cas de la thermo dépendance, on calcul les grandeurs
    if (E_temperature != NULL) E = E_temperature->Valeur(*temperature);
    if (mu_temperature != NULL) mu = mu_temperature->Valeur(*temperature);
    if (xn_temperature != NULL) xn = xn_temperature->Valeur(*temperature);

    Tenseur1BH  DepsBH =  DepsBB * gijHH;
    // recup de l'incrément de temps
	   double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
    // le calcul de la contrainte correspond à l'intégration d'une équation différencielle
    // D=1/E dSig/dt +sig/mu  
    double co1 = deltat * E * mu /(mu + E * deltat);
    switch (type_derive)
     {case -1:
        {Tenseur1BH sigBH_n = gijBB * sigHH_n;
         Tenseur1BH deltaSigBH = co1 * DepsBB * gijHH - (co1/mu) * sigBH_n;
         // ici comme il y a qu'une seule composante, AHB=ABH on ne calcul donc qu'en BH
         sigHH = gijHH * deltaSigBH  + sigHH_n;
         break;}
      case 0: // cas d'une dérivée de Liedeux fois covariantes
        {Tenseur1BB sigBB_n = gijBB * sigHH_n * gijBB;
         Tenseur1BB deltaSigBB = co1 * DepsBB - (co1/mu) * sigBB_n;
         sigHH = gijHH * deltaSigBB * gijHH + sigHH_n; 
         break;}
      case 1:
        {Tenseur1HH deltaSigHH = co1 * (gijHH * DepsBB * gijHH) - (co1/mu) * sigHH_n;
         sigHH = deltaSigHH  + sigHH_n; 
         break;}      
      };
  // ---- traitement des énergies ----
    energ.Inita(0.); double ener_elas=0.; double ener_visqueux = energ_t.DissipationVisqueuse();
    // on peut raisonner en supposant une décomposition de la vitesse de déformation
    //  sig/E=eps_elas -> E_elas = sig*eps_elas/2=sig^2/2/E
    //  sig/mu = D_visqueux -> E_vis = sig^2/mu * delta t
    Tenseur1BH sigBH = gijBB * sigHH;
    double x_inter = (sigBH && sigBH) ;
    ener_elas = x_inter / (2. * E);
    ener_visqueux += x_inter /mu * deltat;
    // ...  mise à jour de energ
    energ.ChangeEnergieElastique(ener_elas); 
    energ.ChangeDissipationVisqueuse(ener_visqueux);
    
    // on libère les tenseurs intermédiaires        
    LibereTenseur();        
 };
        
        // calcul des contraintes a t+dt et de ses variations  
void Loi_maxwell1D::Calcul_DsigmaHH_tdt (TenseurHH& sigHH_t,TenseurBB& DepsBB_,DdlElement & tab_ddl
            ,BaseB& ,TenseurBB & gijBB_t,TenseurHH & gijHH_t,
            BaseB& ,Tableau <BaseB> & ,BaseH& ,Tableau <BaseH> & ,
            TenseurBB & ,Tableau <TenseurBB *>& d_epsBB,TenseurBB & ,
		    TenseurBB & gijBB_tdt,TenseurHH & gijHH_tdt,
		    Tableau <TenseurBB *>& d_gijBB_tdt,
		  	Tableau <TenseurHH *>& d_gijHH_tdt,double& , double& ,
		  	Vecteur& ,TenseurHH& sigHH_tdt,Tableau <TenseurHH *>& d_sigHH
		  	,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double&  module_cisaillement
		  	,const Met_abstraite::Impli& )
 {
   #ifdef MISE_AU_POINT	 	 
	 if (DepsBB_.Dimension() != 1)
	    { cout << "\nErreur : la dimension devrait etre 1 !\n";
		  cout << " Loi_maxwell1D::Calcul_DsigmaHH_tdt\n";
		  Sortie(1);
		};
	 if (tab_ddl.NbDdl() != d_gijBB_tdt.Taille())
	    { cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_tdt  !\n";
		  cout << " Loi_maxwell1D::Calcul_SDsigmaHH_tdt\n";
		  Sortie(1);
		};
    #endif
    const Tenseur1BB & DepsBB = *((Tenseur1BB*) &DepsBB_); // passage en dim 1
    const Tenseur1HH & gijHH = *((Tenseur1HH*) &gijHH_tdt); //   "      "  "  "
    const Tenseur1BB & gijBB = *((Tenseur1BB*) &gijBB_tdt); //   "      "  "  "
    Tenseur1HH & sigHH = *((Tenseur1HH*) &sigHH_tdt);  //   "      "  "  "
    Tenseur1HH & sigHH_n = *((Tenseur1HH*) &sigHH_t);  //   "      "  "  "

    // cas de la thermo dépendance, on calcul les grandeurs
    if (E_temperature != NULL) E = E_temperature->Valeur(*temperature);
    if (mu_temperature != NULL) mu = mu_temperature->Valeur(*temperature);
    if (xn_temperature != NULL) xn = xn_temperature->Valeur(*temperature);

    Tenseur1BH  DepsBH =  DepsBB * gijHH;
    // recup de l'incrément de temps
	double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
    // le calcul de la contrainte correspond à l'intégration d'une équation différencielle
    // D=1/E dSig/dt +sig/mu  

    double unSurDeltat;    
	if (Abs(deltat) >= ConstMath::trespetit)
     {unSurDeltat = 1./deltat;}
    else
    // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
    // {unSurDeltat = Signe(deltat)*ConstMath::tresgrand;}
     { // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
       if (unSurDeltat < 0)
        { cout << "\n le pas de temps est négatif !! "; };
       unSurDeltat = ConstMath::tresgrand;
      }; 

    int nbddl = d_gijBB_tdt.Taille();
    double co1 = deltat * E * mu /(mu + E * deltat);
    switch (type_derive)
     {case -1:
        {Tenseur1BH sigBH_n = gijBB * sigHH_n;
         Tenseur1BH deltaSigBH = co1 * DepsBB * gijHH - (co1/mu) * sigBH_n;
         Tenseur1HB deltaSigHB(deltaSigBH(1,1));
         
         // ici comme il y a qu'une seule composante, AHB=ABH on ne calcul donc qu'en BH
         sigHH = 0.5*(gijHH * deltaSigBH + deltaSigHB * gijHH) + sigHH_n;
         // maintenant calcul de l'opérateur tangent
         for (int i = 1; i<= nbddl; i++)
           {// on fait de faire uniquement une égalité d'adresse et de ne pas utiliser
            // le constructeur d'ou la profusion d'* et de ()
            Tenseur1HH & dsigHH = *((Tenseur1HH*) (d_sigHH(i)));  // passage en dim 1
            const Tenseur1BB & dgijBB =  *((Tenseur1BB*)(d_gijBB_tdt(i)));  // passage en dim 1
            const Tenseur1HH & dgijHH = *((Tenseur1HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
            const Tenseur1BB & depsBB = *((Tenseur1BB *) (d_epsBB(i))); //    "
            Tenseur1BH d_sigBH_n = dgijBB * sigHH_n;
            Tenseur1BH d_deltaSigBH = co1 * ((unSurDeltat) * depsBB * gijHH + DepsBB * dgijHH) 
                                      - (co1/mu) * d_sigBH_n;
            dsigHH = dgijHH * deltaSigBH  + gijHH * d_deltaSigBH; 
            }  
         break;}
      case 0: // cas d'une dérivée de Lie deux fois covariantes
        {Tenseur1BB sigBB_n = gijBB * sigHH_n * gijBB;
         Tenseur1BB deltaSigBB = co1 * DepsBB - (co1/mu) * sigBB_n;
         sigHH = gijHH * deltaSigBB * gijHH + sigHH_n; 
         // maintenant calcul de l'opérateur tangent
         for (int i = 1; i<= nbddl; i++)
           {// on fait de faire uniquement une égalité d'adresse et de ne pas utiliser
            // le constructeur d'ou la profusion d'* et de ()
            Tenseur1HH & dsigHH = *((Tenseur1HH*) (d_sigHH(i)));  // passage en dim 1
            const Tenseur1BB & dgijBB =  *((Tenseur1BB*)(d_gijBB_tdt(i)));  // passage en dim 1
            const Tenseur1HH & dgijHH = *((Tenseur1HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
            const Tenseur1BB & depsBB = *((Tenseur1BB *) (d_epsBB(i))); //    "
            Tenseur1BB d_sigBB_n = dgijBB * sigHH_n * gijBB + gijBB * sigHH_n * dgijBB;
            Tenseur1BB d_deltaSigBB = (co1 * unSurDeltat) * depsBB - (co1/mu) * d_sigBB_n;
            dsigHH = dgijHH * deltaSigBB * gijHH + gijHH * d_deltaSigBB * gijHH 
                     + gijHH * deltaSigBB * dgijHH ; 
            }  
         break;}
      case 1:
        {Tenseur1HH deltaSigHH = co1 * (gijHH * DepsBB * gijHH) - (co1/mu) * sigHH_n;
         sigHH = deltaSigHH  + sigHH_n; 
         // maintenant calcul de l'opérateur tangent
         for (int i = 1; i<= nbddl; i++)
           {// on fait de faire uniquement une égalité d'adresse et de ne pas utiliser
            // le constructeur d'ou la profusion d'* et de ()
            Tenseur1HH & dsigHH = *((Tenseur1HH*) (d_sigHH(i)));  // passage en dim 1
            const Tenseur1BB & dgijBB =  *((Tenseur1BB*)(d_gijBB_tdt(i)));  // passage en dim 1
            const Tenseur1HH & dgijHH = *((Tenseur1HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
            const Tenseur1BB & depsBB = *((Tenseur1BB *) (d_epsBB(i))); //    "
            dsigHH = co1 * (dgijHH * DepsBB * gijHH + gijHH * (( unSurDeltat)*depsBB) * gijHH + gijHH * DepsBB * dgijHH);
            }  
         break;}      
      };  
  // ---- traitement des énergies ----
    energ.Inita(0.); double ener_elas=0.; double ener_visqueux = energ_t.DissipationVisqueuse();
    // on peut raisonner en supposant une décomposition de la vitesse de déformation
    //  sig/E=eps_elas -> E_elas = sig*eps_elas/2=sig^2/2/E
    //  sig/mu = D_visqueux -> E_vis = sig^2/mu * delta t
    Tenseur1BH sigBH = gijBB * sigHH;
    double x_inter = (sigBH && sigBH) ;
    ener_elas = x_inter / (2. * E);
    ener_visqueux += x_inter /mu * deltat;
    // ...  mise à jour de energ
    energ.ChangeEnergieElastique(ener_elas); 
    energ.ChangeDissipationVisqueuse(ener_visqueux);
    
    // on libère les tenseurs intermédiaires        
    LibereTenseur();        
 };