//************************************************************************
// --- Logiciel MEF++
// ---
// --- Copyright 1996-2002 GIREF, Université Laval
// --- TOUS DROITS RÉSERVÉS
// --- ALL RIGHTS RESERVED
// ---
// --- Ce logiciel est couvert par la Loi sur le droit d'auteur.
// --- L'utilisation, ou la modification de ce logiciel sous toutes ses
// --- formes, que ce soit code source ou code compilé, à des fins
// --- personnelles et non commerciales sont autorisées sans frais pour
// --- autant que la présente notice de copyright ainsi que cette
// --- permission apparaissent dans toutes les copies ainsi que dans la
// --- documentation.
// --- Le GIREF et l'Université Laval ne prétendent en aucune façon que
// --- ce logiciel convient à un emploi quelconque. Celui-ci est
// --- distribué sans aucune garantie implicite ou explicite.
//************************************************************************

//********************************************************************
//
//  Resolution du probleme viscoelastique en utilisant le model de Maxwell generalise:
//
//  C_TRho dT/dt - div(K_T grad(T)) = 0                       (1) température
//
//  C_MRho dM/dt - div(K_M grad(M) + K_MT grad(T) ) = 0       (2) teneur en humidité
//
//  -div(sigma(U)) = F                                        (3) déplacements
//
//  avec possiblement un contact sur la partie "entite_en_contact" du domaine (typiquement
//  de la forme UZ >= 0 sur une partie de la frontiere)
//
//  Les équations (1) et (3) sont optionnelles auxquels cas leur traitement est completement
//  ignorées dans le code. Dans (3) le membre de droite F est optionnel.
//
//  USAGE:<pre> ThermoHygroMecanikNewtonViscoElasticiteMaxwellGeneraliseHPP.[dev,opt]
//              préfixe_des_fichiers_d_entrée [entite_en_contact]
//
//  Exemple:
//
//
//********************************************************************

#include "MEFPPUtil.h"
#include "ChaineCar.h"
#include "ChampCelluleViscoElastiqueMaxwellHPP.h"
#include "ChampContraintesHPP.h"
#include "ChampContraintesHPPMateriauVieillissant.h"
#include "ChampElastiqueViscoElastiqueMaxwellHPP.h"
#include "ChampExponentielleTempsRelaxation.h"
#include "ChampGeometrique.h"
#include "ChampPrecalculeDiscPtsIntegration.h"
#include "ChampPseudoContraintesO4.h"
#include "ChampScalaire.h"
#include "ChampScalContinu.h"
#include "ChampScalAConstant.h"
#include "ChampTensO2SymContinu.h"
#include "ChampTensO2SymLinElem.h"
#include "ChampTensorielO4Sym.h"
#include "ChampTensorielO4SymMin.h"
#include "ChampVect3DContinu.h"
#include "ChampContraintesHPP.h"
#include "ExportGIREF.h"
#include "ExportVU.h"
#include "fonctionsChaineCar.h"
#include "Geometrie.h"
#include "GestionFichierChamps.h"
#include "InitialisationArgvArgc.h"
#include "InitialisationTraitementSignal.h"
#include "ListeConditionsLimites.h"
#include "ListEntitesGeometrique.h"
#include "Maillage.h"
#include "MaillagePeau.h"
#include "OptionsSolveurLinPETSc.h"
#include "PETScInitialisation.h"
#include "PPAfficheValeur.h"
#include "PPCallBack.h"
#include "PPCalculChampaXpbY.h"
#include "PPEvalueChampXYZ.h"
#include "PPResolutionProbleme.h"
#include "PPReinterpole.h"
#include "PPCallBack2Args.h"
#include "PPContraintesEL3D.h"
#include "ProblemeEF.h"
#include "ProblemeEFInstationnaire.h"
#include "PtIntgElem1D.h"
#include "PtIntgHexaedre.h"
#include "PtIntgPrismeTri.h"
#include "PtIntgQuadrangle.h"
#include "PtIntgTetraedre.h"
#include "PtIntgTriangle.h"
#include "SchemaIntg.h"
#include "SolveurInstNlinPETSc.h"
#include "SolveurLinPETSc.h"
#include "SolveurGCPAvecCorpsRigide.h"
#include "SolveurQuasiStatiqueNlinPETSc.h"
#include "TFCelluleViscoElastiqueMaxwellHPPEpsilonLin.h"
#include "TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin.h"
#include "TFDiffusion.h"
#include "TFDiffusionNewton.h"
#include "TFMasseGenScal.h"
#include "TFContraintesInitiales.h"
#include "TFCouplageAvecFctsConnues.h"
#include "TFContactV3DUxNVScal.h"
#include "TFSigmaUEpsilonV3D.h"
#include "TFSourceFVVect3D.h"
#include "VerifieElementsNormalises.h"
#include "traducteurERMsg.h"
#include <fstream>
#include <iomanip>

using namespace std;
struct ProprietePurementElastique {

  ChampTensorielO4Sym* Elasticite;
  bool                 MateriauVieillissant;
  ChampScalaire*       DeltaTempsReduit;
  DReel*               LongueurPasDeTemps;

  ChampScalaire* IndicateurEtat;

  ChampTensorielO4Sym* TenseurS;
  ChampTensorielO2Sym* ContraintesCourantes;
  ChampTensorielO2Sym* ContraintesPrecedentes;

};

struct ProprieteViscoElastique {

  ChampTensorielO4Sym*    Elasticite;
  ChampTensorielO4Sym*    Viscosite;
  bool                    MateriauVieillissant;
  ChampScalaire*          DeltaTempsReduit;
  ChampScalaire*          Multiplicateur_b;
  ChampScalaire*          Multiplicateur_bPrecedent;
  ChampScalaire*          IndicateurEtat;

  ChampTensorielO4Sym*    TenseurS;
  ChampTensorielO4Sym*    TenseurI;
  ChampTensorielO4SymMin* PseudoContraintesCourantes;
  ChampTensorielO4SymMin* PseudoContraintesPrecedentes;

};

template<class PTTypeChampSource,class PTTypeChampDestination>
ERMsg ajouteListeReinterpole(PTTypeChampSource& pChampSource, PTTypeChampDestination& pChampDestination, std::vector<PPTraitementGen*>& pVecteurPPReinterpolation){

  ERMsg lMsg;

  Entier lDimensionDuPTTypeChamp1 = pChampSource.reqNbDDLsParComposanteDeMaillage().reqDimension();

  if (lMsg && (1 != lDimensionDuPTTypeChamp1 && 2 != lDimensionDuPTTypeChamp1 && 3 != lDimensionDuPTTypeChamp1
	       && 6 != lDimensionDuPTTypeChamp1 && 9 != lDimensionDuPTTypeChamp1 && 21 != lDimensionDuPTTypeChamp1
               && 36 != lDimensionDuPTTypeChamp1)) {
    lMsg.ajoute(ERMsg(ERMsg::ERREUR, "ERR_DIMENSION_CHAMP_INCONNU"));
    lMsg.ajoute(pChampSource.reqNom());
  }

  if(lMsg){
    switch (lDimensionDuPTTypeChamp1){

    case 1:
    {
      ChampScalaire *lSource      = dynamic_cast<ChampScalaire*>(&pChampSource);
      ChampScalaire *lDestination = dynamic_cast<ChampScalaire*>(&pChampDestination);
      if(lSource == 0 || lDestination == 0) {
        lMsg.ajoute(ERMsg(ERMsg::ERREUR, "ERR_CHAMP_INCOMPATIBLE_AVEC_1D"));
        if(lSource == 0) lMsg.ajoute(pChampSource.reqNom());
        if(lDestination == 0) lMsg.ajoute(pChampDestination.reqNom());
      }
      if(lMsg){
        PPReinterpole<ChampScalaire>* lPPReinterpole1D = new PPReinterpole<ChampScalaire>(*lSource,*lDestination);
        pVecteurPPReinterpolation.push_back(lPPReinterpole1D);
      }
    }
    break;
    case 2:
    {
      ChampVectoriel2D* lSource      = dynamic_cast<ChampVectoriel2D*>(&pChampSource);
      ChampVectoriel2D* lDestination = dynamic_cast<ChampVectoriel2D*>(&pChampDestination);
      if(lSource == 0 || lDestination == 0) {
        lMsg.ajoute(ERMsg(ERMsg::ERREUR, "ERR_CHAMP_INCOMPATIBLE_AVEC_2D"));
        if(lSource == 0) lMsg.ajoute(pChampSource.reqNom());
        if(lDestination == 0) lMsg.ajoute(pChampDestination.reqNom());
      }
      if(lMsg){
        PPReinterpole<ChampVectoriel2D>* lPPReinterpole2D = new PPReinterpole<ChampVectoriel2D>(*lSource,*lDestination);
        pVecteurPPReinterpolation.push_back(lPPReinterpole2D);
      }
    }
    break;
    case 3:
    {
      ChampVectoriel3D* lSource      = dynamic_cast<ChampVectoriel3D*>(&pChampSource);
      ChampVectoriel3D* lDestination = dynamic_cast<ChampVectoriel3D*>(&pChampDestination);
      if(lSource == 0 || lDestination == 0) {
        lMsg.ajoute(ERMsg(ERMsg::ERREUR, "ERR_CHAMP_INCOMPATIBLE_AVEC_3D"));
        if(lSource == 0) lMsg.ajoute(pChampSource.reqNom());
        if(lDestination == 0) lMsg.ajoute(pChampDestination.reqNom());
      }
      if(lMsg){
        PPReinterpole<ChampVectoriel3D>* lPPReinterpole3D = new PPReinterpole<ChampVectoriel3D>(*lSource,*lDestination);
        pVecteurPPReinterpolation.push_back(lPPReinterpole3D);
      }
    }
    break;
    case 6:
    {
      ChampTensorielO2Sym* lSource      = dynamic_cast<ChampTensorielO2Sym*>(&pChampSource);
      ChampTensorielO2Sym* lDestination = dynamic_cast<ChampTensorielO2Sym*>(&pChampDestination);
      if(lSource == 0 || lDestination == 0) {
        lMsg.ajoute(ERMsg(ERMsg::ERREUR, "ERR_CHAMP_INCOMPATIBLE_AVEC_6D"));
        if(lSource == 0) lMsg.ajoute(pChampSource.reqNom());
        if(lDestination == 0) lMsg.ajoute(pChampDestination.reqNom());
      }
      if(lMsg){
        PPReinterpole<ChampTensorielO2Sym>* lPPReinterpole6D = new PPReinterpole<ChampTensorielO2Sym>(*lSource,*lDestination);
        pVecteurPPReinterpolation.push_back(lPPReinterpole6D);
      }
    }
    break;
    case 9:
    {
      ChampTensorielO2NSym* lSource      = dynamic_cast<ChampTensorielO2NSym*>(&pChampSource);
      ChampTensorielO2NSym* lDestination = dynamic_cast<ChampTensorielO2NSym*>(&pChampDestination);
      if(lSource == 0 || lDestination == 0) {
        lMsg.ajoute(ERMsg(ERMsg::ERREUR, "ERR_CHAMP_INCOMPATIBLE_AVEC_9D"));
        if(lSource == 0) lMsg.ajoute(pChampSource.reqNom());
        if(lDestination == 0) lMsg.ajoute(pChampDestination.reqNom());
      }
      if(lMsg){
        PPReinterpole<ChampTensorielO2NSym>* lPPReinterpole9D = new PPReinterpole<ChampTensorielO2NSym>(*lSource,*lDestination);
        pVecteurPPReinterpolation.push_back(lPPReinterpole9D);
      }
    }
    break;
    case 21:
    {
      ChampTensorielO4Sym* lSource      = dynamic_cast<ChampTensorielO4Sym*>(&pChampSource);
      ChampTensorielO4Sym* lDestination = dynamic_cast<ChampTensorielO4Sym*>(&pChampDestination);
      if(lSource == 0 || lDestination == 0) {
        lMsg.ajoute(ERMsg(ERMsg::ERREUR, "ERR_CHAMP_INCOMPATIBLE_AVEC_21D"));
        if(lSource == 0) lMsg.ajoute(pChampSource.reqNom());
        if(lDestination == 0) lMsg.ajoute(pChampDestination.reqNom());
      }
      if(lMsg){
        PPReinterpole<ChampTensorielO4Sym>* lPPReinterpole21D = new PPReinterpole<ChampTensorielO4Sym>(*lSource,*lDestination);
        pVecteurPPReinterpolation.push_back(lPPReinterpole21D);
      }
    }
    break;
    case 36:
    {
      ChampTensorielO4SymMin* lSource      = dynamic_cast<ChampTensorielO4SymMin*>(&pChampSource);
      ChampTensorielO4SymMin* lDestination = dynamic_cast<ChampTensorielO4SymMin*>(&pChampDestination);
      if(lSource == 0 || lDestination == 0) {
        lMsg.ajoute(ERMsg(ERMsg::ERREUR, "ERR_CHAMP_INCOMPATIBLE_AVEC_36D"));
        if(lSource == 0) lMsg.ajoute(pChampSource.reqNom());
        if(lDestination == 0) lMsg.ajoute(pChampDestination.reqNom());
      }
      if(lMsg){
        PPReinterpole<ChampTensorielO4SymMin>* lPPReinterpole36D = new PPReinterpole<ChampTensorielO4SymMin>(*lSource,*lDestination);
        pVecteurPPReinterpolation.push_back(lPPReinterpole36D);
      }
    }
    break;
    }
  }

  return lMsg;
}


ERMsg lisString(const char*& pCharPtrAModifier,std::string& pStringRetour);

int main(int argc, char *argv[]) {

// ***********************************
// * BLOC: Initialisations générales *
// ***********************************
  PETScInitialisation lPETScInit(&argc, &argv);
  sInitialisationArgvArgc.asgnArguments(argc, argv);
  initialiseTraitementSignal();

  const vector<string> lArgv = sInitialisationArgvArgc.reqArgv();

  // On déclare une variable qui contiendra potentiellement un message d'erreur
  ERMsg lMsg = ERMsg(ERMsg::OK);

  // On recupere le GroupeProcessus associe a PETSC_COMM_WORLD
  const PAGroupeProcessus lGroupeGlobal = PETScInitialisation::reqPetscCommWorld();

  // On déclare les variables qui seront extraites de la ligne de commande
  ChaineCar lNomMaillage, lNomFichierChamp, lNomEntite;
  bool lContakt(false);               // le probleme n'est pas un pb de contact

  // Vérification des paramètres passés à la ligne de commande.
  const Entier lNbArg = lArgv.size();
  if (lNbArg < 2 || lNbArg > 3) {
    ChaineCar lNomProgramme = sInitialisationArgvArgc.reqNomProgramme();
    ChaineCar lMessageUsage = "USAGE: \n";
    lMessageUsage += enleveCheminAccesFichier(lNomProgramme);
    lMessageUsage += " Prefixe_fichiers_champs Nom_entite_de_contact (optionnel)";
    lMsg.ajoute(ERMsg(ERMsg::ERREUR, lMessageUsage));
  }
  else {
    lNomFichierChamp = lArgv[1];
    if(lNbArg>2){
      lNomEntite  = lArgv[2];
      lContakt    = true;
    }
  }

  lNomMaillage = lNomFichierChamp;

  // Pre-traitement des donnees du fichier champs:
  {
    const EntierN lTailleBuffer = 65536;
    char *lLigne = new char[lTailleBuffer];
    std::ifstream lFichierChamp((lNomFichierChamp + std::string(".champs")).c_str());
    bool lNomMaillageTrouve(false);
    bool lNomEntiteTrouve(false);

    if(lFichierChamp) {
      lFichierChamp.getline(lLigne,lTailleBuffer);

      while(lFichierChamp && lMsg && (!lNomMaillageTrouve || !lNomEntiteTrouve)){
	if(retournePremierMotSansEspace(lLigne) == "chainecar"){
	  ChaineCar lNomBidon;
	  const char* lPositionDansLigne = &lLigne[0];
	  lMsg = lisString(lPositionDansLigne, lNomBidon);
	  if(lMsg) lMsg = lisString(lPositionDansLigne, lNomBidon);
	  if(!lNomMaillageTrouve && lNomBidon == "NOM_DU_MAILLAGE") {
	    if(lMsg) {
	      lMsg = lisString(lPositionDansLigne, lNomMaillage);
	    }
	    lNomMaillageTrouve = true;
	  }

	  if( !lNomEntiteTrouve && lNomBidon == "NOM_ENTITE_DE_CONTACT") {
	    if(lMsg) {
	      lMsg = lisString(lPositionDansLigne, lNomEntite);
	      lContakt    = true;
	    }
	    lNomEntiteTrouve = true;
	  }
	}
	lFichierChamp.getline(lLigne,lTailleBuffer);
      }
      //lFichierChamp.seekg (0, std::ios::beg);
    }
    lFichierChamp.close();
  }
  // Fin du pre-traitement

  // *****************************************************
  // BLOC: Création des objets utilisés dans le calcul *
  // *****************************************************
  // On crée: la géométrie, la liste d'entités géométriques, le maillage et le champ géométrique.
  Geometrie              lGeo;
  ListEntitesGeometrique lListeEntite;
  ListEntitesGeometrique lListeEntitePeau;

  Maillage               lMail;
  MaillagePeau           lMailPeau;


  // Creation des champs requis pour le probleme :
  ChampGeometrique       *lChampGeo= 0;        // champ geometrique du domnaine entier

  ChampScalContinu       *lM       = 0;        // humidite
  ChampScalContinu       *lMP      = 0;        // humidite au pas precedent (necessaire en cas d'hysterese seulement)
  ChampScalContinu       *lMt      = 0;        // dérivée en temps de M
  ChampScalContinu        *lC_MRho = 0;        // densite * "chaleur specifique"
  ChampTensO2SymContinu   *lK_M    = 0;        // diffusivité effective d_b*D
  ChampTensO2SymContinu   *ldK_MdM = 0;        // derivee de la diffusivite effective en M

 //  ChampScalContinu       *lMtB      = 0;        // dérivée en temps de M
//   ChampScalContinu        *lC_MRhoB = 0;        // densite * "chaleur specifique"
//   ChampTensO2SymContinu   *lK_MB    = 0;        // diffusivité effective d_b*D
//   ChampTensO2SymContinu   *ldK_MdMB = 0;        // derivee de la diffusivite effective en M

  ChampTensO2SymContinu   *lK_MT   = 0;        // thermo-diffusion
  ChampTensO2SymContinu   *lK_T    = 0;        // diffusion thermique
  ChampScalContinu        *lC_TRho = 0;        // chaleur spécifique multiplier par la densité
  ChampScalContinu          *lT    = 0;        // temperature

  ChampVect3DContinu        *lU    = 0;        // les déplacements
  ChampVect3DContinu       *lUPrec = 0;        // les déplacements
  ChampVect3DContinu*lDeplacements = 0;        // les déplacements

  ChampTensO2SymContinu    *lBetaM = 0;        // tenseur O2 des dilatation hydrique
  ChampTensO2SymContinu    *lBetaMB = 0;        // tenseur O2 des dilatation hydrique
  ChampTensO2SymContinu    *lBetaT = 0;        // tenseur O2 des dilatation thermique
  ChampVectoriel3D         *lF     = 0 ;       // les charges volumiques

  // Creation des champs requis pour le probleme Peau
  ChampGeometrique        *lChampGeoPeau=0;      // champ geometrique de la peau
  ChampVect3DContinu      *lUPeau = 0;           // Une version "peau" de U
  ChampScalaire           *lLambdaCtc = 0;       // Ce champ contiendra le multiplicateur
  ChampVectoriel3D        *lNormaleExtPeau = 0;  // Ce champ doit contenir la normale extérieure
  ChampScalContinu        *lGapPeau = 0;         // Ce champ pour Bu <= g.
  ChampScalaire           *lMPeau = 0;           // Une version "peau" de M
  ChampScalaire           *lTPeau = 0;           // Une version "peau" de T

  // Variable lues dans le fichier:

  bool   lThermik(false);               // le problème a une composante thermique
  bool   lMecanik(true);                // le problème a  une composante mecanique
  bool   lEmploiNewtonEnM(false);       // On resout avec le terme dK_M/DM (methode de Newton au lieu de pt fixe)#
  bool   lHysterese(true);              // le tenseur de dilatation hydrique tiens compte d'effet d'hysterese
  Entier lPasDeCalcul;                  // frequence du calcul de mecanique: 0 -> pas de calcul
  Entier lPasImpression;                // frequence d'impression des résultats
  Entier lNbTempsImpose;                // Nombre de points du temps ou on force le calcul
  SolveurInstNlinPETSc  *lSolveurHyd;   // le solveur instationnaire pour M (et T possiblement)

  SolveurQuasiStatiqueNlinPETSc * lSolveurQuasiStatique = 0;  // le solveur quasi statique
  SolveurLinPETSc               *lSolveurMec;                // le solveur stationnaire pour U
  DReel           *lPointeurVersLePasDeTempsDuSolveur = 0;     // Un Pointeur Vers Le Pas De Temps Du Solveur

  ChaineCar lNomEntiteControle;          // le nom de l'entite servant au controle des mouvements de corps rigide
  Entier    lDirectionDesCorps;          // Numero de l'axe paralelle aux corps rigide

  VectDyn<ProprieteViscoElastique> lVectProprieteBrancheViscoElastique;  // un vecteur contenant les proprietes viscoelastique pour chaque branche

  std::vector<PPTraitementGen*> lVectChampAReinterpoler;   // un vecteur des champs qu'on devra reinterpoler AVANT la resolution

  std::vector<Champ*> lVectChampADeleter; // un vecteur des champs qui devront etre detruit en fin de programme

  // Variable lues dans le fichier:

  Entier lNombreBrancheViscoElastique = 0;  // le nombre de composantes/cellules/branches viscoélastique du problème

  bool   lReassemble(true);               // reassemblage possible de la matrice de rigidite.

  VectDyn<DReel> lTempsImposes;           // les points dans le temps ou on veut imposer une approximation des champs
  VectDyn<Vecteur3D> lCoordonnees;        // les coordonnées dont on fait une "trace" pendant l'execution

  bool lVrai(true);

  Entier lPoints1D(4);

  SchemaIntg     lSchemaIntg;
  SchemaIntg     lSchemaIntgNewtonCote;

  if (lMsg) {
    // Création du schéma d'intégration des termes de formulation

    lSchemaIntg.asgnIntegration(PtIntgTriangle(PtIntgTriangle::Hammer_DouzePtsInternes_Degre6));

    lSchemaIntg.asgnIntegrationFace3Sommets(PtIntgTriangle(PtIntgTriangle::Hammer_DouzePtsInternes_Degre6));
    lSchemaIntg.asgnIntegrationFace4Sommets(PtIntgQuadrangle(PtIntgElem1D(lPoints1D),PtIntgElem1D(lPoints1D)));

    lSchemaIntg.asgnIntegration(PtIntgQuadrangle(PtIntgElem1D(lPoints1D),PtIntgElem1D(lPoints1D)));
    lSchemaIntg.asgnIntegration(PtIntgElem1D(lPoints1D));
    lSchemaIntg.asgnIntegrationArete(PtIntgElem1D(lPoints1D));

    lSchemaIntg.asgnIntegration(PtIntgHexaedre(PtIntgElem1D(lPoints1D), PtIntgElem1D(lPoints1D), PtIntgElem1D(lPoints1D)));
    lSchemaIntg.asgnIntegration(PtIntgPrismeTri(PtIntgTriangle(PtIntgTriangle::Hammer_DouzePtsInternes_Degre6), PtIntgElem1D(lPoints1D)));
    lSchemaIntg.asgnIntegration(PtIntgTetraedre(5));

    // On crée un schéma pour la résolution du problème de contact si nécessaire
    if(lMsg && lContakt) {

      PtIntgElem1D      lPtIntgElem1D  (PtIntgElem1D::Newton_Cotes_DeuxPoints_Degre1);
      PtIntgTriangle    lPtIntgTriangle(PtIntgTriangle::Newton_Cotes_TroisPoints_Degre1);
      PtIntgQuadrangle  lPtIntgQuad    (lPtIntgElem1D,lPtIntgElem1D);
      PtIntgTetraedre   lPtIntgTetra   (PtIntgTetraedre::Hammer_QuinzePtsInternes_Degre5);
      PtIntgHexaedre    lPtIntgHexa    (PtIntgElem1D(2),PtIntgElem1D(2),PtIntgElem1D(2));

      lSchemaIntgNewtonCote.asgnIntegration(lPtIntgElem1D);
      lSchemaIntgNewtonCote.asgnIntegration(lPtIntgTriangle);
      lSchemaIntgNewtonCote.asgnIntegration(lPtIntgQuad);
      lSchemaIntgNewtonCote.asgnIntegration(lPtIntgTetra);

      lSchemaIntgNewtonCote.asgnIntegration(lPtIntgHexa);
      lSchemaIntgNewtonCote.asgnIntegrationArete(lPtIntgElem1D);
      lSchemaIntgNewtonCote.asgnIntegrationFace3Sommets(lPtIntgTriangle);
      lSchemaIntgNewtonCote.asgnIntegrationFace4Sommets(lPtIntgQuad);
    }

  }

  //On déclare les problèmes "vides"
  ProblemeEFInstationnaire lProbHyd;
  lProbHyd.asgnNom("Composante hydrique/thermique");
  ProblemeEF lProbMec;
  lProbMec.asgnNom("Composante mecanique");
  ProblemeEF lProbPeauGCP;
  ProblemeEF lProbPeauGCP2;

  lProbHyd.asgnGroupeProcessus(lGroupeGlobal);
  lProbMec.asgnGroupeProcessus(lGroupeGlobal);
  lProbPeauGCP.asgnGroupeProcessus(lGroupeGlobal);
  lProbPeauGCP2.asgnGroupeProcessus(lGroupeGlobal);

  if(lContakt && lMsg){
    lMsg = lProbHyd.lisDonnees(lNomMaillage,lMail,lMailPeau,lNomEntite,lGeo,lListeEntite,lListeEntitePeau);
  } else if(lMsg) {
    lMsg = lProbHyd.lisDonnees(lNomMaillage,lMail, lGeo,lListeEntite);
  }

  VerifieElementsNormalises lVerifieElementsNormalises;
  if(!lVerifieElementsNormalises(lMail,lVrai)) lMsg.ajoute(ERMsg(ERMsg::ERREUR,"Le maillage DOIT etre normalise!"));

  //Déclaration de l'objet qui fait la gestion des champs pour le problème.
  GestionFichierChamps lFichierChamp;
  GestionFichierChamps lFichierChampPeau;

  if (lMsg) {
    lFichierChamp.asgnDonneesBase(&lMail);
    lFichierChamp.asgnDonneesGeometriques(&lGeo,&lListeEntite);
  }

  //On demande au gestionnaire de lire le fichier de champ
  if(lMsg) lMsg = lFichierChamp.lireFichier((lNomFichierChamp + std::string(".champs")).c_str());

  // On choisit s'il faut reassembler la matrice de rigidite a chaque pas de temps
  if(lMsg) lMsg = lFichierChamp.reqBooleen("Reassemble",lReassemble,true);

  //
  // Le solveur quasi-statique, possiblement on voudra un acces au pas de temps du solveur ce
  // qui nous oblige a définir le solveur ici.
  //
  bool lDefinitionDansFichierChamps(true);
  if(lMsg) lMsg = lFichierChamp.reqSolveur<SolveurQuasiStatiqueNlinPETSc>("SolveurQS",lSolveurQuasiStatique,false);
  if(lSolveurQuasiStatique==0) {
    lSolveurQuasiStatique = new SolveurQuasiStatiqueNlinPETSc();
    lDefinitionDansFichierChamps=false;
  }

  // Requête des données du fichier .champs nécessaires au calcul
  // On determine le type de probleme:

  if(lMsg) lMsg = lFichierChamp.reqBooleen("Thermique",lThermik,false);
  if(lMsg) lMsg = lFichierChamp.reqBooleen("Mecanique",lMecanik,false);

  if(lMsg && lContakt && !lMecanik) {
    std::cout<<"Pas de contact sans mecanique! On remet le contact à faux"<<std::endl;
    lContakt = false;
  }

  // On demande le champ de transformation géométrique au fichier de champs
  if (lMsg) {
    std::cout<< "Entre dans le fichier champ "<<std::endl;
    lMsg = lFichierChamp.reqChamp<ChampGeometrique>("ChampGeo",lChampGeo);
    std::cout<< "Sort du fichier champ "<<std::endl;
  }
  if(lMsg) lProbHyd.asgnChampGeometrique(*lChampGeo);

  //if(lMecanik && lMsg) lProbMec.asgnDonnees(lMail,lGeo,*lChampGeo,lListeEntite);

  // ***************************
  // * Pour la partie hydrique:
  // ***************************

  //On demande  un champ qui se nomme "M" pour l'humidité
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("M",lM);

  //On demande  un champ qui se nomme "M_t" pour la valeur difference finie de M (hystérèse présente dans le problème)
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("M_t",lMt,false);
  // if(lMsg){
//     lMsg = lFichierChamp.reqChamp<ChampScalContinu>("M_tB",lMtB,false);
//      std::cout<< "fffffffffM_tB "<<std::endl;
//   }

  //On demande  un champ qui se nomme "d_b" pour la densité de base divisée par 100 (kg m^3 / 100)
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("C_MRho",lC_MRho);
  // if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("C_MRhoB",lC_MRhoB);
   //if(lMsg && lC_MRhoB == 0) lC_MRhoB = lC_MRho->newChampCopie();

  //On demande  un champ qui se nomme "K_M" pour le tenseur de "diffusivité effective" (d_b * D)
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("K_M",lK_M);
  //if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("K_MB",lK_MB);
  //if(lMsg && lK_MB == 0) lK_MB = lK_M->newChampCopie();

  //On demande  un champ qui se nomme "dK_MdM" pour la derivee du tenseur de "diffusivité effective" en M
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("dK_MdM",ldK_MdM,false);
  // if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("dK_MdMB",ldK_MdMB,false);
   //if(lMsg && ldK_MdMB == 0) ldK_MdMB = ldK_MdM->newChampCopie();
  if(lMsg && ldK_MdM) {
    lEmploiNewtonEnM = true;
     std::cout<< "Fin Pour la partie hydrique "<<std::endl;}


  // ***************************
  // * Pour la partie thermique
  // ***************************

  if(lThermik && lMsg){
    //On demande  un champ qui se nomme "K_MT" pour la composante thermique de la "diffusion hydrique"
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("K_MT",lK_MT);

    //On demande  un champ qui se nomme "T" pour la température
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("T",lT);

    //On demande  un champ qui se nomme "K_T" pour la "diffusion thermique"
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("K_T",lK_T);

    //On demande  un champ qui se nomme "C" pour la chaleur spécifique multiplie par la densité de base
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("C_TRho",lC_TRho);
  }

  // ***************************
  // * pour la partie mécanique:
  // ***************************

  // Frequence du calcul mecanique
  lPasDeCalcul = 1;
  if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("freqCalcul",lPasDeCalcul,false);

  if(lMsg && lPasDeCalcul == 0){
    std::cout<<"ATTENTION: pas de sauvegarde des resultats de mecanique/contact alors je ne les fait pas!"<<std::endl;
    lMecanik = false;
  }

  ChampScalAConstant* lUn = 0;
  if(lMecanik && lMsg){
    // On demande  un champ qui se nomme "U" pour le champ des deplacements.
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampVect3DContinu>("U",lU);
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampVect3DContinu>("Deplacements",lDeplacements,false);
    if(!lDeplacements && lContakt)
      lMsg.ajoute(ERMsg(ERMsg::ERREUR, "En cas de contact on doit definir Deplacements dans le fichier champs et dans celui de peau."));

    if(lMsg && !lDeplacements) lDeplacements = lU->newChampCopie();

    // On demande  un champ qui se nomme "UPrec" pour le deplacement au temps precedent.
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampVect3DContinu>("UPrec",lUPrec,false);
    if(lMsg && lUPrec == 0) lUPrec = lU->newChampCopie();

    // On demande  un champ qui se nomme "F" pour le chargement volumique
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampVectoriel3D>("F",lF,false);

    // On demande un champs qui se nomme "Beta" pour le tenseur de dilatation hydrique
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("BetaM",lBetaM);
    //if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("BetaMB",lBetaMB);


    // Y'a-t'il une hysterese dans la description du tenseur BetaM?
    if(lMsg && lBetaM){
      lMsg = lFichierChamp.reqBooleen("Hysterese",lHysterese,false);
      if(lHysterese) lFichierChamp.reqChamp<ChampScalContinu>("MP",lMP);
    }

    // On demande un champs qui se nomme "BetaT" pour le tenseur de dilatation thermique
    if(lMsg && lThermik) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("BetaT",lBetaT,false);
    if(lMsg && lBetaT){
      lUn = new ChampScalAConstant(*lChampGeo,1.0, "Un");
    }

    if(lContakt && lMsg){
      // Champs sur la peau

      lFichierChampPeau.asgnDonneesBase(&lMailPeau);
      lFichierChampPeau.asgnDonneesGeometriques(&lGeo,&lListeEntitePeau);

      if (lMsg) {
	//On crée un champ d'extension de peau pour tous les champs du maillage d'origine.
	lMsg = lFichierChamp.creeTousChampsExtensionPeau(lFichierChamp,
							 lFichierChampPeau,
							 lMailPeau);
      }

      if (lMsg) {
	lMsg = lFichierChampPeau.lireFichier(lNomFichierChamp + std::string("-peau.champs"));
      }

      // On demande le champ de transformation géométrique au fichier de champs Peau
      if (lMsg) lMsg = lFichierChampPeau.reqChamp<ChampGeometrique>("ChampGeo",lChampGeoPeau);

      //       // Cree une extension de U sur la peau

      if(lMsg) lMsg = lFichierChampPeau.reqChamp("U",lUPeau);
      if(lMsg) lMsg = lFichierChampPeau.reqChamp("M",lMPeau);
      if(lMsg && lThermik) lMsg = lFichierChampPeau.reqChamp("T",lTPeau);
      if(lMsg) lMsg = lFichierChampPeau.reqChamp("LambdaCtc",lLambdaCtc);
      if(lMsg) lMsg = lFichierChampPeau.reqChamp("normale_ext",lNormaleExtPeau);
      if(lMsg) lMsg = lFichierChampPeau.reqChamp("GapPeau",lGapPeau);

      if(lMsg) lProbPeauGCP.asgnDonnees( lMailPeau, lGeo, *lChampGeoPeau, lListeEntitePeau);
      if(lMsg) lProbPeauGCP2.asgnDonnees( lMailPeau, lGeo, *lChampGeoPeau, lListeEntitePeau);
    }

    std::cout<< "Fin Pour la partie Mecanique "<<std::endl;
  }

  bool lCorpsParGCP(lContakt);    // SolveurGCP s'occupe du corps rigide

  if(lMsg && lCorpsParGCP){
    lCorpsParGCP = false;
    lMsg = lFichierChamp.reqBooleen("CorpsRigideGCP",lCorpsParGCP,false);
  }

  if(lMsg && lCorpsParGCP){
    lDirectionDesCorps = 2;
    lMsg = lFichierChamp.reqScalaireConstant("DirectionRigide",lDirectionDesCorps);//lReelBidon);
    lMsg = lFichierChamp.reqChaineCar("lNomEntiteControle",lNomEntiteControle);
    //    lDirectionDesCorps = (Entier) lReelBidon;
     std::cout<< "lMsg && lCorpsParGCP "<<std::endl;
  }

  //
  // Composante purement elastique: ce qui rend le modele "generalise", si absent du
  // fichier champ on suppose que l'on utilise un modele de Maxwell
  //

  ProprietePurementElastique                              lProprietePurementElastique = {0};

  ChampTensorielO4Sym*                                    lERef_Inf             = 0;
  bool                                                    lMateriauVieillissant = false;
  ChampScalaire*                                          lDeltaTempsReduit     = 0;
  ChampScalaire*                                          lIndicateurEtat       = 0;
  ChampElastiqueViscoElastiqueMaxwellHPP*                 lChampTenseurS        = 0;
  ChampContraintesHPPMateriauVieillissant*                lChampSigmaCourant    = 0;
  ChampPrecalculeDiscPtsIntegration<ChampTensorielO2Sym>* lChampSigmaPrecedent  = 0;

  //On demande un champ qui se nomme "ERef_Inf" pour le tenseur d'élasticité de la composante purement elastique
  //On suppose que E_Inf = b(t)*ERef_Inf
  if(lMsg) lMsg =lFichierChamp.reqChamp<ChampTensorielO4Sym>("ERef_Infijkl", lERef_Inf,false);

  if(lMsg && lERef_Inf){// si ERef_Inf
    std::cout<< "lMsg && lERef_Inf "<<std::endl;

    //On demande si l'elasticite depend du temps
    if(lMsg) lMsg =lFichierChamp.reqBooleen("E_InfVieillissant",lMateriauVieillissant);

    //On demande un champ pour le pseudo temps reduit: l'integrale (sur un pas de temps) du multiplicateur b
    if(lMsg && lMateriauVieillissant) lMsg =lFichierChamp.reqChamp<ChampScalaire>("Integrale_bInf", lDeltaTempsReduit);

    //On demande un champ pour l'état du solide dans la branche courante
    if(lMsg) lMsg =lFichierChamp.reqChamp<ChampScalaire>("IndicateurEtat_Inf", lIndicateurEtat);

    if(lMsg){
      //
      //On construit les tenseurs qui sont vraiement utilisés dans le modèle
      //
      lChampTenseurS       = new ChampElastiqueViscoElastiqueMaxwellHPP;
      lChampSigmaCourant   = new ChampContraintesHPPMateriauVieillissant;
      lChampSigmaPrecedent = new ChampPrecalculeDiscPtsIntegration<ChampTensorielO2Sym>;

      lVectChampADeleter.push_back(lChampTenseurS);
      lVectChampADeleter.push_back(lChampSigmaCourant);
      lVectChampADeleter.push_back(lChampSigmaPrecedent);
    }

    if(lMsg && lMateriauVieillissant) lMsg = ajouteListeReinterpole(*lChampSigmaCourant,*lChampSigmaPrecedent,lVectChampAReinterpoler);

    // On assigne le champ geometrique
    if(lMsg) {
      lChampTenseurS      ->asgnChampGeometrique(*lChampGeo);
      lChampSigmaCourant  ->asgnChampGeometrique(*lChampGeo);
      lChampSigmaPrecedent->asgnChampGeometrique(*lChampGeo);
    }

    if(lMsg) lSolveurQuasiStatique->reqPtrDelta(lPointeurVersLePasDeTempsDuSolveur);

    if(lPointeurVersLePasDeTempsDuSolveur == 0){
      lMsg.ajoute(ERMsg(ERMsg::ERREUR,"ERREUR Le pointeur vers le pas de temps du solveur QS est nul"));
    }

    // On assigne les parametres aux differents champs tensoriels

    if (lMsg) {

      if(lMateriauVieillissant){
        lChampSigmaPrecedent                    ->asgnChamp(*lChampSigmaCourant);

        lChampTenseurS->asgnChampsAvecVieillissement(*lERef_Inf,
                                                     *lDeltaTempsReduit,
                                                     lPointeurVersLePasDeTempsDuSolveur);
      }

      if(!lMateriauVieillissant){
        lChampTenseurS->asgnChampsSansVieillissement(*lERef_Inf);
        lChampSigmaPrecedent    ->asgnChamp(*lChampSigmaCourant);
      }


      lChampSigmaCourant                       ->asgnChamps(*lU,
                                                            *lUPrec,
                                                            *lChampTenseurS,
                                                            *lChampSigmaPrecedent,
                                                            *lIndicateurEtat);
      std::cout<< "Fin On assigne les parametres aux differents champs  "<<std::endl;
    }

    if(lMsg) lChampSigmaPrecedent->ajouteSchemaIntegration(lSchemaIntg);

  }

  if(lMsg) {
    lProprietePurementElastique.Elasticite             = lERef_Inf;
    lProprietePurementElastique.MateriauVieillissant   = lMateriauVieillissant;
    lProprietePurementElastique.DeltaTempsReduit       = lDeltaTempsReduit;
    lProprietePurementElastique.LongueurPasDeTemps     = lPointeurVersLePasDeTempsDuSolveur;
    lProprietePurementElastique.IndicateurEtat         = lIndicateurEtat;

    lProprietePurementElastique.TenseurS               = lChampTenseurS;
    lProprietePurementElastique.ContraintesCourantes   = lChampSigmaCourant;
    lProprietePurementElastique.ContraintesPrecedentes = lChampSigmaPrecedent;
  }

  //
  // On passe au traitement des branche/cellules/composantes viscoelastique
  //
  // les donnees relatives a chaque branche sont données dans le fichier champs de la facon suivante:
  //
  //      NombreBrancheViscoElastique <= 100
  //
  //      {  Branche_I
  //         ERef_ijkl
  //         ViscoRef_ijkl              <- I = 00,01,02,...,NombreBrancheViscoElastique-1
  //         MateriauVieillissant
  //         DeltaTempsReduit
  //         Multiplicateur_b          (est necessaire SEULEMENT si MateriauVieillissant = vrai)
  //         Multiplicateur_bPrecedent (optionnel si absent et necessaire sera cree avec newChampCopie)
  //         IndicateurEtat            (est necessaire SEULEMENT si MateriauVieillissant = vrai)
  //      }


  if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("NombreBrancheViscoElastique",lNombreBrancheViscoElastique,false);

  if(lMsg && lERef_Inf == 0 && lNombreBrancheViscoElastique == 0) {
    lMsg.ajoute(ERMsg(ERMsg::ERREUR, "Modele de Maxwell sans composante viscoelastique ????"));
  }

  if(lMsg && lNombreBrancheViscoElastique) {
    lVectProprieteBrancheViscoElastique.asgnLongueur(lNombreBrancheViscoElastique);
  }

  ChaineCar lNomBranche= "Branche_";
  for(Entier lNumeroDeBranche=0;lNumeroDeBranche<lNombreBrancheViscoElastique && lMsg;++lNumeroDeBranche){
    // on cree le nom du namescape et les noms des propriétés exigées pour chaque branches
    ChaineCar lNomBranche_i         = ajouteNumeroALaFin(lNomBranche,lNumeroDeBranche,2);

    ChaineCar lNomEijkl             = lNomBranche_i + "::ERef_ijkl";
    ChaineCar lNomViscoijkl         = lNomBranche_i + "::ViscoRef_ijkl";
    ChaineCar lNomVieillissant      = lNomBranche_i + "::MateriauVieillissant";
    ChaineCar lNomb                 = lNomBranche_i + "::Multiplicateur_b";
    ChaineCar lNombPrecedent        = lNomBranche_i + "::Multiplicateur_bPrecedent";
    ChaineCar lNomDeltaTempsReduit  = lNomBranche_i + "::DeltaTempsReduit";
    ChaineCar lNomEtat              = lNomBranche_i + "::IndicateurEtat";

    //On demande un champ pour le tenseur des coefficients d'élasticité de reference
    ChampTensorielO4Sym* lElasticite=0;
    if(lMsg) lMsg =lFichierChamp.reqChamp<ChampTensorielO4Sym>(lNomEijkl,lElasticite);

    //On demande un champ pour le tenseur des coefficients de viscosite de reference
    ChampTensorielO4Sym* lViscosite=0;
    if(lMsg) lMsg =lFichierChamp.reqChamp<ChampTensorielO4Sym>(lNomViscoijkl,lViscosite);

    //On demande un champs pour le temps reduit de la branche courante
    ChampScalaire* lDeltaTempsReduit=0;
    if(lMsg) lMsg =lFichierChamp.reqChamp<ChampScalaire>(lNomDeltaTempsReduit, lDeltaTempsReduit);

    //On demande si l'elasticite depend du temps
    bool lMateriauVieillissant(false);
    if(lMsg) lMsg =lFichierChamp.reqBooleen(lNomVieillissant,lMateriauVieillissant);

    //On demande un champ pour le multiplicateur du tenseur des coefficients d'élasticité de reference
    ChampScalaire* lMultiplicateur_b = 0;
    if(lMsg && lMateriauVieillissant) lMsg =lFichierChamp.reqChamp<ChampScalaire>(lNomb, lMultiplicateur_b);

    //On demande  un champ pour le multiplicateur du tenseur des coefficients d'élasticité de reference au temps precedent
    ChampScalaire* lMultiplicateur_bPrecedent = 0;
    if(lMsg && lMateriauVieillissant){
      lMsg = lFichierChamp.reqChamp<ChampScalaire>(lNombPrecedent,lMultiplicateur_bPrecedent,false);
      if(lMsg && lMultiplicateur_bPrecedent == 0) lMultiplicateur_bPrecedent = lMultiplicateur_b->newChampCopie();
    }

    if(lMsg && lMateriauVieillissant) lMsg = ajouteListeReinterpole(*lMultiplicateur_b,*lMultiplicateur_bPrecedent,
                                                                    lVectChampAReinterpoler);

    //On demande un champ pour l'état du solide dans la branche courante
    ChampScalaire* lIndicateurEtat=0;
    if(lMsg && lMateriauVieillissant) lMsg =lFichierChamp.reqChamp<ChampScalaire>(lNomEtat, lIndicateurEtat);

    //
    // On declare les tenseurs qui sont vraiement utilisés dans le modèle
    //
    ChampExponentielleTempsRelaxation*          lChampTenseurI       = 0;
    ChampCelluleViscoElastiqueMaxwellHPP*       lChampTenseurS       = 0;
    ChampPseudoContraintesO4*                   lChampSigmaO4Courant = 0;

    //On construit le champ Precalculer aux pts d'integration
    ChampPrecalculeDiscPtsIntegration<ChampTensorielO4SymMin>* lChampSigmaO4Precedent = 0;

    if(lMsg){
      //
      //On construit les tenseurs qui sont vraiement utilisés dans le modèle
      //
      lChampTenseurI       = new ChampExponentielleTempsRelaxation;
      lChampTenseurS       = new ChampCelluleViscoElastiqueMaxwellHPP;
      lChampSigmaO4Courant = new ChampPseudoContraintesO4;

      //On construit le champ Precalculer aux pts d'integration
      lChampSigmaO4Precedent = new ChampPrecalculeDiscPtsIntegration<ChampTensorielO4SymMin>;

      lVectChampADeleter.push_back(lChampTenseurI);
      lVectChampADeleter.push_back(lChampTenseurS);
      lVectChampADeleter.push_back(lChampSigmaO4Courant);
      lVectChampADeleter.push_back(lChampSigmaO4Precedent);
    }

    if(lMsg) lMsg = ajouteListeReinterpole(*lChampSigmaO4Courant,*lChampSigmaO4Precedent,lVectChampAReinterpoler);

    // On assigne le champ geometrique
    if(lMsg) {
      lChampTenseurS        ->asgnChampGeometrique(*lChampGeo);
      lChampTenseurI        ->asgnChampGeometrique(*lChampGeo);
      lChampSigmaO4Courant  ->asgnChampGeometrique(*lChampGeo);
      lChampSigmaO4Precedent->asgnChampGeometrique(*lChampGeo);
    }

    // On assigne les parametres aux differents champs tensoriels

    if (lMsg) {
      lChampSigmaO4Precedent                  ->asgnChamp(*lChampSigmaO4Courant);

      lChampTenseurI                          ->asgnChamps(*lElasticite,
                                                           *lViscosite,
                                                           *lDeltaTempsReduit);

      if(lMateriauVieillissant) lChampTenseurS->asgnChampsAvecVieillissement(*lElasticite,
                                                                             *lViscosite,
                                                                             *lDeltaTempsReduit,
                                                                             *lMultiplicateur_b,
                                                                             *lMultiplicateur_bPrecedent,
                                                                             *lIndicateurEtat);

      if(!lMateriauVieillissant) lChampTenseurS->asgnChampsSansVieillissement(*lElasticite,
                                                                              *lViscosite,
                                                                              *lDeltaTempsReduit);

      lChampSigmaO4Courant                     ->asgnChamps(*lChampSigmaO4Precedent,
                                                            *lChampTenseurI,
                                                            *lChampTenseurS,
                                                            *lU,
                                                            *lUPrec);
    }

    if(lMsg) lChampSigmaPrecedent->ajouteSchemaIntegration(lSchemaIntg);

    if(lMsg) {
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].Elasticite                   = lElasticite;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].Viscosite                    = lViscosite;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].MateriauVieillissant         = lMateriauVieillissant;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].DeltaTempsReduit             = lDeltaTempsReduit;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].Multiplicateur_b             = lMultiplicateur_b;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].Multiplicateur_bPrecedent    = lMultiplicateur_bPrecedent;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].IndicateurEtat               = lIndicateurEtat;

      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].TenseurS                     = lChampTenseurS;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].TenseurI                     = lChampTenseurI;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].PseudoContraintesCourantes   = lChampSigmaO4Courant;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].PseudoContraintesPrecedentes = lChampSigmaO4Precedent;
    }

    std::cout<<"On cree le namespace"<<std::endl;

  }

  // Frequence d'impression
  //lReelBidon=1.;
  lPasImpression = 1;
  if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("freqImpression",lPasImpression,false);//lReelBidon,false);
  //  lPasImpression = (Entier) lReelBidon;

  //
  // On permet d'imposer des points dans le temps ou le calcul doit etre fait
  // utile dans le cas où on a des discontinuités en temps dans les données
  // que l'on traite en imposant un calcul sur la discontinuité
  //

  bool lTemps(false);

  if(lMsg) lMsg = lFichierChamp.reqBooleen("TempsImpose",lTemps,false);

  if(lTemps && lMsg){
    // Deux choix:
    // a) les temps sont donnés dans le fichier champs:
    //      scalaire NbTempsImpose N
    //      scalaire TempsImpose_0 t_0
    //      ...
    //      scalaire TempsImpose_N-1 t_N-1
    // b) les temps sont placés dans un fichier "lNom".temps (un temps par ligne avec le nombre de temps impose sur la premiere ligne)

    //    lReelBidon = -1.;
    lNbTempsImpose = 0;
    if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("NbTempsImpose",lNbTempsImpose,false);//lReelBidon,false);
    //    lNbTempsImpose = (Entier) lReelBidon;

    if(lMsg){
      // les temps impose sont dans le fichier champs:
      if(lNbTempsImpose >=0){
	lTempsImposes.asgnLongueur(lNbTempsImpose,0.);
	ChaineCar lTImpose= "TempsImpose_";
	for (Entier i=0; i<lNbTempsImpose && lMsg; ++i){
	  ChaineCar lTempsImposei = ajouteNumeroALaFin(lTImpose,i,1);
	  lMsg = lFichierChamp.reqScalaireConstant(lTempsImposei,lTempsImposes[i]);
	}
      } else {
	// On lit dans un fichier d'extension ".temps"
	ChaineCar lNomTemps = lNomFichierChamp + std::string(".temps");
	std::ifstream lFichier(lNomTemps.c_str());
	if(!lFichier){
	  std::string lMessageUsage1 = "Probleme de temps imposes: pas de points dans fichier champs et";
	  std::string lMessageUsage = "Erreur lors de l'ouverture du fichier: " + lNomTemps;
	  lMsg.ajoute(ERMsg(ERMsg::ERREUR, lMessageUsage));
	}

	if(lMsg) lFichier >> lNbTempsImpose;
	if(lMsg && lNbTempsImpose < 0){
	  std::cout<<"NbTempsImpose est negatif je le remet a 0!"<<std::endl;
	  lNbTempsImpose = 0;
	}

	lTempsImposes.asgnLongueur(lNbTempsImpose,0.);
	if(lTemps && lMsg){
	  for (Entier i=0; i<lNbTempsImpose && lMsg; ++i){
	    lFichier >> lTempsImposes[i];
	  }
	}
      }
    }
  }

  // on demande les coordonnées pour la trace des champs T,M et U:

  Entier lNombreEval(0);
  if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("NbEval",lNombreEval,false);
  if(lNombreEval < 0){
    std::cout<<"Le nombre d'evaluation est negatif je le remet a 0!"<<std::endl;
    lNombreEval = 0;
  }

  if(lNombreEval>0) {
    std::cout<<lNombreEval<<" évaluations des variables d'état aux points: "<<std::endl;
    lCoordonnees.asgnLongueur(lNombreEval);
  }

  if(lNombreEval>0 && lMsg){
    ChaineCar lCoor= "Coord_";
    for (Entier i=0; i<lNombreEval && lMsg; ++i){
      ChaineCar lCoordi = ajouteNumeroALaFin(lCoor,i,1);
      lMsg = lFichierChamp.reqVecteur3DConstant(lCoordi,lCoordonnees[i]);
      std::cout<<lCoordi<<" "<<lCoordonnees[i]<<std::endl;
    }
    std::cout<<std::endl;
  }

  if(lMecanik && lMsg) lProbMec.asgnDonnees(lMail,lGeo,*lChampGeo,lListeEntite);

  //
  // On sépare les données sur types de fichiers: hyd (humidité) et mec (déplacements)
  ChaineCar lNomHyd = lNomMaillage + std::string("_hyd");
  ChaineCar lNomMec = lNomMaillage + std::string("_mec");

  // *********************************************
  // * BLOC: Déclaration des termes de formulation
  // *       et conditions aux limites
  // *********************************************
  // On crée l'objet qui contiendra la liste des conditions limites
  ListeConditionsLimites lListeCondLimHyd, lListeCondLimMec;

  if (lMsg) {
    // On assigne l'objet lFichierChamp aux conditions aux limites, permettant ainsi aux champs
    // définits dans le fichier ".champs" d'être utilisés dans le fichier ".CL"
    lListeCondLimHyd.asgnGestionFichierChamps(lFichierChamp);
    lListeCondLimMec.asgnGestionFichierChamps(lFichierChamp);

    lListeCondLimHyd.asgnChampGeoEtListeEntitesGeo(*lChampGeo, lListeEntite);
    lListeCondLimMec.asgnChampGeoEtListeEntitesGeo(*lChampGeo, lListeEntite);

    // On assigne le schéma d'intégration aux conditions aux limites
    lListeCondLimHyd.asgnSchemaIntg(lSchemaIntg);
    lListeCondLimMec.asgnSchemaIntg(lSchemaIntg);
    std::cout<<"conditions aux limites"<<std::endl;
  }

  // On déclare les termes de formulation que l'on veut utiliser.

  // hydrique
  TFMasseGenScal            lTFMasseM;
  TFDiffusion               lTFDiffusionM;
  TFDiffusion               lTFThermoDiffusion;
  TFDiffusionNewton         lTFDiffusionNewton;

  // thermique
  TFDiffusion               lTFDiffusionT;
  TFMasseGenScal            lTFMasseT;

  //mécanique
  TFCouplageAvecFctsConnues lTFCouplageHydMElas, lTFCouplageHydMVisc;
  TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin        lTFRigiditeComposanteElastiquePure;
  VectDyn<TFCelluleViscoElastiqueMaxwellHPPEpsilonLin> lVectTFRigiditeBrancheViscoElastique(lNombreBrancheViscoElastique);
  TFSourceFVVect3D                                     lTFChargeVolumique;

  // contact
  TFContactV3DUxNVScal  lTFContactGCP;
  TFContactV3DUxNVScal  lTFContactGCP2;

  //On assigne les champs aux termes de formulation.
  if(lMsg) lMsg = lTFDiffusionM. asgnChamps(*lChampGeo, *lM, *lM, *lK_M, false,   lSchemaIntg);
  if(lMsg) lMsg = lTFMasseM.     asgnChamps(*lChampGeo, *lM, *lM, *lC_MRho,lSchemaIntg);
  if(lMsg && lEmploiNewtonEnM) lMsg = lTFDiffusionNewton.asgnChamps(*lChampGeo, *lM, *lM, *lM, *ldK_MdM, lSchemaIntg);

  if(lThermik && lMsg){
    lMsg = lTFThermoDiffusion. asgnChamps(*lChampGeo, *lT, *lM, *lK_MT, false,   lSchemaIntg);
    if(lMsg) lMsg = lTFDiffusionT.  asgnChamps(*lChampGeo, *lT, *lT, *lK_T, false,   lSchemaIntg);
    if(lMsg) lMsg = lTFMasseT.      asgnChamps(*lChampGeo, *lT, *lT, *lC_TRho,lSchemaIntg);
  }

  if(lMecanik && lMsg){
     if(lMsg) lMsg = lTFCouplageHydMElas.asgnChamps(*lChampGeo, *lM, *lMP, *lU, *(lProprietePurementElastique.TenseurS),
				      *lBetaM, lSchemaIntg);
    if (lMsg && lERef_Inf) {
      lMsg = lTFRigiditeComposanteElastiquePure.asgnChamps(*lChampGeo,
                                                           *lU,
                                                           *lU,
                                                           *lUPrec,
                                                           *(lProprietePurementElastique.ContraintesPrecedentes),
                                                           *(lProprietePurementElastique.TenseurS),
                                                           *(lProprietePurementElastique.IndicateurEtat),
                                                           lSchemaIntg);
      if(lMsg) lProbMec.ajouteTF(lTFRigiditeComposanteElastiquePure);
      if(lMsg) lProbMec.ajouteTF(lTFCouplageHydMElas);
       std::cout<<"lProbMec.ajouteTF(lTFRigiditeComposanteElastiquePure)"<<std::endl;
    }

    for(Entier lNumeroDeBranche=0;lNumeroDeBranche<lNombreBrancheViscoElastique && lMsg; ++lNumeroDeBranche){
      ProprieteViscoElastique& lBranche = lVectProprieteBrancheViscoElastique[lNumeroDeBranche];
      TFCelluleViscoElastiqueMaxwellHPPEpsilonLin& lTFRigiditeBranche = lVectTFRigiditeBrancheViscoElastique[lNumeroDeBranche];
      lMsg = lTFRigiditeBranche.asgnChamps(*lChampGeo,
                                           *lU,
                                           *lU,
                                           *lUPrec,
                                           *(lBranche.PseudoContraintesPrecedentes),
                                           *(lBranche.TenseurS),
                                           *(lBranche.TenseurS),
                                           *(lBranche.TenseurI),
                                           lSchemaIntg);

      lTFCouplageHydMVisc.asgnChamps(*lChampGeo, *lM, *lMP, *lU, *(lBranche.TenseurS),*lBetaM, lSchemaIntg);

      if(lMsg) lProbMec.ajouteTF(lTFRigiditeBranche);
      if(lMsg) lProbMec.ajouteTF(lTFCouplageHydMVisc);
      std::cout<<"lProbMec.ajouteTF(lTFRigiditeBranche)"<<std::endl;
    }

    if(lMsg && lF){
      lMsg = lTFChargeVolumique.asgnChamps(*lChampGeo,*lF,*lU,lSchemaIntg);
    }

    if(lMsg && lContakt) lMsg = lTFContactGCP.asgnChamps(*lChampGeoPeau, *lUPeau, *lLambdaCtc,
							 *lNormaleExtPeau,*lGapPeau, lSchemaIntgNewtonCote);
    if(lMsg && lContakt) lMsg = lTFContactGCP2.asgnChamps(*lChampGeoPeau, *lUPeau, *lLambdaCtc,
							  *lNormaleExtPeau,*lGapPeau, lSchemaIntg);
  }

  if(lMsg){
    //On ajoute nos termes de formulations à notre problème
    lProbHyd.ajouteTF    (lTFDiffusionM);
    lProbHyd.ajouteTFInst(lTFMasseM);
    if(lEmploiNewtonEnM) lProbHyd.ajouteTF    (lTFDiffusionNewton);

    if(lThermik){
      lProbHyd.ajouteTFInst(lTFMasseT);
      lProbHyd.ajouteTF    (lTFThermoDiffusion);
      lProbHyd.ajouteTF    (lTFDiffusionT);
    }
    if(lMecanik){
      if(lF){
	lProbMec.ajouteTF(lTFChargeVolumique);
      }

      if(lContakt) lProbPeauGCP.ajouteTF(lTFContactGCP);
      if(lContakt) lProbPeauGCP2.ajouteTF(lTFContactGCP2);
    }
     std::cout<<"On ajoute nos termes de formulations à notre problème"<<std::endl;
  }

  //On lit et initialise nos conditions aux limites
  if(lMsg) lMsg = lProbHyd.lisEtInitialiseConditionsLimites(lNomHyd, lListeCondLimHyd);
  if(lMsg && lMecanik){
    lMsg = lProbMec.lisEtInitialiseConditionsLimites(lNomMec, lListeCondLimMec);
    if(lMsg && lContakt) lProbPeauGCP.initialiseNumerotationPeau(lProbMec.reqNumerotation());
    if(lMsg && lContakt) lProbPeauGCP2.initialiseNumerotationPeau(lProbMec.reqNumerotation());
     std::cout<<"lProbMec.lisEtInitialiseConditionsLimites"<<std::endl;
  }

  // lecture des solveurs
  if(lMsg) {
    lMsg = lFichierChamp.reqSolveur<SolveurInstNlinPETSc>("SolveurInst",lSolveurHyd);
     std::cout<<"END reqSolveur<SolveurInstNlinPETSc"<<std::endl;
  }

   if(lMsg) lMsg = lFichierChamp.reqSolveur<SolveurLinPETSc>("SolveurMec",lSolveurMec);
   //if(lMsg) lMsg = lFichierChamp.reqSolveur<SolveurLinPETSc>("SolveurLinMec",lSolveurMec);
  //if(lMsg) lMsg = lFichierChamp.reqSolveur<SolveurLinPETSc>("SolveurMec",lSolveurMec);
  // On assigne les temps ou l'on force le calcul
  if(lMsg && lTemps){
    lSolveurHyd->asgnTempsImposes(lTempsImposes);
     std::cout<<"END  lSolveurHyd->asgnTempsImposes"<<std::endl;
  }

  SolveurGCPAvecCorpsRigide      *lSolveurGCP = 0;

  Entier lComtpeurIteration = 0;
  DReel lValeurRef(0.);
  DReel lCritere0(1.e-12);

  EntiteGeometrique* lEntite = 0;
  if(lMsg && lContakt) {
     std::cout<<"lMsg && lContakt"<<std::endl;
    if(lMsg && lCorpsParGCP) lEntite =  &(lListeEntite.reqEntiteGeometrique(lNomEntiteControle));

    lSolveurGCP = new SolveurGCPAvecCorpsRigide();
    lSolveurGCP->asgnDonnees(lProbPeauGCP,true);
    lSolveurGCP->asgnDonneesLambda(lProbPeauGCP2);
    if(lMsg && lCorpsParGCP) lSolveurGCP->corrigeCorpsRigide(lProbMec.reqNumerotation(),*lDeplacements,*lU,
							     lDirectionDesCorps,lCritere0,lValeurRef,lEntite);
    lSolveurGCP->asgnParametres (5000,
				 1.e-13,
				 true,
				 false,
				 lComtpeurIteration);
    lSolveurGCP->relaxeTolerance(2,150,1e-10);
    lSolveurMec = lSolveurGCP;
     std::cout<<"lSolveurGCP->asgnDonneesLambda"<<std::endl;
  }

  if(lMecanik && lMsg){
    // A priori la matrice de rigidité depend de M et à besoin d'être reassemblée à chaque fois:
    bool lforceAssemblageMatrice(true);
    lMsg = lFichierChamp.reqBooleen("forceAssemblageMatriceMec",lforceAssemblageMatrice,false);
    lSolveurMec->forceAssemblageMatrice(lforceAssemblageMatrice);
    lSolveurMec->forceAssemblageResidu(lVrai);

    // pour les sauvegardes mecanique

    Entier lEntierBidon(1);
    if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("NbPasDeTemps",lEntierBidon,false);
    if(lMsg) lSolveurMec->asgnLongueurNumeroPourExport(lEntierBidon+lNbTempsImpose);
     std::cout<<"lforceAssemblageMatrice(true)"<<std::endl;
  }

  // Bloc pour l'impression
  ExportGIREF lExportGIREF;
  if(lPasImpression && lMsg){
    lExportGIREF.asgnChampGeometrique(*lChampGeo);
    lExportGIREF.asgnExporteMaillage(false);
    lProbHyd.ajouteExportation(lExportGIREF, lPasImpression, lNomFichierChamp);
    if(lMecanik){
      lMsg = lExportGIREF.ajouteChamp(*lDeplacements, "U");
      lMsg = lExportGIREF.ajouteChamp(*lU, "deltaU");
    }
  }

  ExportVU lExportVU;
  if(lPasImpression && lMsg){
    lExportVU.asgnChampGeometrique(*lChampGeo);
    lProbHyd.ajouteExportation(lExportVU, lPasImpression, lNomHyd);
    if(lMecanik){
      lMsg = lExportVU.ajouteChamp(*lDeplacements, "U");
    }
  }

   // lPPReqTempsHyd       Pour savoir la valeur du temps (pour la passer au probleme mecanique) on demande
  //                      le temps au probleme hydrique.
  // lPPAsgnTempsMec      Le problème mecanique est en post-traitement de l'hydrique alors il nous faut
  //                      lui passer la valeur du temps au moment du calcul.
  // lPPAsgnTempsCtc      Si on a un probleme avec contact on veut s'assurer que le temps est assigne au champs de peau
  //

  PPTraitementGen* lPPReqTempsHyd = 0;
  PPTraitementGen* lPPAsgnTempsMec = 0;
  PPTraitementGen* lPPAsgnTempsCtc = 0;
  PPTraitementGen* lPPAsgnTempsDeplacements= 0;
  DReel lTempsMec;

   if (lMsg) {
    lPPReqTempsHyd  = newPPCallBack(lProbHyd, mem_fun_ref(&ProblemeEF::reqTemps),lTempsMec);
    lPPAsgnTempsMec = newPPCallBack2Args(lProbMec, lTempsMec, mem_fun_ref(&ProblemeEF::asgnTemps));

    if(lContakt) lPPAsgnTempsCtc = newPPCallBack2Args(*lUPeau, lTempsMec, mem_fun_ref(&Champ::asgnTemps));
    lPPAsgnTempsDeplacements = newPPCallBack2Args(*lDeplacements, lTempsMec, mem_fun_ref(&Champ::asgnTemps));
  }

  // Points d'interet dans le domaine:

  // pour M
  PPEvalueChampXYZ<ChampScalContinu> lEvalueM;
  VectDyn<DReel> lMAuxPointsDonnes(lNombreEval,0.);

  // pour T
  PPEvalueChampXYZ<ChampScalContinu> lEvalueT;
  VectDyn<DReel> lTAuxPointsDonnes(lNombreEval,0.);

  //pour U
  PPEvalueChampXYZ<ChampVect3DContinu> lEvalueU;
  VectDyn<Vecteur3D> lUAuxPointsDonnes(lNombreEval);


  // pour l'ecriture des evaluations ponctuelles:
  PPAfficheValeur<VectDyn<DReel> > lAfficheM;
  ofstream lFichierM;

  PPAfficheValeur<VectDyn<DReel> > lAfficheT;
  ofstream lFichierT;

  PPAfficheValeur<VectDyn<Vecteur3D> > lAfficheUz;
  ofstream lFichierUz;

  // On a donné dans le fichier champs des points ou on veut une évaluation:
  if(lMsg && lNombreEval){
    lFichierM.open((lNomMaillage + ".M").c_str());
    lEvalueM.asgnParametres(*lChampGeo,*lM,lCoordonnees,lMAuxPointsDonnes);
    lAfficheM.afficheLeTemps(lProbHyd,lVrai);
    lAfficheM.asgnParametres("M",lMAuxPointsDonnes,lFichierM);
    lFichierM << std::setprecision(8);
    lFichierM << std::scientific;
  }

  if(lThermik && lNombreEval && lMsg){
    lFichierT.open((lNomMaillage + ".T").c_str());
    lEvalueT.asgnParametres(*lChampGeo,*lT,lCoordonnees,lTAuxPointsDonnes);
    lAfficheT.afficheLeTemps(lProbHyd,lVrai);
    lAfficheT.asgnParametres("T",lTAuxPointsDonnes,lFichierT);
    lFichierT << std::setprecision(8);
    lFichierT << std::scientific;
  }

  if(lMecanik && lNombreEval && lMsg){
    lFichierUz.open((lNomMaillage + ".U").c_str());
    lEvalueU.asgnParametres(*lChampGeo,*lU,lCoordonnees,lUAuxPointsDonnes);
    lAfficheUz.afficheLeTemps(lProbHyd,lVrai);
    lAfficheUz.asgnParametres("U",lUAuxPointsDonnes,lFichierUz);
    lFichierUz << std::setprecision(8);
    lFichierUz << std::scientific;
  }

  // On sauvegarde la valeur du pas de temps à chaque pas
  // dans un fichier prefixe.DT utilise pour l'adaptation en pas variable.
  ofstream lFichierDT((lNomMaillage + ".DT").c_str());
  lFichierDT << std::setprecision(12);
  lFichierDT << std::scientific;

  PPAfficheValeur <DReel> lAfficheTemps;
  DReel *lDeltaT = 0;
  if (lMsg) {
    lSolveurHyd->reqPtrDelta(lDeltaT);
    lAfficheTemps.afficheLeTemps(lProbHyd,lVrai);
    lAfficheTemps.asgnParametres("",*lDeltaT, lFichierDT); // affiche le pas de temps dans le fichier
    lProbHyd.ajoutePostTraitement(lAfficheTemps);   // On ajoute l'affichage au post traitement du probleme hydrique
  }

  // Le probleme mecanique est resolu en post-traitement du probleme hydrique:
  PPResolutionProbleme lPPResolution;
  if(lMsg && lMecanik) lPPResolution.asgnParametres(lProbMec,*lSolveurMec);

  // en cas d'hysterese on doit conserver la valeur de M au pas precedent
  PPReinterpole<ChampScalContinu> lStockeMAuPasPrecedent;
  if (lMsg && lMecanik && lBetaM && lHysterese) lStockeMAuPasPrecedent.asgnParametres(*lM,*lMP);

   // ajout de l'incrément U aux Deplacements
  PPCalculChampaXpbY<ChampVect3DContinu> lMAJU;
  DReel lUnReel(1.);
  if(lMecanik && lMsg) lMAJU.asgnParametres(lUnReel,lUnReel,*lDeplacements,*lU,*lDeplacements);


  // On ajoute en Pretraitement la reinterpolation de U dans UPrecedent
  PPReinterpole<ChampVectoriel3D> lReinterpoleU;
  if(lMsg) {
    lReinterpoleU.asgnParametres(*lU,*lUPrec);
  }


  // ********************
  // * BLOC: Résolution *
  // ********************

  //
  // On a déja défini le solveur quasi statique maintenant on optimise, si le
  // solveur était défini dans le fichier champs on a rien a faire de tout cela...
  //

  if(!lDefinitionDansFichierChamps && lMsg){
    // Le solveur lineaire de base
    if(lMsg) lMsg = lFichierChamp.reqSolveur<SolveurLinPETSc>("SolveurMec",lSolveurMec,false);

    if(lSolveurMec && lMsg) lSolveurQuasiStatique->asgnSolveurLin(*lSolveurMec);

    if(lMsg){
      // A priori la matrice de rigidité depend de M et à besoin d'être reassemblée à chaque fois:
      lSolveurMec->forceAssemblageMatrice(lReassemble);
      lSolveurMec->forceAssemblageResidu(lVrai);
    }

    //On demande le temps initial
    DReel lTempsInitial = 0.;
    if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("TempsInitial",lTempsInitial,false);

    //On demande la longueur du pas de temps
    DReel lLongueurPasDeTemps = 0.;
    if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("LongueurPasDeTemps",lLongueurPasDeTemps);

    //On demande le nombre de pas de temps
    Entier lNombreDePasDeTemps = 0;
    if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("NombreDePasDeTemps",lNombreDePasDeTemps);

    // On demande si on veut un suivi à chaque pas de temps
    bool lVerboseInstationnaire(false);
    if(lMsg) lMsg = lFichierChamp.reqBooleen("VerboseInstationnaire",lVerboseInstationnaire,false);

    if(lMsg) lSolveurQuasiStatique->asgnParametresInst(lLongueurPasDeTemps,lNombreDePasDeTemps,lTempsInitial,lVerboseInstationnaire);

    // On demande si on veut faire l'adaptation de la longueur de pas de temps
    bool lUtilisePasAdapte(false);
    if(lMsg) lMsg = lFichierChamp.reqBooleen("UtilisePasAdapte",lUtilisePasAdapte,false);

    if(lUtilisePasAdapte && lMsg) {

      //On demande la longueur du pas de temps
      DReel lLongueurMaximalePasDeTemps = 0.;
      if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("LongueurMaximalePasDeTemps",lLongueurMaximalePasDeTemps);

      //On demande le temps initial
      DReel lTempsFinal = 0.;
      if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("TempsFinal",lTempsFinal);

      //On demande le temps initial
      DReel lToleranceAdaptation = 0.;
      if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("ToleranceAdaptation",lToleranceAdaptation);

      //
      // On utilise un pas de temps adapté
      //
      if(lMsg) lSolveurQuasiStatique->utilisePasDeTempsVariable(lLongueurPasDeTemps,lLongueurMaximalePasDeTemps,
                                                                lTempsFinal,lToleranceAdaptation);
        std::cout<<"lDefinitionDansFichierChamps && lMsg)"<<std::endl;
    }

    //On demande la tolerance pour le point fixe
    DReel lTolerancePointFixe = 1.e-6;
    if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("TolerancePointFixe",lTolerancePointFixe,false);

    //On demande le nombre maximal d'iterations de point fixe
    Entier lNombreMaximalIterationsPointFixe = 250;
    if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("NombreMaximalIterationsPointFixe",lNombreMaximalIterationsPointFixe,false);

    // On demande si on veut un suivi du point fixe
    bool lVerbosePointFixe(false);
    if(lMsg) lMsg = lFichierChamp.reqBooleen("VerbosePointFixe",lVerbosePointFixe,false);

    if(lMsg) lSolveurQuasiStatique->asgnParametresNlin(lTolerancePointFixe,lNombreMaximalIterationsPointFixe,lVerbosePointFixe);

    // On demande les proprietes de la matrice
    bool lMatriceDependDuTemps(true);
    bool lMatriceDependDeLaSolution(true);
    bool lMatriceDependDesCL(true);
    if(lMsg) lMsg = lFichierChamp.reqBooleen("MatriceDependTemps",lMatriceDependDuTemps,false);
    if(lMsg) lMsg = lFichierChamp.reqBooleen("MatriceDependSolution",lMatriceDependDeLaSolution,false);
    if(lMsg) lMsg = lFichierChamp.reqBooleen("MatriceDependCL",lMatriceDependDesCL,false);

    if(lMsg) lSolveurQuasiStatique->asgnParametresMatrice(lMatriceDependDuTemps,lMatriceDependDeLaSolution,lMatriceDependDesCL);
  }

  //
  // on ajoute les pre/post traitement ATTENTION la chronologie est importante!
  //

  if(lMecanik && lMsg){
    // on calcule U apres M:
    lProbHyd.ajoutePostTraitement(*lPPReqTempsHyd);
    lProbHyd.ajoutePostTraitement(*lPPAsgnTempsMec);
    if(lContakt) lProbHyd.ajoutePostTraitement(*lPPAsgnTempsCtc);
    lProbHyd.ajoutePostTraitement(lPPResolution,lPasDeCalcul);

     lProbMec.ajoutePostTraitement(lMAJU); // ajoute l'increment aux deplacements
     lProbMec.ajoutePostTraitement(*lPPAsgnTempsDeplacements);

    if(lMsg && lNombreEval>0){
      lProbHyd.ajoutePostTraitement(lEvalueU);
      lProbHyd.ajoutePostTraitement(lAfficheUz);
    }

  }

  if(lMsg && lNombreEval>0){
    lProbHyd.ajoutePostTraitement(lEvalueM);
    lProbHyd.ajoutePostTraitement(lAfficheM);
    if(lThermik){
      lProbHyd.ajoutePostTraitement(lEvalueT);
      lProbHyd.ajoutePostTraitement(lAfficheT);
    }
  }

  // mise a jour de MP  pour le "prochain" pas dans le cas ou on a hysterese
  if (lMsg && lMecanik && lBetaM && lHysterese) lProbHyd.ajoutePostTraitement(lStockeMAuPasPrecedent);
  if (lMsg&& lMecanik) lProbMec.ajoutePreTraitement(lReinterpoleU);

  //
  // Résolution du problème en T, M et en U
  //

  if (lMsg) {
    // On résoud
    if(!lThermik && !lMecanik) std::cout<<"Resolution du probleme d'humidité"<<std::endl;
    if(lThermik && !lMecanik) std::cout<<"Resolution du probleme température/humidité"<<std::endl;
    if(!lThermik && lMecanik) std::cout<<"Resolution du probleme hygro-mécanique";
    if(lThermik && lMecanik) std::cout<<"Resolution du probleme thermo-hygro-mécanique";
    if(lContakt) std::cout<<" avec contact"<<std::endl;
    std::cout << std::endl << "Ecriture des solutions aux différents pas de temps dans les fichiers de préfixe ";
    std::cout << lNomFichierChamp<<std::endl;
    if(lNombreEval>0){
      std::cout << std::endl;
      std::cout << "Ecriture des traces des solutions aux différents pas de temps dans les fichiers de préfixe ";
      std::cout << lNomMaillage<<std::endl;
    }
    lMsg = lSolveurHyd->resoudre(lProbHyd);


    // On assigne les temps ou l'on force le calcul
    if(lMsg && lTemps) lSolveurQuasiStatique->asgnTempsImposes(lTempsImposes);

     //
    // On s'assure d'avoir tout les champs initialisé en faisant une reinterpolation AVANT l'appel au solveur puis
    // on ajoute les reinterpolations en post traitement pour s'assurer d'avoir les bonnes valeurs de temps
    //

    std::vector<PPTraitementGen* >::iterator lIterReinterpole = lVectChampAReinterpoler.begin();
    const std::vector<PPTraitementGen* >::iterator lIterReinterpoleFin = lVectChampAReinterpoler.end();
    while(lIterReinterpole != lIterReinterpoleFin && lMsg){
      (*lIterReinterpole)->effectueCalcul();
      lProbMec.ajoutePostTraitement(**lIterReinterpole);
      //lProbHyd.ajoutePostTraitement(**lIterReinterpole);
      ++lIterReinterpole;
    }

    //  lMsg = lSolveurHyd->resoudre(lProbHyd);

    //
    // On fait le menage, on detruit tout les ptr crees manuellement
    //
    lIterReinterpole = lVectChampAReinterpoler.begin();
    while(lIterReinterpole != lIterReinterpoleFin && lMsg){
      delete *lIterReinterpole;
    }

    std::vector<Champ*>::iterator       lIterDelete = lVectChampADeleter.begin();
    const std::vector<Champ*>::iterator lIterDeleteFin = lVectChampADeleter.end();

    while(lIterDelete != lIterDeleteFin && lMsg){
      delete *lIterDelete;
    }

  }

  if(lNombreEval>0) lFichierM.close();

  if(lMecanik && lNombreEval>0){
    lFichierUz.close();
  }

  if(lThermik && lNombreEval>0){
    lFichierT.close();
  }

  //Traitement d'erreur: si une erreur s'est produite, on l'affiche.
  if(!lMsg) {
    std::cerr << traducteurERMsg(lMsg) << endl;
    return 1;
  }

  return 0;
}

ERMsg lisString(const char*& pCharPtrAModifier,std::string& pStringRetour)
{
  ERMsg lMsg;

  while(espaceBlanc(*pCharPtrAModifier)) {
    ++pCharPtrAModifier;
  }

  if(*pCharPtrAModifier == '\0') {
    lMsg = ERMsg(ERMsg::ERREUR,"ERR_LECTURE_STRING");
    lMsg.ajoute("FIN_CHAINE");
  }
  else if(*pCharPtrAModifier == '\"') {
    ++pCharPtrAModifier;
    while(*pCharPtrAModifier != '\0' && *pCharPtrAModifier != '\"') {
      pStringRetour += *pCharPtrAModifier;
      ++pCharPtrAModifier;
    }
    if(*pCharPtrAModifier == '\0') {
      lMsg = ERMsg(ERMsg::ERREUR,"ERR_DEFINITION_CHAINE");
    }
    else {
      pCharPtrAModifier++;
    }
  }
  else {
    pStringRetour = "";
    while(*pCharPtrAModifier != '\0' && !espaceBlanc(*pCharPtrAModifier)) {
      pStringRetour += *pCharPtrAModifier;
      ++pCharPtrAModifier;
    }
  }
  return lMsg;
}