// ************************************************************************
// --- 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.
// ************************************************************************

//
//  Résolution du problème:
//
//  (1) température
//  C_TRho dT/dt - div(K_T grad(T)) = 0
//
//  (2) teneur en humidité
//  C_MRho dM/dt - div(K_M grad(M) + K_MT grad(T) ) = 0
//
//  (3) déplacements
//  -div(E:(epsilon(U) - BetaM*(M-MInitiale) - BetaT)) = F
//
//  Les équations (1) et (3) sont optionnelles auxquels cas leur
//  traitement est complètement ignorées dans le code. Dans (3) le
//  membre de droite F est optionnel.
//
//  U est obtenu de manière incrémentale i.e. par correction à chaque
//  pas de temps. Il faut donc dans les fichiers champs imposé les CL,
//  CI et efforts incrémentaux!!!
//
//  USAGE: ThermoHygroMecanik.[dev,opt] préfixe_des_fichiers_d_entrée

#include "ChaineCar.h"
#include "ChampContraintesHPP.h"
#include "ChampGeometrique.h"
#include "ChampScalAConstant.h"
#include "ChampScalConstElem.h"
#include "ChampScalaire.h"
#include "ChampTensO2SymContinu.h"
#include "ChampTensorielO4Sym.h"
#include "ChampVect3DContinu.h"
#include "ExportVTK.h"
#include "Geometrie.h"
#include "GestionFichierChamps.h"
#include "InitialisationArgvArgc.h"
#include "InitialisationTraitementSignal.h"
#include "ListEntitesGeometrique.h"
#include "ListeConditionsLimites.h"
#include "MEFPPUtil.h"
#include "Maillage.h"
#include "OptionsSolveurLinPETSc.h"
#include "PETScInitialisation.h"
#include "PPAfficheValeur.h"
#include "PPCalculChampaXpbY.h"
#include "PPCallBack.h"
#include "PPCallBack2Args.h"
#include "PPEvalueChampSurCourbe.h"
#include "PPEvalueChampXYZ.h"
#include "PPIntegre.h"
#include "PPProjectionL2.h"
#include "PPReinterpole.h"
#include "PPResolutionProbleme.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 "SolveurGCP.h"
#include "SolveurGCPAvecCorpsRigide.h"
#include "SolveurInstNlinPETSc.h"
#include "SolveurLinPETSc.h"
#include "TFContraintesInitiales.h"
#include "TFCouplageAvecFctsConnues.h"
#include "TFDiffusion.h"
#include "TFDiffusionNewton.h"
#include "TFMasseGenScal.h"
#include "TFMasseGenV3D.h"
#include "TFSigmaUEpsilonV3D.h"
#include "TFSigmaUEpsilonV3DIncrementalVieillissant.h"
#include "TFSourceFVVect3D.h"
#include "TenseurO4Sym.h"
#include "VerifieElementsNormalises.h"
#include "fonctionsChaineCar.h"
#include "fonctionsTenseurs.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 récupère le GroupeProcessus associé à 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, lNomSortie;

  // 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_sortie (optionnel)";
    lMsg.ajoute(ERMsg(ERMsg::ERREUR, lMessageUsage));
  }
  else {
    lNomFichierChamp = lArgv[1];
    lNomSortie = lNomFichierChamp;
    if (lNbArg > 2) lNomSortie = lArgv[2];
  }

  lNomMaillage = lNomFichierChamp;

  // Pré-traitement des données du fichier champs:
  if (lMsg) {
    const EntierN lTailleBuffer = 65536;
    char *lLigne = new char[lTailleBuffer];
    ifstream lFichierChamp((lNomFichierChamp + 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;
          }
        }
        lFichierChamp.getline(lLigne,lTailleBuffer);
      }
    }
    lFichierChamp.close();
  }
  // Fin du pré-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;

  Maillage lMail;

  // Création des champs requis pour le problème:
  // champ géometrique du domaine entier
  ChampGeometrique *lChampGeo=0;
  // humidité
  ChampScalContinu *lM = 0;
  // humidité au pas précedent
  ChampScalContinu *lMInitiale = 0;
  // humidité au pas précedent
  ChampScalContinu *lMP = 0;
  // dérivée en temps de M
  ChampScalContinu *lMt = 0;
  // densité * "chaleur spécifique"
  ChampScalContinu *lC_MRho = 0;
  // diffusivité effective d_b*D
  ChampTensO2SymContinu *lK_M = 0;
  // derivée de la diffusivité effective en M
  ChampTensO2SymContinu *ldK_MdM = 0;
  // thermo-diffusion
  ChampTensO2SymContinu *lK_MT = 0;
  // diffusion thermique
  ChampTensO2SymContinu *lK_T = 0;
  // chaleur spécifique multipliée par la densité
  ChampScalContinu *lC_TRho = 0;
  // température
  ChampScalContinu *lT = 0;
  // l'incrément des déplacements
  ChampVect3DContinu *lU = 0;
  // les déplacements
  ChampVect3DContinu *lDeplacements = 0;
  // tenseur des coefficients d'élasticité
  ChampTensorielO4Sym *lE = 0;
  // tenseur des coefficients d'élasticité au pas de temps précédent
  ChampTensorielO4Sym *lEP = 0;
  // tenseur des dilatations hydriques
  ChampTensO2SymContinu *lBetaM = 0;
  // tenseur des dilatations hydriques  au pas de temps précédent
  ChampTensO2SymContinu *lBetaMP = 0;
  // tenseur des dilatations thermiques
  ChampTensO2SymContinu *lBetaT = 0;
  // les charges volumiques
  ChampVectoriel3D *lF = 0;

  // Variables lues dans le fichier:
  // le problème a une composante thermique
  bool lThermik(false);
  // le problème a une composante mécanique
  bool lMecanik(true);
  // formulation mécanique est incrémentale
  bool lFormulationIncrementale(true);
  // On résout avec le terme dK_M/DM (méthode de Newton au lieu de pt
  // fixe)
  bool lEmploiNewtonEnM(false);
  // le tenseur de dilatation hydrique tiens compte d'effet d'hystérèse
  bool lHysterese(true);
  // fréquence du calcul de mécanique: 0 -> pas de calcul
  Entier lPasDeCalcul;
  // fréquence d'impression des résultats
  Entier lPasImpression;
  // Nombre de points du temps ou on force le calcul
  Entier lNbTempsImpose;
  // le solveur instationnaire pour M (et T possiblement)
  SolveurInstNlinPETSc *lSolveurHyd = 0;
  // le solveur stationnaire pour U
  SolveurLinPETSc *lSolveurMec = 0;

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

  bool lVrai(true), lFaux(false);

  Entier lPoints1D(4);

  SchemaIntg lSchemaIntg;

  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 déclare les problèmes "vides"
  ProblemeEFInstationnaire lProbHyd;
  if (lMsg)
    lProbHyd.asgnNom("Composante hydrique/thermique");
  ProblemeEF lProbMec;
  if (lMsg)
    lProbMec.asgnNom("Composante mécanique");

  if (lMsg) {
    lProbHyd.asgnGroupeProcessus(lGroupeGlobal);
    lProbMec.asgnGroupeProcessus(lGroupeGlobal);
  }

  if (lMsg) {
    lMsg = lProbHyd.lisDonnees(lNomMaillage,lMail, lGeo,lListeEntite);
  }

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

  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
  if (lMsg) lMsg = lFichierChamp.reqBooleen("Thermique",lThermik,false);
  if (lMsg) lMsg = lFichierChamp.reqBooleen("Mecanique",lMecanik,false);

  // On demande le champ de transformation géométrique au fichier de champs
  if (lMsg) lMsg = lFichierChamp.reqChamp("ChampGeo",lChampGeo);
  if (lMsg) lProbHyd.asgnChampGeometrique(*lChampGeo);

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

  PPProjectionL2 lPPProjectionL2;
  if (lMsg) lPPProjectionL2.asgnChampGeometrique(*lChampGeo);

  PPCalculChampaXpbY<ChampTensO2SymLin> lMAjSigma;
  ChampContraintesHPP lDeltaSigma;

  ChampTensO2SymLin lDeltaSigmaLin;
  if (lMsg) lDeltaSigmaLin.asgnChampGeometrique(*lChampGeo);

  ChampTensO2SymLin lSigma;
  if (lMsg) lSigma.asgnChampGeometrique(*lChampGeo);

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

  // On demande un champ qui se nomme "M_t" pour la valeur différence
  // finie de M (hystérèse présente dans le problème)
  if (lMsg) lMsg = lFichierChamp.reqChamp("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("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("K_M",lK_M);

  // On demande un champ qui se nomme "dK_MdM" pour la derivée du
  // tenseur de "diffusivité effective" en M
  if (lMsg) lMsg = lFichierChamp.reqChamp("dK_MdM",ldK_MdM,false);
  if (lMsg && ldK_MdM) lEmploiNewtonEnM = true;

  // ChampTensorielO2NSym *lTransformation = 0;
  // if (lMsg) lMsg = lFichierChamp.reqChamp("transformation",lTransformation);

  // 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("K_MT",lK_MT);

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

    // On demande un champ qui se nomme "K_T" pour la "diffusion thermique"
    if (lMsg) lMsg = lFichierChamp.reqChamp("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("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 résultats de mécanique alors je ne les fait pas!"<<std::endl;
    lMecanik = false;
  }

  // On demande une constante qui se nomme Epsilon pour la matrice
  // masse Epsilon*M
  DReel lEpsilon = 0;
  if (lMsg)
    lMsg = lFichierChamp.reqScalaireConstant("Epsilon",lEpsilon,false);

  ChampScalAConstant* lUn = 0;
  if (lMecanik && lMsg) {
    // On demande si la formluation en déplacements est incrémentale
    if (lMsg)
      lMsg = lFichierChamp.reqBooleen("Formulation_incrementale",lFormulationIncrementale);
    // On demande un champ qui se nomme "U" pour le déplacement
    // engendré par l'humidité.
    if (lMsg)
      lMsg = lFichierChamp.reqChamp("U",lU);
    if (lMsg) {
      if (lFormulationIncrementale) {
        lMsg = lFichierChamp.reqChamp("Deplacements",lDeplacements,false);
        if (lMsg && !lDeplacements)
          lDeplacements = lU->newChampCopie();
      } else {
        lDeplacements = lU;
      }
    }

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

    // Un champ "EijklP" est possible si la formulation est incrémentale
    if (lFormulationIncrementale && lMsg)
      lMsg = lFichierChamp.reqChamp("EijklP",lEP,false);
    // if (lMsg) lMsg = lFichierChamp.reqChamp("EL",lEL);
    // if (lMsg) lMsg = lFichierChamp.reqChamp("ER",lER);
    // if (lMsg) lMsg = lFichierChamp.reqChamp("ET",lET);
    // if (lMsg) lMsg = lFichierChamp.reqChamp("GLR",lGLR);
    // if (lMsg) lMsg = lFichierChamp.reqChamp("GLT",lGLT);
    // if (lMsg) lMsg = lFichierChamp.reqChamp("GRT",lGRT);
    // if (lMsg) lMsg = lFichierChamp.reqChamp("VTR",lVTR);
    // if (lMsg) lMsg = lFichierChamp.reqChamp("VLR",lVLR);
    // if (lMsg) lMsg = lFichierChamp.reqChamp("VLT",lVLT);

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

    // Un champ "BetaP" est possible si la formulation est incrémentale
    if (lFormulationIncrementale && lMsg)
      lMsg = lFichierChamp.reqChamp("BetaMP",lBetaMP,false);

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

    // On demande un champs qui se nomme "BetaT" pour le tenseur de
    // dilatation thermique
    if (lMsg && lThermik) lMsg = lFichierChamp.reqChamp("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("F",lF,false);

  }

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

  // On permet d'imposer des points dans le temps ou le calcul doit être 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 imposés 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 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 négatif 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 lNombrePointsSurCourbe(100);
  if (lMsg) lMsg = lFichierChamp.reqScalaireConstant("NbPointsSurCourbe",lNombrePointsSurCourbe,false);

  Entier lNombreEval(0);
  if (lMsg) lMsg = lFichierChamp.reqScalaireConstant("NbEval",lNombreEval,false);
  if (lNombreEval < 0 && lMsg) lNombreEval = 0;
  if (lNombreEval > 0 && lMsg) lCoordonnees.asgnLongueur(lNombreEval);
  if (lNombreEval > 0 && lMsg) {
    ChaineCar lCoor = "Coord_";
    for (Entier i = 0; i < lNombreEval && lMsg; ++i) {
      ChaineCar lCoordi = ajouteNumeroALaFin(lCoor,i,3);
      lMsg = lFichierChamp.reqVecteur3DConstant(lCoordi,lCoordonnees[i]);
    }
  }

  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 + string("_hyd");
  ChaineCar lNomMec = lNomMaillage + 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 lTFCouplageHydT;
  TFSourceFVVect3D lTFChargeVolumique;
  TFMasseGenV3D lTFMasseGenV3D;

  // mécanique formulation non incrementale
  TFCouplageAvecFctsConnues lTFCouplageHydM;
  TFSigmaUEpsilonV3D lTFRigidite;

  // mécanique formulation incrementale
  TFSigmaUEpsilonV3DIncrementalVieillissant lTFRigiditeIncremental;

  // 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 (!lFormulationIncrementale) {
      lMsg = lTFRigidite.asgnChamps(*lChampGeo,*lU,*lU,*lE,lSchemaIntg);
      if (lMsg)
        lMsg = lTFCouplageHydM.asgnChamps(*lChampGeo,*lM,*lMP,*lU,*lE,*lBetaM,lSchemaIntg);
    } else {
      // Est-ce que E dépend du temps?
      if(lEP) {
        lMsg = lTFRigiditeIncremental.asgnChamps(*lChampGeo,*lU,*lU,*lDeplacements,*lE,*lEP,lSchemaIntg);
      } else {
        lMsg = lTFRigiditeIncremental.asgnChamps(*lChampGeo,*lU,*lU,*lDeplacements,*lE,lSchemaIntg);
      }
      if (lMsg){
        lMInitiale = lM->newChampCopie();
        // Estc-ce que Beta dépend du temps ?
        if(lBetaMP){
          lTFRigiditeIncremental.asgnCouplage(*lM,*lMP,*lMInitiale,*lBetaM,*lBetaMP);
        } else {
          lTFRigiditeIncremental.asgnCouplage(*lM,*lMP,*lMInitiale,*lBetaM);
        }
      }
    }

    if (lMsg)
      lMsg = lTFMasseGenV3D.asgnChamps(*lChampGeo,*lU,*lU,lEpsilon,lSchemaIntg,false);
    if (lMsg && lBetaT)
      lMsg = lTFCouplageHydT.asgnChamps(*lChampGeo,*lUn,*lU,*lE,*lBetaT,lSchemaIntg);
    if (lMsg && lF)
      lMsg = lTFChargeVolumique.asgnChamps(*lChampGeo,*lF,*lU,lSchemaIntg);
  }

  if (lMsg) {
    // Ajout des termes de formulations au 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 (!lFormulationIncrementale){
        lProbMec.ajouteTF(lTFRigidite);
        lProbMec.ajouteTF(lTFCouplageHydM);
      }else {
        lProbMec.ajouteTF(lTFRigiditeIncremental);
      }

      lProbMec.ajouteTF(lTFMasseGenV3D);
      if (lF)
        lProbMec.ajouteTF(lTFChargeVolumique);
      if (lBetaT)
        lProbMec.ajouteTF(lTFCouplageHydT);
    }
  }

  // Lecture et initialisation des conditions limites
  if (lMsg)
    lMsg = lProbHyd.lisEtInitialiseConditionsLimites(lNomHyd,lListeCondLimHyd);
  if (lMsg && lMecanik) {
    lMsg = lProbMec.lisEtInitialiseConditionsLimites(lNomMec,lListeCondLimMec);
  }

  // Lecture des solveurs
  if (lMsg) lMsg = lFichierChamp.reqSolveur<SolveurInstNlinPETSc>("SolveurInst",lSolveurHyd);

  SolveurLinPETSc* lSolveurTempo=0;
  if (lMsg) lMsg = lFichierChamp.reqSolveur<SolveurLinPETSc>("SolveurLinmecanique",lSolveurTempo);
  lSolveurMec = lSolveurTempo;

  // On assigne les temps où l'on force le calcul
  if (lMsg && lTemps) lMsg = lSolveurHyd->asgnTempsImposes(lTempsImposes);

  if (lMecanik && lMsg) {
    // A priori la matrice de rigidité dépend 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 mécanique
    Entier lEntierBidon(1);
    if (lMsg) lMsg = lFichierChamp.reqScalaireConstant("NbPasDeTemps",lEntierBidon,false);
    if (lMsg) lSolveurMec->asgnLongueurNumeroPourExport(lEntierBidon+lNbTempsImpose);
  }

  ExportVTK lExport;
  if (lPasImpression && lMsg) {
    lExport.asgnChampGeometrique(*lChampGeo);
    lProbHyd.ajouteExportation(lExport,lPasImpression,lNomSortie);
    if (lMecanik) {
      if (lMsg) lMsg = lExport.ajouteChamp(*lDeplacements,"U");
    }
  }

  // Début des pre/post traitements

  // lPPReqTempsHyd       Pour savoir la valeur du temps (pour la
  //                      passer au problème mécanique) on demande le
  //                      temps au problème hydrique.

  // lPPAsgnTempsMec      Le problème mécanique est en post-traitement
  //                      de l'hydrique alors il nous faut lui passer
  //                      la valeur du temps au moment du calcul.

  PPTraitementGen* lPPReqTempsHyd = 0;
  PPTraitementGen* lPPAsgnTempsMec = 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));

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

  PPEvalueChampSurCourbe<ChampScalContinu> lEvalueMSurCourbe;
  VectDyn<DReel> lMAuxPointsDelaCourbe;

  // PPEvalueChampSurCourbe<ChampScalContinu> lEvalueMSurCourbe1;
  VectDyn<DReel> lMAuxPointsDelaCourbe1;

  // Points d'intérêt dans le domaine:

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

  // pour M moyenne
  DReel lIntegraleM = 0;
  DReel lMmoyen = 0;

  PPIntegre<ChampScalContinu, std::identity<DReel> > lEvalueMmoyen;
  if (lMsg)
    lMsg = lEvalueMmoyen.asgnParametres(*lM,std::identity<DReel>(),0,lIntegraleM,lMmoyen,false,false,lSchemaIntg,0);

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

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

  // pour l'écriture des évaluations ponctuelles:
  ofstream lFichierM;

  PPAfficheValeur<DReel> lAfficheMmoyen;
  ofstream lFichierMmoyen;

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

  ofstream lFichierU;

  // On a donné dans le fichier champs des points où on veut une évaluation:
  if (lMsg && lNombreEval) {
    lFichierM.open((lNomSortie + ".M").c_str());
    lEvalueM.asgnParametres(*lChampGeo,*lM,lCoordonnees,lMAuxPointsDonnes);
    lEvalueM.afficheEvaluations(lFichierM);
  }

  if (lMsg && lNombreEval) {
    lFichierMmoyen.open((lNomSortie + ".Mmoyen").c_str());
    lAfficheMmoyen.afficheLeTemps(lProbHyd,lVrai);
    lAfficheMmoyen.asgnParametres("Mmoyen",lMmoyen,lFichierMmoyen);
    lFichierMmoyen << setprecision(8);
    lFichierMmoyen << scientific;
  }

  if (lMecanik && lNombreEval && lMsg) {
    lFichierU.open((lNomSortie + ".U").c_str());
    lEvalueU.asgnParametres(*lChampGeo,*lDeplacements,lCoordonnees,lUAuxPointsDonnes);
    lEvalueU.afficheEvaluations(lFichierU);
  }

  // Le problème mécanique est résolu en post-traitement du problème
  // hydrique:
  PPResolutionProbleme lPPResolution;
  if (lMsg && lMecanik)
    lPPResolution.asgnParametres(lProbMec,*lSolveurMec);

  // En cas d'hystérèse on doit conserver la valeur de M au pas précédent
  PPReinterpole lStockeMAuPasPrecedent;
  if (lMecanik && lFormulationIncrementale && lMsg)
    lMsg = lStockeMAuPasPrecedent.asgnParametres(*lM,*lMP);

  // Ajout de l'incrément U aux déplacements
  PPCalculChampaXpbY<ChampVect3DContinu> lMAJU;
  DReel lUnReel(1.);
  if (lMecanik && lFormulationIncrementale && lMsg)
    lMAJU.asgnParametres(lUnReel,lUnReel,*lDeplacements,*lU,*lDeplacements);
  if (lMsg)
    lPPProjectionL2.asgnParametres(&lGeo,lSchemaIntg,lDeltaSigma,lDeltaSigmaLin);
  if (lMsg)
    lMAjSigma.asgnParametres(lUnReel,lUnReel,lSigma,lDeltaSigmaLin,lSigma);
  if (lMsg) {
    lDeltaSigma.asgnChamps(*lU,*lE);
    lDeltaSigma.asgnCouplage(*lM,*lMP,*lBetaM);
  }

  // On ajoute les pré/post traitements ATTENTION la chronologie est
  // importante!

  if (lMecanik && lMsg) {
    // on calcule U après M:
    lProbHyd.ajoutePostTraitement(*lPPReqTempsHyd);
    lProbHyd.ajoutePostTraitement(*lPPAsgnTempsMec);
    // lProbHyd.ajoutePostTraitement(lPPContraintesEL3D);
    lProbHyd.ajoutePostTraitement(lPPResolution,lPasDeCalcul);

    // Ajoute l'incrément aux déplacements
    if (lFormulationIncrementale)
      lProbMec.ajoutePostTraitement(lMAJU);
    lProbMec.ajoutePostTraitement(*lPPAsgnTempsDeplacements);

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

  if (lMsg && lNombreEval>0) {
    lProbHyd.ajoutePostTraitement(lEvalueM);
    lProbHyd.ajoutePostTraitement(lEvalueMmoyen);
    lProbHyd.ajoutePostTraitement(lAfficheMmoyen);
    if (lThermik) {
      lProbHyd.ajoutePostTraitement(lEvalueT);
      lProbHyd.ajoutePostTraitement(lAfficheT);
    }
    if (lMsg && lNombrePointsSurCourbe > 0 )
      lProbHyd.ajoutePostTraitement(lEvalueMSurCourbe);
  }

  // Mise a jour de MP pour le "prochain" pas dans le cas avec hystérèse
  if (lMsg && lMecanik && lBetaM && lFormulationIncrementale)
    lProbMec.ajoutePostTraitement(lStockeMAuPasPrecedent);

  // Résolution du problème en T, M et en U
  if (lMsg) lMsg = lSolveurHyd->resoudre(lProbHyd);

  if (lNombreEval>0) {
    lFichierM.close();
    lFichierMmoyen.close();
  }
  if (lMecanik && lNombreEval>0) lFichierU.close();
  if (lThermik && lNombreEval>0) lFichierT.close();

  // Traitement d'erreur: si une erreur s'est produite, on l'affiche.
  if (!lMsg) {
    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;
}