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

//********************************************************************
//  Nom du fichier: 
//  
//  USAGE:<pre> bois.[dev,opt]  préfixe_des_fichiers_d_entrée
// 
//  Exemple:  
//  
//
//********************************************************************

#include "MEFPPUtil.h"
#include "ChaineCar.h"
#include "ChampContraintesHPP.h"
#include "ChampGeometrique.h"
#include "ChampPrecalculeContinuPtsIntegration.h"
#include "ChampScalaire.h"
#include "ChampTensorielO2Sym.h"
#include "ChampTensO2SymLin.h"
#include "ChampTensO2SymQuad.h"
#include "ChampTensO4SymContinu.h"
#include "ChampSolideHookeenEmpilement.h"
#include "ChampScalQuad.h"
#include "ChampVect3DQuad.h"
#include "ChampScalConstElem.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 "OptionsSolveurLinPETSc.h"
#include "PETScInitialisation.h"
#include "PPAfficheValeur.h"
#include "PPCallBack.h"
#include "PPEvalueChampXYZ.h"
#include "PPDiffFinieChamp.h"
#include "PPContraintesEL3D.h"
#include "PPResolutionProbleme.h"
#include "PPReinterpole.h"
#include "PPCalculChampaXpbY.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 "TFDiffusion.h"
#include "TFMasseGenScal.h"
#include "TFContraintesInitiales.h"
#include "TFCouplageAvecFctsConnues.h"
#include "TFSigmaUEpsilonV3D.h"
#include "traducteurERMsg.h"
#include "VecteurPETSc.h"
#include "ChampScalCrouRav.h"
#include "ChampVect3DCrouRav.h"
#include "PPNormeL2.h"
#include <fstream>

using namespace std;

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

  // 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_d_entree";
    lMsg.ajoute(ERMsg(ERMsg::ERREUR, lMessageUsage));
  }
  else {
    lNom          = lArgv[1];
  }

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

  SchemaIntg     lSchemaIntg;

  bool lCalculMec(true); // si on veut faire le calcul des deplacements
  //lCalculMec=false; // si on ne veut pas faire le calcul des deplacements

  Entier lTrois(3);
  
  if (lMsg) {
    // Création du schéma d'intégration des termes de formulation (avec valeurs par défaut)
    // Attention, on ne veut peut-être pas utiliser le schéma par défaut.
    
    lSchemaIntg.asgnIntegration(PtIntgTriangle(PtIntgTriangle::DouzePtsInternes));
    lSchemaIntg.asgnIntegrationFace3Sommets(PtIntgTriangle(PtIntgTriangle::DouzePtsInternes));
    lSchemaIntg.asgnIntegrationFace4Sommets(PtIntgQuadrangle(PtIntgElem1D(lTrois),PtIntgElem1D(lTrois)));
    
    lSchemaIntg.asgnIntegration(PtIntgElem1D(lTrois));
    lSchemaIntg.asgnIntegrationArete(PtIntgElem1D(lTrois));

    //     lSchemaIntg.asgnIntegration(PtIntgTetraedre(6));
    //     lSchemaIntg.asgnIntegration(PtIntgQuadrangle(PtIntgElem1D(3), PtIntgElem1D(3)));
    lSchemaIntg.asgnIntegration(PtIntgHexaedre(PtIntgElem1D(lTrois), PtIntgElem1D(lTrois), PtIntgElem1D(lTrois)));
    lSchemaIntg.asgnIntegration(PtIntgPrismeTri(PtIntgTriangle(PtIntgTriangle::DouzePtsInternes), PtIntgElem1D(3)));

  }

  //On déclare les problèmes "vides"
  ProblemeEFInstationnaire lProbHyd;
  ProblemeEF lProbMec;
  lProbHyd.asgnGroupeProcessus(lGroupeGlobal);
  lProbMec.asgnGroupeProcessus(lGroupeGlobal);
  lProbHyd.lisDonnees(lNom,lMail,lGeo,lListeEntite);
  //Déclaration de l'objet qui fait la gestion des champs pour le problème.
  GestionFichierChamps lFichierChamp;

  lFichierChamp.asgnDonneesBase(&lMail);
  lFichierChamp.asgnDonneesGeometriques(&lGeo,&lListeEntite);

//   // le champs du tenseur des contraintes courant
//   ChampContraintesHPP lSigmaU0("SigmaU0");
//   if(lMsg) lMsg = lFichierChamp.asgnChamp(&lSigmaU0, TypesDonneesMEFPP::ChampTensO2SymContinu);	
//   // le champs du tenseur des contraintes au pas de temps precedent
//   ChampContraintesHPP lSigmaU1("SigmaU1");
//   if(lMsg) lMsg = lFichierChamp.asgnChamp(&lSigmaU1, TypesDonneesMEFPP::ChampTensO2SymContinu);	


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

  // *****Requête des données du fichier .champs nécessaires au calcul*****

  // On demande le champ de transformation géométrique au fichier de champs
  ChampGeometrique *lChampGeo=0;
  if (lMsg) {
    lMsg = lFichierChamp.reqChamp<ChampGeometrique>("ChampGeo",lChampGeo);
  }
  if(lMsg) lProbHyd.asgnChampGeometrique(*lChampGeo);
  // Assignation de la géométrie, entité, maillage, champgéométrique  au problème
  if(lCalculMec && lMsg) lProbMec.asgnDonnees(lMail,lGeo,*lChampGeo,lListeEntite);

  //
  // Pour la partie hydrique:
  //

  //On demande  un champ qui se nomme "M" pour l'humidité
  ChampScalContinu *lM = 0;
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("M",lM);
  //
  ChampScalContinu *lMPrec = 0;
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("MPrec",lMPrec);
  ChampScalContinu *lMInit = 0;
  if(lMsg){
    lMInit = lM->newChampCopie();
  }

  //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)
  ChampScalContinu *lMt = 0;
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("M_t",lMt,false);
//   ChampScalContinu *lMtPrec = 0;
//   if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("M_tPrec",lMtPrec,false);

  // le calcul de la difference finie depend de l'interpolation de M
  bool lMtQuad(true);
  bool lMtLin(false);
  if(lMsg) lMsg = lFichierChamp.reqBooleen("MtQuad",lMtQuad,false);
  if(lMsg) lMsg = lFichierChamp.reqBooleen("MtLin",lMtLin,false);

  //On demande  un champ qui se nomme "d_b" pour la densité de base divisée par 100 (kg m^3 / 100)
  ChampScalContinu *ld_b = 0;
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampScalContinu>("densite",ld_b);

  //On demande  un champ qui se nomme "K_M" pour le tenseur de "diffusivité effective" (d_b * D)
  ChampTensO2SymContinu *lK = 0;
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("K_M",lK);
//   ChampTensO2SymContinu *lKPrec = 0;
//   if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("K_MPrec",lKPrec);

  //
  // Pour la partie mécanique:
  //

    std::cerr << traducteurERMsg(lMsg) << endl;
  //On demande  un champ qui se nomme "U" pour le deplacement engendré par l'humidité.
  ChampVect3DContinu      *lU     = 0;
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampVect3DContinu>("U",lU);
  ChampVect3DContinu      *lUPrec = 0;
  if(lMsg) lUPrec = lU->newChampCopie();

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

  // On demande un champs qui se nomme "Beta" pour le tenseur de dilatation
  ChampTensO2SymContinu   *lBeta     = 0;
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("Beta",lBeta);
//   ChampTensO2SymContinu   *lBetaPrec = 0;
//   if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("BetaPrec",lBetaPrec);

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

//   //
//   // Preparation des contraintes initiales
//   //
  // Champ des contraintes 
  ChampContraintesHPP lSigmaU0("SigmaU0");
  if(lMsg){
    lSigmaU0.asgnChamps(*lUPrec,*lE);
    lSigmaU0.asgnChampGeometrique(*lChampGeo);
    lSigmaU0 *= -1.;
  }
//   // Champ des contraintes au temps n-1
//   if(lMsg){
//     lSigmaU1.asgnChamps(*lUPrec,*lEPrec);
//     lSigmaU1.asgnCouplage(*lMPrec,*lMInit,*lBetaPrec);
//     lSigmaU1.asgnChampGeometrique(*lChampGeo);
//   }
 
//   // Pour finir on lit le champs faisant la combinaison lDeltaSigma = SigmaU0 - SigmaU1
//   ChampTensO2SymContinu   *lDeltaSigma     = 0;
//   if(lMsg) lMsg = lFichierChamp.reqChamp<ChampTensO2SymContinu>("DeltaSigma",lDeltaSigma);

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

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

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

  DReel lNbreDiscontinuiteEnTemps(0);
  if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("NbTempsImpose",lNbreDiscontinuiteEnTemps,false);
  Entier lNbTempsImpose = (Entier) lNbreDiscontinuiteEnTemps;
  bool lTemps(lNbTempsImpose);

  VectDyn<DReel> lTempsImposes(lNbTempsImpose,0.);
  if(lTemps){
    ChaineCar lTImpose= "TempsImpose_";
    for (Entier i=0; i<lNbTempsImpose && lMsg; ++i){
      ChaineCar lTempsImposei = ajouteNumeroALaFin(lTImpose,i,1);
      cout<<lTempsImposei<<endl;
      lMsg = lFichierChamp.reqScalaireConstant(lTempsImposei,lTempsImposes[i]);
    }
  }

  // On demande les proprietes de la matrice masse
  bool lMasseDependDuTemps(true);
  bool lMasseDependDeLaSolution(true);
  if(lMsg) lMsg = lFichierChamp.reqBooleen("MDependTemps",lMasseDependDuTemps,false);
  if(lMsg) lMsg = lFichierChamp.reqBooleen("MDependSolution",lMasseDependDeLaSolution,false);

  // On demande les proprietes de la matrice de rigidité
  bool lRigiditeDependDuTemps(true);
  bool lRigiditeDependDeLaSolution(true);
  bool lRigiditeDependDesCL(true);
  if(lMsg) lMsg = lFichierChamp.reqBooleen("KDependTemps",lRigiditeDependDuTemps,false);
  if(lMsg) lMsg = lFichierChamp.reqBooleen("KDependSolution",lRigiditeDependDeLaSolution,false);
  if(lMsg) lMsg = lFichierChamp.reqBooleen("KDependCL",lRigiditeDependDesCL,false);

  // lecture des solveurs
  SolveurInstNlinPETSc* lSolveurHyd;
  SolveurLinPETSc* lSolveurMec;      
  if(lMsg){
    if(lMsg) lMsg = lFichierChamp.reqSolveur<SolveurInstNlinPETSc>("SolveurInst",lSolveurHyd);
  }
    
  if(lMsg){
    if(lTemps) lSolveurHyd->asgnTempsImposes(lTempsImposes);
    lSolveurHyd->asgnParametresMatriceMasse(lMasseDependDuTemps,lMasseDependDeLaSolution);
    lSolveurHyd->asgnParametresMatriceRigidite(lRigiditeDependDuTemps,lRigiditeDependDeLaSolution,lRigiditeDependDesCL);
    
    if(lMsg) lMsg = lFichierChamp.reqSolveur<SolveurLinPETSc>("SolveurMec",lSolveurMec);
  }

  // *****FIN de la requête des données du fichier .champs*****

  //
  // On sépare les données sur types de fichiers: hyd (humidité) et mec (déplacements)
  ChaineCar lNomHyd = lNom + std::string("_hyd");
  ChaineCar lNomMec = lNom + 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.
  TFDiffusion               lTFDiffusion;
  TFMasseGenScal            lTFMasse;
  TFCouplageAvecFctsConnues lTFCouplageHyd;
  TFContraintesInitiales    lTFContraintesInitiales;
  TFSigmaUEpsilonV3D        lTFRigidite;

  //On assigne les champs aux termes de formulation.
  if(lMsg) lMsg = lTFDiffusion.  asgnChamps(*lChampGeo, *lM, *lM, *lK,    lSchemaIntg);
  if(lMsg) lMsg = lTFMasse.      asgnChamps(*lChampGeo, *lM, *lM, *ld_b,lSchemaIntg);
  if(lCalculMec && lMsg){
    lMsg = lTFContraintesInitiales.asgnChamps(*lChampGeo,lSigmaU0,*lU,lSchemaIntg);
//     if(lMsg) lMsg = lTFCouplageHyd.asgnChamps(*lChampGeo, *lM, *lMInit, *lU, *lE, *lBeta, lSchemaIntg);
    if(lMsg) lMsg = lTFCouplageHyd.asgnChamps(*lChampGeo, *lM, *lMPrec, *lU, *lE, *lBeta, lSchemaIntg);
    if(lMsg) lMsg = lTFRigidite.   asgnChamps(*lChampGeo, *lU, *lU, *lF, *lE, lSchemaIntg);
  }

  if(lMsg){
    //On ajoute nos termes de formulations à notre problème
    lProbHyd.ajouteTF    (lTFDiffusion);
    lProbHyd.ajouteTFInst(lTFMasse);
    if(lCalculMec){
      lProbMec.ajouteTF    (lTFCouplageHyd);
      lProbMec.ajouteTF    (lTFContraintesInitiales);
      lProbMec.ajouteTF    (lTFRigidite);
    }
  }

  //On lit et initialise nos conditions aux limites
  if(lMsg) lMsg = lProbHyd.lisEtInitialiseConditionsLimites(lNomHyd, lListeCondLimHyd);
  if(lMsg && lCalculMec) lMsg = lProbMec.lisEtInitialiseConditionsLimites(lNomMec, lListeCondLimMec);
    
  // Bloc pour l'impression
  ExportGIREF lExportGIREF;
  if(lPasImpression && lMsg){
    lExportGIREF.asgnChampGeometrique(*lChampGeo);
    lExportGIREF.asgnExporteMaillage(false);
    if (lMsg) {
      lProbHyd.ajouteExportation(lExportGIREF, lPasImpression, lNomHyd);
      if(lCalculMec) lExportGIREF.ajouteChampVect3DContinu(*lU, "U");
    }
      
    // On conserve les déplacements à tous les pas de temps calcule
//    if (lMsg && lCalculMec) {
//      lProbMec.ajouteExportation(lExportGIREF, lPasImpression, lNomMec);
//    }
  }

  // Debut des pre/post traitements

  // Points d'interet dans le domaine:
  VectDyn<Vecteur3D> lCoordonnees(3);
  lCoordonnees[0].asgnXYZ(0.28,0.28,0.0);
  lCoordonnees[1].asgnXYZ(0.28,0.28,0.00525);
  lCoordonnees[2].asgnXYZ(0.28,0.28,0.0105);

  // pour M
  PPEvalueChampXYZ<ChampScalContinu> lEvalueM;
  VectDyn<DReel> lMAuxPointsDonnes(3,0.);
  const DReel& lM0 = lMAuxPointsDonnes[0];
  const DReel& lM1 = lMAuxPointsDonnes[1];
  const DReel& lM2 = lMAuxPointsDonnes[2];

  //pour U
  Vecteur3D lZero3D(0.,0.,0.);
  PPEvalueChampXYZ<ChampVect3DContinu> lEvalueU;
  VectDyn<Vecteur3D> lUAuxPointsDonnes(3,lZero3D);
  const DReel& lUz0 = lUAuxPointsDonnes[0].reqZ();
  const DReel& lUz1 = lUAuxPointsDonnes[1].reqZ();
  const DReel& lUz2 = lUAuxPointsDonnes[2].reqZ();

  // Pour les deformations
  TenseurO2Sym lZeroTO2(0.,0.,0.,0.,0.,0);
  PPEvalueChampXYZ<ChampTensO2SymLin> lEvalueS;
  VectDyn<TenseurO2Sym> lSAuxPointsDonnes(3,lZeroTO2);
  const DReel& lSigmaX0 = lSAuxPointsDonnes[0].reqComposante00();
  const DReel& lSigmaX1 = lSAuxPointsDonnes[1].reqComposante00();
  const DReel& lSigmaX2 = lSAuxPointsDonnes[2].reqComposante00();

  // pour l'ecriture des evaluations ponctuelles:
  PPAfficheValeur<DReel> lAfficheM0,lAfficheM1,lAfficheM2;
  ofstream lFichierM0((lNom + ".M0").c_str());
  ofstream lFichierM1((lNom + ".M1").c_str());
  ofstream lFichierM2((lNom + ".M2").c_str());

  PPAfficheValeur<DReel> lAfficheUz0,lAfficheUz1,lAfficheUz2;
  ofstream lFichierUz0((lNom + ".Uz0").c_str());
  ofstream lFichierUz1((lNom + ".Uz1").c_str());
  ofstream lFichierUz2((lNom + ".Uz2").c_str());

  PPAfficheValeur<DReel> lAfficheSX0,lAfficheSX1,lAfficheSX2;
  ofstream lFichierSx0((lNom + ".Sx0").c_str());
  ofstream lFichierSx1((lNom + ".Sx1").c_str());
  ofstream lFichierSx2((lNom + ".Sx2").c_str());

  // Le probleme mecanique est resolu en post-traitement du probleme hydrique:
  PPResolutionProbleme lPPResolution;

  // calcul de dM/dt
  PPDiffFinieChamp<ChampScalLin> lDeriveeML;
  PPDiffFinieChamp<ChampScalQuad> lDeriveeMQ;
  if(lMtQuad && lMt){
    lDeriveeMQ.asgnChampGeo(*lChampGeo);
    ChampScalQuad* lMSQ = (ChampScalQuad*) lM;
    ChampScalQuad* lMtSQ = (ChampScalQuad*) lMt;
    lDeriveeMQ.asgnParametres(*lMSQ,*lMtSQ);
  }
  if(lMtLin && lMt){
    lDeriveeML.asgnChampGeo(*lChampGeo);
    ChampScalLin* lMSL = (ChampScalLin*) lM;
    ChampScalLin* lMtSL = (ChampScalLin*) lMt;
    lDeriveeML.asgnParametres(*lMSL,*lMtSL);
  }

  // calcul de epsilon(u)
  PPContraintesEL3D lPPContraintesEL3D(lMail);
  ChampTensO2SymLin lEpsilon(*lChampGeo);
  lPPContraintesEL3D.asgnParametres(*lChampGeo,*lU,*lE,&lSchemaIntg);
  if(lMsg) lMsg = lPPContraintesEL3D.calculChampDeformation(lEpsilon);

  // calcul de la contrainte initiale pour le pas de temps suivant
  PPReinterpole<ChampScalContinu> lReinterpoleM;
  lReinterpoleM.asgnParametres(*lM,*lMPrec);
//   PPReinterpole<ChampScalContinu> lReinterpoleMt;
//   lReinterpoleMt.asgnParametres(*lMt,*lChampGeo,*lMtPrec,*lChampGeo);
  PPReinterpole<ChampVect3DContinu> lReinterpoleU;
  lReinterpoleU.asgnParametres(*lU,*lUPrec);


  lEvalueM.asgnParametres(*lChampGeo,*lM,lCoordonnees,lMAuxPointsDonnes);
  lProbHyd.ajoutePostTraitement(lEvalueM);
  lAfficheM0.afficheLeTemps(lProbHyd);
  lAfficheM1.afficheLeTemps(lProbHyd);
  lAfficheM2.afficheLeTemps(lProbHyd);
  lAfficheM0.asgnParametres("M",lM0,lFichierM0);
  lAfficheM1.asgnParametres("M",lM1,lFichierM1);
  lAfficheM2.asgnParametres("M",lM2,lFichierM2);
  lProbHyd.ajoutePostTraitement(lAfficheM0);
  lProbHyd.ajoutePostTraitement(lAfficheM1);
  lProbHyd.ajoutePostTraitement(lAfficheM2);
  
  // on ajoute les pre/post traitement ATTENTION la chronologie est importante!

  if(lMt && lMsg && lMtLin && lCalculMec){
    // On s'assure d'avoir calculer la derivee de M a chaque iteration non-lineaire:
    lProbHyd.ajoutePostTraitementIteration(lDeriveeML);
  }

  if(lMt && lMsg && lMtQuad && lCalculMec){
    // On s'assure d'avoir calculer la derivee de M a chaque iteration non-lineaire:
    lProbHyd.ajoutePostTraitementIteration(lDeriveeMQ);
  }

  if(lCalculMec){
    // on calcule U apres M:
    lPPResolution.asgnParametres(lProbMec,*lSolveurMec);
    lProbHyd.ajoutePostTraitement(lPPResolution,lPasDeCalcul);

    // on sauvegarde la valeur précédente de Sigma
    lProbMec.ajoutePostTraitement(lReinterpoleM);
    lProbMec.ajoutePostTraitement(lReinterpoleU);
//     lProbMec.ajoutePostTraitement(lReinterpoleMt);

    // on evalue les depacements ponctuels:
    lEvalueU.asgnParametres(*lChampGeo,*lU,lCoordonnees,lUAuxPointsDonnes);
    lProbMec.ajoutePostTraitement(lEvalueU);

    // on affiche les valeurs évaluées
    lAfficheUz0.afficheLeTemps(lProbHyd);
    lAfficheUz1.afficheLeTemps(lProbHyd);
    lAfficheUz2.afficheLeTemps(lProbHyd);
    lAfficheUz0.asgnParametres("Deplacement vertical",lUz0,lFichierUz0);
    lAfficheUz1.asgnParametres("Deplacement vertical",lUz1,lFichierUz1);
    lAfficheUz2.asgnParametres("Deplacement vertical",lUz2,lFichierUz2);
    lProbMec.ajoutePostTraitement(lAfficheUz0);
    lProbMec.ajoutePostTraitement(lAfficheUz1);
    lProbMec.ajoutePostTraitement(lAfficheUz2);
  
    // post traitement stress:
    // calcul du tenseur des contraintes
    if(lMsg){
      lProbMec.ajoutePostTraitement(lPPContraintesEL3D);

      // evaluations ponctuelles:
      lEvalueS.asgnParametres(*lChampGeo,lEpsilon,lCoordonnees,lSAuxPointsDonnes);
      lProbMec.ajoutePostTraitement(lEvalueS);

      // sauvegarde des valeurs ponctuelles:
      lAfficheSX0.afficheLeTemps(lProbHyd);
      lAfficheSX1.afficheLeTemps(lProbHyd);
      lAfficheSX2.afficheLeTemps(lProbHyd);
      lAfficheSX0.asgnParametres("Epsilon x",lSigmaX0,lFichierSx0);
      lAfficheSX1.asgnParametres("Epsilon x",lSigmaX1,lFichierSx1);
      lAfficheSX2.asgnParametres("Epsilon x",lSigmaX2,lFichierSx2);
      lProbMec.ajoutePostTraitement(lAfficheSX0);
      lProbMec.ajoutePostTraitement(lAfficheSX1);
      lProbMec.ajoutePostTraitement(lAfficheSX2);
    }

    // ecriture du stress dans le fichier solution
    if(lPasImpression) lExportGIREF.ajouteChampTensO2SymContinu(lEpsilon, "epsilon");
  }
  //
  // Résolution du problème en M et en U
  //

  if (lMsg) {
    // On résoud
    lMsg = lSolveurHyd->resoudre(lProbHyd);
  }
  lFichierM0.close();
  lFichierM1.close();
  lFichierM2.close();
  lFichierUz0.close();
  lFichierUz1.close();
  lFichierUz2.close();
  lFichierSx0.close();
  lFichierSx1.close();
  lFichierSx2.close();

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

  return 0;
}