//************************************************************************
// --- 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.
//
//  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 "ChampScalaire.h"
#include "ChampScalContinu.h"
#include "ChampScalAConstant.h"
#include "ChampTensO2SymContinu.h"
#include "ChampTensorielO4Sym.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 "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 "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;

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 *lMInitiale = 0;          // humidite initiale
  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;                // les déplacements
//  ChampContraintesHPP* lSigmaU=0;            // les contraintes associées à lU
  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
  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)
  SolveurLinPETSc       *lSolveurMec;   // 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

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

  // *****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) {
    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>("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);

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

    if(lMsg){
      lMsg = lFichierChamp.reqChamp<ChampScalContinu>("MInitiale",lMInitiale,false);
      if(lMsg && lMInitiale == 0) lMInitiale = lM->newChampCopie();
    }

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

    // 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) 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
  //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);
  lNbTempsImpose = 0;

  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.;
    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;
  TFSigmaUEpsilonV3D        lTFRigidite;

  // contact
  TFContactV3DUxNVScal  lTFContactGCP;
  TFContactV3DUxNVScal  lTFContactGCP2;

  //On assigne les champs aux termes de formulation.
  if(lMsg) lMsg = lTFDiffusionM. asgnChamps(*lChampGeo, *lM, *lM, *lK_M,    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,    lSchemaIntg);
    if(lMsg) lMsg = lTFDiffusionT.  asgnChamps(*lChampGeo, *lT, *lT, *lK_T,    lSchemaIntg);
    if(lMsg) lMsg = lTFMasseT.      asgnChamps(*lChampGeo, *lT, *lT, *lC_TRho,lSchemaIntg);
  }

  if(lMecanik && lMsg){

    lMsg = lTFRigidite.asgnChamps(*lChampGeo, *lU, *lU, *lE, lSchemaIntg);

    if(lMsg) lMsg = lTFCouplageHydM.asgnChamps(*lChampGeo, *lM, *lMInitiale, *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);
    }

    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){
      lProbMec.ajouteTF    (lTFRigidite);
      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>("SolveurMec",lSolveurMec);
  // On assigne les temps ou l'on force le calcul
  if(lMsg && lTemps) lSolveurHyd->asgnTempsImposes(lTempsImposes);

  SolveurGCPAvecCorpsRigide      *lSolveurGCP = 0;

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

  EntiteGeometrique* lEntite = 0;
  if(lMsg && lContakt) {
    if(lMsg && lCorpsParGCP) lEntite =  &(lListeEntite.reqEntiteGeometrique(lNomEntiteControle));

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

  bool lVrai(true);

  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
  ExportGIREF lExportGIREF;
  if(lPasImpression && lMsg){
    lExportGIREF.asgnChampGeometrique(*lChampGeo);
    lExportGIREF.asgnExporteMaillage(false);
    lProbHyd.ajouteExportation(lExportGIREF, lPasImpression, lNomFichierChamp);
    if(lMecanik){
      lMsg = lExportGIREF.ajouteChamp(*lU, "U");
    }
  }

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

  // 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(8);
   lFichierDT << std::scientific;

   PPAfficheValeur <DReel> lAfficheTemps;
   DReel *lDeltaT = 0;
   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);



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

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

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