//************************************************************************
// --- 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:
//
//  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(E:(epsilon(U) - BetaM*(M-MInitiale) - BetaT)) = 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.
//
//  U est obtenu de manière incrémental i.e. par correction à chaque pas de temps. Il faut donc
//  dans les fichiers champs imposé les CL, CI et efforts incrémentaux!!!
//
//  USAGE:<pre> ThermoHygroMecanik.[dev,opt]  préfixe_des_fichiers_d_entrée [entite_en_contact]
//
//  Exemple:
//
//
//********************************************************************

#include "MEFPPUtil.h"
#include "ChaineCar.h"
#include "ChampGeometrique.h"
#include "ChampGeoQuad.h"
#include "ChampGeoSAFE.h"
#include "ChampScalaire.h"
#include "ChampScalContinu.h"
#include "ChampScalAConstant.h"
#include "ChampTensO2SymContinu.h"
#include "ChampTensO2SymConstElem.h"
#include "ChampTensO2SymLinElem.h"
#include "ChampTensO2SymLin.h"
#include "ChampTensorielO4Sym.h"
#include "ChampVect3DContinu.h"
#include "ChampContraintesHPP.h"
#include "ExportGIREF.h"
#include "ExportVU.h"
#include "ExportVTK.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 "PPContraintesEL3D.h"
#include "PPEvalueChampXYZ.h"
#include "PPEvalueChampSurCourbe.h"
#include "PPResolutionProbleme.h"
#include "PPReinterpole.h"
#include "PPCallBack2Args.h"
#include "PPProjectionL2.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 "SolveurGCP.h"
#include "TFDiffusion.h"
#include "TFDiffusionNewton.h"
#include "TFElasLinMixteGen.h"
#include "TFMasseGenScal.h"
#include "TFMasseMixteGeneralise.h"
#include "TFContraintesInitiales.h"
#include "TFCouplageAvecFctsConnues.h"
#include "TFContactV3DUxNVScal.h"
#include "TFHPPMixteAnisotropeUU.h"
#include "TFHPPMixteAnisotropeUP.h"
#include "TFHPPMixteAnisotropePU.h"
#include "TFSigmaUEpsilonV3D.h"
#include "TFSourceFVVect3D.h"
#include "VerifieElementsNormalises.h"
#include "traducteurERMsg.h"
#include <fstream>
#include <iomanip>

using namespace std;

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;

  // 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 ) {
    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 relative
  ChampScalContinu *lVraiM = 0;                  // teneur en humidite
  ChampScalContinu *lVraiMP = 0;                 // teneur en 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

  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;                // l'increment des déplacements
  ChampScalContinu *lP = 0;                  // l'increment de la pression
  ChampVect3DContinu *lDeplacements = 0;     // les déplacements
  ChampScalContinu *lPression = 0;           // la pression
  ChampTensorielO4Sym *lE     = 0;           // tenseur O4 des coeff. d'élasticité

  ChampTensO2SymContinu *lBetaM = 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
  ChampScalContinu        *lPPeau = 0;           // Une version "peau" de P
  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   lFormulationIncrementale(true); // la formulation mecanique est incrementale
  bool   lHPPMixte(false);               // la composante mecanique esr résolue avec une formulation mixte ou non
  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 = 0; // le solveur instationnaire pour M (et T possiblement)
  SolveurLinPETSc      *lSolveurMec = 0; // le solveur stationnaire pour U
  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<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), lFaux(false);

  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) {
    lMsg = lFichierChamp.asgnDonneesBase(&lMail);
    if (lMsg) 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());

  // *****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 && lMecanik) lMsg = lFichierChamp.reqBooleen("HPPMixte",lHPPMixte,false);

  if(lMsg && lContakt && !lMecanik) lMsg.ajoute(ERMsg(ERMsg::ERREUR,"Pas de contact sans mecanique!"));

  // On demande le champ de transformation géométrique au fichier de champs
  if (lMsg) {
    lMsg = lFichierChamp.reqChamp<ChampGeometrique>("ChampGeo",lChampGeo);
  }
  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>("RH",lM);
  //On demande  un champ qui se nomme "M" pour l'humidité
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("M",lVraiM);

  //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);

  //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);

  //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);

  //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 && ldK_MdM) lEmploiNewtonEnM = true;

  //
  // 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
  //  DReel lReelBidon(1.);
  lPasDeCalcul = 1;
  if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("freqCalcul",lPasDeCalcul,false); //lReelBidon,false);
  //  lPasDeCalcul = (Entier) lReelBidon;

  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  si la formluation en deplacements est incrementale
    if(lMsg) lMsg = lFichierChamp.reqBooleen("Formulation_incrementale",lFormulationIncrementale);

    // On demande  un champ qui se nomme "U" pour le deplacement engendré par l'humidité.
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampVect3DContinu>("U",lU);
    if(lMsg){
      if(lFormulationIncrementale){
        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();
      } else {
        lDeplacements = lU;
      }
    }

    // On demande un champs qui se nomme "E" pour le tenseur des coefficients d'élasticité
    if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensorielO4Sym>("Eijkl",lE);

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

    // On est en mixte il faut la pression
    if(lMsg && lHPPMixte){
      lMsg = lFichierChamp.reqChamp("P",lP);
      if(lMsg) lMsg = lFichierChamp.reqChamp("Pression",lPression,false);
      if(!lPression && lMsg) lPression = lP->newChampCopie();
    }

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

    // 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");
    }

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

    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 && lHPPMixte) lMsg = lFichierChampPeau.reqChamp("P",lPPeau);
      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);
    }
  }

  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;
  }


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

  // Bloc pour l'impression
  Exportation lExporteur;
  Exportation lExporteurPropMat;
  PiloteExportationGIREF  lPiloteGIREF;
  lPiloteGIREF.asgnExporteUneFoisMaillage(false);
  PiloteExportationVTK    lPiloteVTK;
  PiloteExportationVU     lPiloteVU;

  if(lPasImpression > 0) {
    //
    // On demande le format d'exportation
    ChaineCar lTypeExportation("GIREF");
    if (lMsg) {

       lMsg = lFichierChamp.reqChaineCar("Format_Exportation",lTypeExportation,false);
       lTypeExportation = mettreEnMajuscule(lTypeExportation);

       bool lAtLeastOneExportation = false;
       if(lTypeExportation.find("GIREF") != ChaineCar::npos){
         lExporteur.ajoutePiloteExportation(lPiloteGIREF);
         lExporteurPropMat.ajoutePiloteExportation(lPiloteGIREF);
         lAtLeastOneExportation = true;
       }

       if(lTypeExportation.find("VU") != ChaineCar::npos){
         lExporteur.ajoutePiloteExportation(lPiloteVU);
         lExporteurPropMat.ajoutePiloteExportation(lPiloteVU);
         lAtLeastOneExportation = true;
       }

       if(lTypeExportation.find("VTK") != ChaineCar::npos){
         lExporteur.ajoutePiloteExportation(lPiloteVTK);
         lExporteurPropMat.ajoutePiloteExportation(lPiloteVTK);
         lAtLeastOneExportation = true;
       }

       if (not lAtLeastOneExportation) {
           lMsg = ERMsg(ERMsg::ERREUR, "No export format chosen. You have the choices: GIREF, VU, VTK. Multiple choices are possible.");
       }

       if (lMsg) {

         lExporteurPropMat.asgnChampGeometrique(*lChampGeo);
         lExporteur.asgnChampGeometrique(*lChampGeo);

         // On vérifie si le degré d'interpolation a été choisi, sinon on en choisi un.
         if (dynamic_cast<const ChampGeoQuad*>(lChampGeo)) {
            lExporteurPropMat.asgnDegreInterpolation(Exportation::eMAILLAGE_QUADRATIQUE);
            lExporteur.asgnDegreInterpolation(Exportation::eMAILLAGE_QUADRATIQUE);
         }
         else if (dynamic_cast<const ChampGeoSAFE*>(lChampGeo)) {
            lExporteurPropMat.asgnDegreInterpolation(Exportation::eMAILLAGE_SAFE);
            lExporteur.asgnDegreInterpolation(Exportation::eMAILLAGE_SAFE);
         }
         else {
            lExporteurPropMat.asgnDegreInterpolation(Exportation::eMAILLAGE_LINEAIRE);
            lExporteur.asgnDegreInterpolation(Exportation::eMAILLAGE_LINEAIRE);
         }
       }
    }
    if(lMsg) std::cout<<"Vous avez choisi une exportation au format "<<lTypeExportation<<std::endl;
  }

  //
  // 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);
  }

  // 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 lTFCouplageHydM, lTFCouplageHydT;
  TFSourceFVVect3D          lTFChargeVolumique;

  //HPP deplacemment
  TFSigmaUEpsilonV3D        lTFRigidite;
  //HPP Mixte
  TFHPPMixteAnisotropeUU lTFHPPMixteUU;
  TFHPPMixteAnisotropeUP lTFHPPMixteUP;
  TFHPPMixteAnisotropePU lTFHPPMixtePU;
  TFMasseMixteGeneralise lTFHPPMixtePP;

  // 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(!lHPPMixte){
      lMsg = lTFRigidite.asgnChamps(*lChampGeo, *lU, *lU, *lE, lSchemaIntg);
    } else {
      lTFHPPMixteUP.optimisationConstructionDuBlocUP(lVrai,lVrai);
      lTFHPPMixtePU.optimisationConstructionDuBlocPU(lFaux,lFaux);

      lTFHPPMixteUU.asgnChamps(*lChampGeo, *lU, *lU, *lE, lSchemaIntg);

      lTFHPPMixteUP.asgnChamps(*lChampGeo, *lP,*lU, *lE, lSchemaIntg);
      lTFHPPMixtePU.asgnChamps(*lChampGeo, *lU,*lP, *lE, lSchemaIntg);

      lTFHPPMixtePP.asgnChamps(*lChampGeo, *lP, *lE, lSchemaIntg);
      lTFHPPMixtePP.asgnAssembleResidu(lVrai);

    }
    if(lMsg) lMsg = lTFCouplageHydM.asgnChamps(*lChampGeo, *lVraiM, *lVraiMP, *lU, *lE,
                      *lBetaM, lSchemaIntg);
    if(lMsg && lBetaT){
      lMsg = lTFCouplageHydT.asgnChamps(*lChampGeo, *lUn, *lU, *lE,
                    *lBetaT, lSchemaIntg);
    }

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

    ChampScalAConstant* lPtrNULL = 0;
    if(lMsg && lContakt) lMsg = lTFContactGCP.asgnChamps(*lChampGeoPeau, *lUPeau, *lLambdaCtc,
                             *lNormaleExtPeau,*lGapPeau,lPtrNULL,lSchemaIntgNewtonCote);
    if(lMsg && lContakt) lMsg = lTFContactGCP2.asgnChamps(*lChampGeoPeau, *lUPeau, *lLambdaCtc,
                              *lNormaleExtPeau,*lGapPeau,lPtrNULL,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(!lHPPMixte){
        lProbMec.ajouteTF    (lTFRigidite);
      } else {
        lProbMec.ajouteTF(lTFHPPMixteUU);
        lProbMec.ajouteTF(lTFHPPMixtePU);
        lProbMec.ajouteTF(lTFHPPMixteUP);
        lProbMec.ajouteTF(lTFHPPMixtePP);
      }
      lProbMec.ajouteTF    (lTFCouplageHydM);
      if(lF){
    lProbMec.ajouteTF(lTFChargeVolumique);
      }
      if(lBetaT){
    lProbMec.ajouteTF    (lTFCouplageHydT);
      }
      if(lContakt) lProbPeauGCP.ajouteTF(lTFContactGCP);
      if(lContakt) lProbPeauGCP2.ajouteTF(lTFContactGCP2);
    }
  }

  //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());
  }

  // lecture des solveurs
  if(lMsg) lMsg = lFichierChamp.reqSolveur<SolveurInstNlinPETSc>("SolveurInst",lSolveurHyd);
  if(lMsg) lMsg = lFichierChamp.reqSolveur<SolveurLinPETSc>("SolveurLinmecanique",lSolveurMec);
  // On assigne les temps ou l'on force le calcul
  if(lMsg && lTemps) lMsg = lSolveurHyd->asgnTempsImposes(lTempsImposes);

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

  EntiteGeometrique* lEntite = 0;

  if(lMsg && lContakt) {
    if(lMsg && lCorpsParGCP) {
      lEntite =  &(lListeEntite.reqEntiteGeometrique(lNomEntiteControle));
      SolveurGCPAvecCorpsRigide* lSolveurGCP = new SolveurGCPAvecCorpsRigide();
      lSolveurGCP->corrigeCorpsRigide(lProbMec.reqNumerotation(),*lDeplacements,*lU,
                      lDirectionDesCorps,lCritere0,lValeurRef,lEntite);

      lSolveurGCP->asgnDonnees(lProbPeauGCP,true);
      lSolveurGCP->asgnDonneesLambda(lProbPeauGCP2);
      lSolveurGCP->relaxeTolerance(2,150,1e-10);

	  lSolveurGCP->asgnParametres (5000,
				   1.e-13,
				   true,
				   false,
				   lComtpeurIteration);
	  lSolveurMec = lSolveurGCP;
	} else {
	  SolveurGCP* lSolveurGCP = new SolveurGCP();

      lSolveurGCP->asgnDonnees(lProbPeauGCP);
      lSolveurGCP->asgnDonneesLambda(lProbPeauGCP2);
      lSolveurGCP->relaxeTolerance(2,150,1e-10);

	  lSolveurGCP->asgnParametres (5000,
				   1.e-13,
				   1.e-13,
				   false,
				   false,
				   lComtpeurIteration);
	  lSolveurMec = lSolveurGCP;
	}

  }

  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);
  }

  // Bloc pour l'impression des propriétés des matériaux
  if(lMsg && 1==0){
    ChampScalaire* lVraia = 0;
    lMsg = lFichierChamp.reqChamp("a",lVraia);
    ChampScalLin la(*lChampGeo);
    if(lMsg){
      PPReinterpole lPPReinterpole0;
      if(lMsg) lMsg = lPPReinterpole0.asgnParametres(*lVraia,la);
      if(lMsg) lMsg = lPPReinterpole0.effectueCalcul();
    }

    if(lMsg) lMsg = lExporteurPropMat.ajouteChamp(la,"a");

    ChampScalLin lCMRho(*lChampGeo);
    PPReinterpole lPPReinterpole;
    if(lMsg) lMsg = lPPReinterpole.asgnParametres(*lC_MRho,lCMRho);
    if(lMsg) lMsg = lPPReinterpole.effectueCalcul();

    if(lMsg) lMsg = lExporteurPropMat.ajouteChamp(lCMRho,"CMRho");

    ChampTensO2SymLin lKM(*lChampGeo);
    PPReinterpole lPPReinterpole1;
    if(lMsg) lMsg = lPPReinterpole1.asgnParametres(*lK_M,lKM);
    if(lMsg) lMsg = lPPReinterpole1.effectueCalcul();

    if(lMsg) lMsg = lExporteurPropMat.ajouteChamp(lKM,"K_M");

    ChampTensO2SymLin lBM(*lChampGeo);
    PPReinterpole lPPReinterpole2;
    if(lMsg && lBetaM) lMsg = lPPReinterpole2.asgnParametres(*lBetaM,lBM);
    if(lMsg) lMsg = lPPReinterpole2.effectueCalcul();

    if(lMsg && lBetaM) lMsg = lExporteurPropMat.ajouteChamp(lBM, "BetaM");

    ChampScalaire* lELAna = 0;
    if(lMsg) lMsg = lFichierChamp.reqChamp("EL",lELAna);
    ChampScalLin lEL(*lChampGeo);
    if(lMsg) lMsg = lPPReinterpole.asgnParametres(*lELAna,lEL);
    if(lMsg) lMsg = lPPReinterpole.effectueCalcul();

    if(lMsg) lMsg = lExporteurPropMat.ajouteChamp(lEL,"EL");

    ChampScalaire* lETAna = 0;
    if(lMsg) lMsg = lFichierChamp.reqChamp("ET",lETAna);
    ChampScalLin lET(*lChampGeo);
    if(lMsg) lMsg = lPPReinterpole.asgnParametres(*lETAna,lET);
    if(lMsg) lMsg = lPPReinterpole.effectueCalcul();

    if(lMsg) lMsg = lExporteurPropMat.ajouteChamp(lET,"ET");

    ChampScalaire* lERAna = 0;
    if(lMsg) lMsg = lFichierChamp.reqChamp("ER",lERAna);
    ChampScalLin lER(*lChampGeo);
    if(lMsg) lMsg = lPPReinterpole.asgnParametres(*lERAna,lER);
    if(lMsg) lMsg = lPPReinterpole.effectueCalcul();

    if(lMsg) lMsg = lExporteurPropMat.ajouteChamp(lER,"ER");


    ChampScalLinElem lE0000(*lChampGeo);
    if(lMsg) lMsg = lPPReinterpole.asgnParametres(lE->reqComposante(0),lE0000);
    if(lMsg) lMsg = lPPReinterpole.effectueCalcul();

    if(lMsg) lMsg = lExporteurPropMat.ajouteChamp(lE0000,"E0000");

    ChampScalLinElem lE1111(*lChampGeo);
   if(lMsg) lMsg =  lPPReinterpole.asgnParametres(lE->reqComposante(5),lE1111);
    if(lMsg) lMsg = lPPReinterpole.effectueCalcul();

    if(lMsg) lMsg = lExporteurPropMat.ajouteChamp(lE1111,"E1111");

    ChampScalLinElem lE2222(*lChampGeo);
    if(lMsg) lMsg = lPPReinterpole.asgnParametres(lE->reqComposante(20),lE2222);
    if(lMsg) lMsg = lPPReinterpole.effectueCalcul();

    if(lMsg) lMsg = lExporteurPropMat.ajouteChamp(lE2222,"E2222");


    if(lMsg) lMsg = lExporteurPropMat.exporteInfo("valeurs_materiaux");
  }

  //return 0;
  ChampScalaire* lDM=0;
  if(lMsg) lMsg = lFichierChamp.reqChamp("DeltaM",lDM);

  std::cout << " no 9 "<<traducteurERMsg(lMsg) << endl;
  // Bloc pour l'impression
  if(lPasImpression && lMsg){
    lProbHyd.ajouteExportation(lExporteur, lPasImpression, lNomFichierChamp+"_");
    lMsg = lExporteur.ajouteChamp(*lDM,"deltaM");
    if(lMsg) lMsg = lExporteur.ajouteChamp(*lVraiM,"vraiM");

    if(lMecanik){
      lMsg = lExporteur.ajouteChamp(*lDeplacements, "U");
      if(lMsg) lMsg = lExporteur.ajouteChamp(*lU, "deltaU");
    }
  }

  std::cout << " no 10 "<<traducteurERMsg(lMsg) << endl;
  // Debut des pre/post traitements

  // 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 && lMecanik) {
    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:

  std::cout << " no 11 "<<traducteurERMsg(lMsg) << endl;
  // 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,*lDeplacements,lCoordonnees,lUAuxPointsDonnes);
    lAfficheUz.afficheLeTemps(lProbHyd,lVrai);
    lAfficheUz.asgnParametres("U",lUAuxPointsDonnes,lFichierUz);
    lFichierUz << std::setprecision(8);
    lFichierUz << std::scientific;
  }

  // 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 lStockeMAuPasPrecedent;
  if(lMecanik && lFormulationIncrementale && lMsg) lStockeMAuPasPrecedent.asgnParametres(*lVraiM,*lVraiMP);

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


  //
  // 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);

    if(lFormulationIncrementale) lProbMec.ajoutePostTraitement(lMAJU); // ajoute l'increment aux deplacements
    if(lHPPMixte && lFormulationIncrementale) lProbMec.ajoutePostTraitement(lMAJP); // ajoute l'increment aux deplacements
    lProbMec.ajoutePostTraitement(*lPPAsgnTempsDeplacements);

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

  }

  std::cout << " no 12 "<<traducteurERMsg(lMsg) << endl;
  if(lMsg && lNombreEval>0){
    lEvalueM.effectueCalcul();
    lAfficheM.effectueCalcul();
    lProbHyd.ajoutePostTraitement(lEvalueM);
    lProbHyd.ajoutePostTraitement(lAfficheM);
    if(lThermik){
    lEvalueT.effectueCalcul();
    lAfficheT.effectueCalcul();
      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 && lFormulationIncrementale && lBetaM) lProbMec.ajoutePostTraitement(lStockeMAuPasPrecedent,lPasDeCalcul);

  // Traitement des contraintes
  ChampTensO2SymConstElem lSigmaConst;
  if(lMsg){
    lSigmaConst.asgnChampGeometrique(*lChampGeo);
    lSigmaConst.asgnNom("SigmaConst");
  }

  ChampContraintesHPP lSigmaAna("SigmaAna");
  PPProjectionL2 lPPProjectionL2SigmaAna;

  if (lMsg && lMecanik){
    lSigmaAna.asgnChampGeometrique(*lChampGeo);
    lSigmaAna.asgnChamps(*lU,*lE);
    lSigmaAna.asgnCouplage(*lVraiM,*lVraiMP,*lBetaM);

    if (1==0){
      std::cout<<"DEBUGJD A Revoir si formulation incrementale "<<std::endl;
      if(lFormulationIncrementale) lSigmaAna.asgnContraintesInitiales(lSigmaConst);
      lPPProjectionL2SigmaAna.asgnChampGeometrique(*lChampGeo);
      lPPProjectionL2SigmaAna.asgnParametres(0, lSchemaIntg,lSigmaAna,lSigmaConst);
      lProbMec.ajoutePostTraitement(lPPProjectionL2SigmaAna); // Fait la mise a jour du tenseur de contraintes lissé (projection L2)
      lMsg = lExporteur.ajouteChamp(lSigmaConst, "SigmaConst");
    }
  }

  std::cout << " no 13 "<<traducteurERMsg(lMsg) << endl;
  // Sigma en lin Elem
  ChampTensO2SymLinElem lSigmaLin;
  if(lMsg){
    lSigmaLin.asgnChampGeometrique(*lChampGeo);
    lSigmaLin.asgnNom("SigmaLin");
  }

  PPProjectionL2 lPPProjectionL2SigmaAnaLin;

  if (lMsg && lMecanik){
    lPPProjectionL2SigmaAnaLin.asgnChampGeometrique(*lChampGeo);
    lPPProjectionL2SigmaAnaLin.asgnParametres(0, lSchemaIntg,lSigmaAna,lSigmaLin);
    lProbMec.ajoutePostTraitement(lPPProjectionL2SigmaAnaLin); // Fait la mise a jour du tenseur de contraintes lissé (projection L2)
    lMsg = lExporteur.ajouteChamp(lSigmaLin, "Sigma");
  }

  ChampTensO2SymLinElem lSigmaLinM;
  if(lMsg){
    lSigmaLinM.asgnChampGeometrique(*lChampGeo);
    lSigmaLinM.asgnNom("SigmaLinM");
  }
  ChampContraintesHPP lSigmaAnaM("SigmaAnaM");
  PPProjectionL2 lPPProjectionL2SigmaAnaM;
  ChampVect3DContinu* lZero3D=0;
  if(lMsg) lMsg = lFichierChamp.reqChamp("Zero_3D",lZero3D);
  if (lMsg && lMecanik){
    lSigmaAnaM.asgnChampGeometrique(*lChampGeo);
    lSigmaAnaM.asgnChamps(*lZero3D,*lE);
    lSigmaAnaM.asgnCouplage(*lVraiM,*lVraiMP,*lBetaM);

    lPPProjectionL2SigmaAnaM.asgnChampGeometrique(*lChampGeo);
    lPPProjectionL2SigmaAnaM.asgnParametres(0, lSchemaIntg,lSigmaAnaM,lSigmaLinM);
    //lProbMec.ajoutePostTraitement(lPPProjectionL2SigmaAnaM); // Fait la mise a jour du tenseur de contraintes lissé (projection L2)
   // lMsg = lExporteur.ajouteChamp(lSigmaLinM, "SigmaM");
    //lMsg = lExporteur.ajouteChampCommeChampPrecalculePtsIntegration(lSigmaAnaM, "SigmaM",lSchemaIntg);
  }

  ChampTensO2SymLinElem lSigmaLinU;
  if(lMsg){
    lSigmaLinU.asgnChampGeometrique(*lChampGeo);
    lSigmaLinU.asgnNom("SigmaLinU");
  }
  ChampContraintesHPP lSigmaAnaU("SigmaAnaU");
  PPProjectionL2 lPPProjectionL2SigmaAnaU;
  if (lMsg && lMecanik){
    lSigmaAnaU.asgnChampGeometrique(*lChampGeo);
    lSigmaAnaU.asgnChamps(*lU,*lE);

    lPPProjectionL2SigmaAnaU.asgnChampGeometrique(*lChampGeo);
    lPPProjectionL2SigmaAnaU.asgnParametres(0, lSchemaIntg,lSigmaAnaU,lSigmaLinU);
    //lProbMec.ajoutePostTraitement(lPPProjectionL2SigmaAnaU); // Fait la mise a jour du tenseur de contraintes lissé (projection L2)
    //lMsg = lExporteur.ajouteChamp(lSigmaLinU, "SigmaU");
    //lMsg = lExporteur.ajouteChampCommeChampPrecalculePtsIntegration(lSigmaAnaU, "SigmaU",lSchemaIntg);
  }

  //
  // Evaluation des contraintes le long d'une ligne dans le domaine
  //
  //Vecteur3D lPointDepart(0.013,0.0,0.0121);
  //Vecteur3D lPointArrivee(0.013,0.0,0.007);
  //Entier lNbEval(250);
  //VectDyn<DReel> lVecteurEvaluations(lNbEval);
  //Entier lNbChiffres(3);
  //ChaineCar lNomSortie("SigmaZZ");
  //PPEvalueChampSurCourbe<ChampScalaire> lPPEvalueSurLigne;
  //lPPEvalueSurLigne.asgnParametres(lSigmaAna.reqComposante(2),lPointDepart,lPointArrivee,lNbEval,lVecteurEvaluations);
  //lPPEvalueSurLigne.afficheEvaluations(lNomSortie,lNbChiffres);
  //lProbMec.ajoutePostTraitement(lPPEvalueSurLigne);

  //
  // 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);
    // Bloc pour l'impression des propriétés des matériaux
    if(lMsg && 1==0){

      ExportGIREF* lExporteurPropMat1 = new ExportGIREF();
      lExporteurPropMat1 ->asgnExporteToujoursMaillage(false);
      lExporteurPropMat1->asgnChampGeometrique(*lChampGeo);

      lExporteurPropMat1->asgnDegreInterpolation(Exportation::eMAILLAGE_LINEAIRE);

      ChampScalaire* lVraia = 0;
      lMsg = lFichierChamp.reqChamp("a",lVraia);
      ChampScalLin la(*lChampGeo);
      if(lMsg){
      PPReinterpole lPPReinterpole0;
      lPPReinterpole0.asgnParametres(*lVraia,la);
      lMsg = lPPReinterpole0.effectueCalcul();
      }

      if(lMsg) lMsg = lExporteurPropMat1->ajouteChamp(la,"a");

      ChampScalLin lCMRho(*lChampGeo);
      PPReinterpole lPPReinterpole;
      lPPReinterpole.asgnParametres(*lC_MRho,lCMRho);
      lMsg = lPPReinterpole.effectueCalcul();

      if(lMsg) lMsg = lExporteurPropMat1->ajouteChamp(lCMRho,"CMRho");

      ChampTensO2SymLin lKM(*lChampGeo);
      PPReinterpole lPPReinterpole1;
      lPPReinterpole1.asgnParametres(*lK_M,lKM);
      lPPReinterpole1.effectueCalcul();

      if(lMsg) lMsg = lExporteurPropMat1->ajouteChamp(lKM,"K_M");

      ChampTensO2SymLin lBM(*lChampGeo);
      PPReinterpole lPPReinterpole2;
      lPPReinterpole2.asgnParametres(*lBetaM,lBM);
      if(lMsg) lMsg = lPPReinterpole2.effectueCalcul();

      if(lMsg) lMsg = lExporteurPropMat1->ajouteChamp(lBM, "BetaM");

      ChampScalaire* lELAna = 0;
      if(lMsg) lMsg = lFichierChamp.reqChamp("EL",lELAna);
      ChampScalLin lEL(*lChampGeo);
      lPPReinterpole.asgnParametres(*lELAna,lEL);
      lMsg = lPPReinterpole.effectueCalcul();

      if(lMsg) lMsg = lExporteurPropMat1->ajouteChamp(lEL,"EL");

      ChampScalaire* lETAna = 0;
      if(lMsg) lMsg = lFichierChamp.reqChamp("ET",lETAna);
      ChampScalLin lET(*lChampGeo);
      lPPReinterpole.asgnParametres(*lETAna,lET);
      lMsg = lPPReinterpole.effectueCalcul();

      if(lMsg) lMsg = lExporteurPropMat1->ajouteChamp(lET,"ET");

      ChampScalaire* lERAna = 0;
      if(lMsg) lMsg = lFichierChamp.reqChamp("ER",lERAna);
      ChampScalLin lER(*lChampGeo);
      lPPReinterpole.asgnParametres(*lERAna,lER);
      lMsg = lPPReinterpole.effectueCalcul();

      if(lMsg) lMsg = lExporteurPropMat1->ajouteChamp(lER,"ER");


      ChampScalLinElem lE0000(*lChampGeo);
      lPPReinterpole.asgnParametres(lE->reqComposante(0),lE0000);
      if(lMsg) lMsg = lPPReinterpole.effectueCalcul();

      if(lMsg) lMsg = lExporteurPropMat1->ajouteChamp(lE0000,"E0000");

      ChampScalLinElem lE1111(*lChampGeo);
      lPPReinterpole.asgnParametres(lE->reqComposante(5),lE1111);
      if(lMsg) lMsg = lPPReinterpole.effectueCalcul();

      if(lMsg) lMsg = lExporteurPropMat1->ajouteChamp(lE1111,"E1111");

      ChampScalLinElem lE2222(*lChampGeo);
      lPPReinterpole.asgnParametres(lE->reqComposante(20),lE2222);
      if(lMsg) lMsg = lPPReinterpole.effectueCalcul();

      if(lMsg) lMsg = lExporteurPropMat1->ajouteChamp(lE2222,"E2222");

      if(lMsg) lMsg = lExporteurPropMat1->exporteInfo("valeurs_finales_materiaux");
    }

  }

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

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

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

  std::cout << " no 15 "<<traducteurERMsg(lMsg) << endl;
  //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;

  pStringRetour = "";

  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;
}