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

// ************************************************************************
//
//  Resolution du probleme viscoelastique en utilisant le model de Maxwell generalise
//  en formulation differentielles
//
//  USAGE:<pre> testViscoElasMaxwellGeneraliseFormulationThetaSchemasHPP.[dev,opt]  préfixe_des_fichiers_d_entrée
//
//  Exemple:
//
//
// ********************************************************************

#include "MEFPPUtil.h"
#include "ChaineCar.h"
#include "ChampQijklCelluleViscoMaxwellHPP.h"
#include "ChampTijklCelluleViscoMaxwellHPP.h"
#include "ChampContraintesHPPFormulationDifferentielleMateriauVieillissant.h"
#include "ChampCelluleViscoElastiqueFormulationThetaSchemasMaxwellHPP.h"
#include "ChampPseudoRelaxationMateriauVieillisantFormulationThetaSchemas.h"
#include "ChampGeometrique.h"
#include "ChampPrecalculeDiscPtsIntegration.h"
#include "ChampScalaire.h"
#include "ChampScalContinu.h"
#include "ChampTensO2SymContinu.h"
#include "ChampTensO2SymLinElem.h"
#include "ChampTensorielO4Sym.h"
#include "ChampTensorielO4SymMin.h"
#include "ChampVect3DContinu.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 "PETScInitialisation.h"
#include "PPAfficheValeur.h"
#include "PPEvalueChampXYZ.h"
#include "PPReinterpole.h"
#include "PPNormeL2.h"
#include "ProblemeEF.h"
#include "PtIntgElem1D.h"
#include "PtIntgHexaedre.h"
#include "PtIntgPrismeTri.h"
#include "PtIntgQuadrangle.h"
#include "PtIntgTetraedre.h"
#include "PtIntgTriangle.h"
#include "SchemaIntg.h"
#include "SolveurLinPETSc.h"
#include "SolveurQuasiStatiqueNlinPETSc.h"
#include "GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin.h"
#include "GeorgesVisco_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin.h"
#include "TFSourceFVVect3D.h"
#include "traducteurERMsg.h"
#include "VerifieElementsNormalises.h"
#include <fstream>
#include <iomanip>

using namespace std;

struct ProprietePurementElastique {

  ChampTensorielO4Sym* Elasticite;
  bool                 MateriauVieillissant;

  ChampScalaire* IndicateurEtat;

  ChampTensorielO2Sym* ContraintesCourantes;
  ChampTensorielO2Sym* ContraintesPrecedentes;

};

struct ProprieteViscoElastique {

  ChampTensorielO4Sym*    Elasticite;
  ChampTensorielO4Sym*    FonctionElasticiteViscosite;
  bool                    MateriauVieillissant;
  ChampScalaire*          Multiplicateur_b;
  ChampScalaire*          Multiplicateur_bPrecedent;
  ChampScalaire*          DeriveeMultiplicateur_b;
  ChampScalaire*          DeriveeMultiplicateur_bPrecedent;
  ChampScalaire*          MultiplicateurLambda_l;
  ChampScalaire*          MultiplicateurLambda_lPrecedent;
  ChampScalaire*          IndicateurEtat;

  ChampTensorielO4Sym*    TenseurS;
  ChampTensorielO4Sym*    TenseurT;
  ChampTensorielO4Sym*    TenseurQ;
  ChampTensorielO4SymMin* PseudoContraintesCourantes;
  ChampTensorielO4SymMin* PseudoContraintesPrecedentes;

};

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

  ERMsg lMsg;

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

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

  if(lMsg){
    switch (lDimensionDuPTTypeChamp1){

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

  return lMsg;
}

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 lNomFichierChamp, lNomMaillage, lNomCL,lNomSortie;


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

  //
  // On peut avoir le meme prefixe pour le maillage, les CL et le fichierChamp
  //
  lNomMaillage = lNomFichierChamp;
  lNomCL       = lNomFichierChamp;
  // On fera la sauvegarde dans le repertoire courant avec le prefixe du fichier champ
  lNomSortie   = enleveCheminAccesFichier(lNomFichierChamp);

  //
  // Pre-traitement des donnees du fichier champs. On pourra ici lire des donnees
  // qui seront utilisé AVANT la declaration du probleme,
  // exemple: le nom du fichier maillage, le nom de l'entite geoemtrique pour un contact, etc
  //

  if(lMsg) {
    const EntierN lTailleBuffer = 65536;
    char *lLigne = new char[lTailleBuffer];
    std::ifstream lFichierChamp((lNomFichierChamp + std::string(".champs")).c_str());
    bool lNomMaillageTrouve(false);
    bool lNomCLTrouve(false);

    //
    // On cherche le nom du fichier maillage
    //
    if(lFichierChamp) {
      lFichierChamp.getline(lLigne,lTailleBuffer);

      while(lFichierChamp && lMsg && (!lNomMaillageTrouve || !lNomCLTrouve)){
	if(retournePremierMotSansEspace(lLigne) == "chainecar"){
	  ChaineCar lNomBidon;
	  const char* lPositionDansLigne = &lLigne[0];
	  lMsg = lisStringDansPtrChar(lPositionDansLigne, lNomBidon);
	  if(lMsg) lMsg = lisStringDansPtrChar(lPositionDansLigne, lNomBidon);
	  if(!lNomMaillageTrouve && lNomBidon == "_GIREF_NOM_DU_MAILLAGE") {
            //
            // On a trouve le nom du fichier maillage
            //
	    if(lMsg) {
	      lMsg = lisStringDansPtrChar(lPositionDansLigne, lNomMaillage);
	    }
	    lNomMaillageTrouve = true;
	  }

	  if( !lNomCLTrouve && lNomBidon == "_GIREF_NOM_FICHIER_CL") {
            //
            // On a trouve le nom du fichier des CL
            //
	    if(lMsg) {
	      lMsg = lisStringDansPtrChar(lPositionDansLigne, lNomCL);
	    }
	    lNomCLTrouve = true;
	  }
	}
	lFichierChamp.getline(lLigne,lTailleBuffer);
      }
      lFichierChamp.close();
    } else {
      lMsg.ajoute(ERMsg(ERMsg::ERREUR, "Erreur lors de l'ouverture du fichier: " + lNomFichierChamp + std::string(".champs")));
    }
  }

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

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

  ChampVect3DContinu    *lU          = 0;  // les déplacements
  ChampVect3DContinu    *lUPrec      = 0;  // les déplacements
  ChampVectoriel3D      *lF          = 0;  // les charges volumiques

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

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

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

  // Variable lues dans le fichier:

  Entier lNombreBrancheViscoElastique=0; // le nombre de composantes/cellules/branches viscoélastique du problème
  Entier lFrequenceImpression(1);        // frequence d'impression des résultats
  Entier lNbTempsImpose;                 // Nombre de points du temps ou on force le calcul

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

  SolveurQuasiStatiqueNlinPETSc* lSolveurQuasiStatique = 0;    // le solveur quasi statique
  SolveurLinPETSc *lSolveurLin = 0;                           // le solveur lineaire utilisé pour le point fixe du solveur QS
  DReel           *lPointeurVersLePasDeTempsDuSolveur = 0;  // Un Pointeur Vers Le Pas De Temps Du Solveur

  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
  VectDyn<DReel> lCoordFleche;           // les coordonnées dont on fait une "trace" pendant l'execution pour le calcul de la norme de la fleche

  bool lVrai(true);

  Entier lPoints1D(4);                   // le nombre de points d'integration pour une arete

  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 déclare les problèmes "vides"
  ProblemeEF lProbMec;

  if (lMsg) {
    lProbMec     .asgnNom("Probleme viscoelastique");
    lProbMec     .asgnGroupeProcessus(lGroupeGlobal);
  }

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

  VerifieElementsNormalises lVerifieElementsNormalises;
  if(lMsg && !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;


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

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

  // On assigne le champ de transformation géométrique au problème.
  if (lMsg) {
    lProbMec.asgnChampGeometrique(*lChampGeo);
  }

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

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

  // On demande  un champ qui se nomme "U" pour le champ des deplacements.
  if(lMsg) lMsg = lFichierChamp.reqChamp<ChampVect3DContinu>("U",lU);

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

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


  // On demande  un nombre reel  qui se nomme "Theta" pour le parametre du theta schema
  DReel lTheta =0;
  if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("Theta",lTheta);


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

  ProprietePurementElastique                              lProprietePurementElastique  = {0};

  ChampTensorielO4Sym*                                    lERef_Inf                    = 0;
  ChampTensorielO4Sym*                                    lE_Inf                       = 0;
  bool                                                    lMateriauVieillissant        = false;
  ChampScalaire*                                          lIndicateurEtat              = 0;
  ChampContraintesHPPFormulationDifferentielleMateriauVieillissant* lChampSigmaCourant = 0;
  ChampPrecalculeDiscPtsIntegration<ChampTensorielO2Sym>* lChampSigmaPrecedent         = 0;

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

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

  if(lMsg) lMsg =lFichierChamp.reqChamp<ChampTensorielO4Sym>("E_Infijkl", lE_Inf,false);

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

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

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

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

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

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

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

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

    // On assigne les parametres aux differents champs tensoriels

    if (lMsg) {
      lChampSigmaPrecedent                    ->asgnChamp(*lChampSigmaCourant);
      lChampSigmaCourant                      ->asgnChamps(*lU,
                                                            *lUPrec,
                                                            *lERef_Inf,
                                                            *lChampSigmaPrecedent,
                                                            *lCoef_b,
                                                            *lIndicateurEtat);
    }

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

  if(lMsg) {
    lProprietePurementElastique.Elasticite             = lERef_Inf;
    lProprietePurementElastique.MateriauVieillissant   = lMateriauVieillissant;
    lProprietePurementElastique.IndicateurEtat         = lIndicateurEtat;

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


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


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

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

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

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

    ChaineCar lNomEijkl                = lNomBranche_i + "::ERef_ijkl";
    ChaineCar lNomFonctElasViscoijkl   = lNomBranche_i + "::FonctElasViscoRef_ijkl";
    ChaineCar lNomVieillissant         = lNomBranche_i + "::MateriauVieillissant";
    ChaineCar lNomb                    = lNomBranche_i + "::Multiplicateur_b";
    ChaineCar lNombPrecedent           = lNomBranche_i + "::Multiplicateur_bPrecedent";
    ChaineCar lNombDerivee             = lNomBranche_i + "::DeriveeMultiplicateur_b";
    ChaineCar lNombDeriveePrecedent    = lNomBranche_i + "::DeriveeMultiplicateur_bPrecedent";
    ChaineCar lNombLambda              = lNomBranche_i + "::MultiplicateurLambda_l";
    ChaineCar lNombLambdaPrecedent     = lNomBranche_i + "::MultiplicateurLambda_lPrecedent";
    ChaineCar lNomEtat                 = lNomBranche_i + "::IndicateurEtat";

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

    //On demande un champ pour le tenseur des coefficients
    ChampTensorielO4Sym* lFonctionElasticiteViscosite=0;
    if(lMsg) lMsg =lFichierChamp.reqChamp<ChampTensorielO4Sym>(lNomFonctElasViscoijkl,lFonctionElasticiteViscosite);


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

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

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

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

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

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

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

    if(lMsg && lMateriauVieillissant) lMsg = ajouteListeReinterpole(*lDeriveeMultiplicateur_b,*lDeriveeMultiplicateur_bPrecedent,lVectChampAReinterpoler);

    if(lMsg && lMateriauVieillissant) lMsg = ajouteListeReinterpole(*lMultiplicateurLambda_l,*lMultiplicateurLambda_lPrecedent,lVectChampAReinterpoler);

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

    //
    // On declare les tenseurs qui sont vraiement utilisés dans le modèle
    //
    ChampTijklCelluleViscoMaxwellHPP * lChampTenseurT       = 0;
    ChampQijklCelluleViscoMaxwellHPP * lChampTenseurQ       = 0;
    ChampCelluleViscoElastiqueFormulationThetaSchemasMaxwellHPP    *  lChampTenseurS      = 0;
    ChampPseudoRelaxationMateriauVieillisantFormulationThetaSchemas*  lChampSigmaO4Courant = 0;

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

    if(lMsg){
      //
      //On construit les tenseurs qui sont vraiement utilisés dans le modèle
      //
      lChampTenseurT       = new ChampTijklCelluleViscoMaxwellHPP;
      lChampTenseurQ       = new ChampQijklCelluleViscoMaxwellHPP;
      lChampTenseurS       = new ChampCelluleViscoElastiqueFormulationThetaSchemasMaxwellHPP;
      lChampSigmaO4Courant = new ChampPseudoRelaxationMateriauVieillisantFormulationThetaSchemas;

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

      lVectChampADeleter.push_back(lChampTenseurT);
      lVectChampADeleter.push_back(lChampTenseurQ);
      lVectChampADeleter.push_back(lChampTenseurS);
      lVectChampADeleter.push_back(lChampSigmaO4Courant);
      lVectChampADeleter.push_back(lChampSigmaO4Precedent);
    }

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

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

    // On assigne les parametres aux differents champs tensoriels

    if (lMsg) {

      if(lMateriauVieillissant) {

        lChampSigmaO4Precedent                  ->asgnChamp(*lChampSigmaO4Courant);

        lChampTenseurT ->asgnChampsAvecVieillissement(*lElasticite,
                                    *lFonctionElasticiteViscosite,
                                    *lMultiplicateur_b,
                                    *lDeriveeMultiplicateur_b,
                                    *lMultiplicateur_bPrecedent,
                                    *lDeriveeMultiplicateur_bPrecedent,
                                    *lMultiplicateurLambda_l,
                                    *lMultiplicateurLambda_lPrecedent,
                                    *lIndicateurEtat,
                                    *lPointeurVersLePasDeTempsDuSolveur,
                                     lTheta);

        lChampTenseurQ ->asgnChampsAvecVieillissement(*lElasticite,
                                                      *lFonctionElasticiteViscosite,
                                                      *lMultiplicateur_b,
                                                      *lDeriveeMultiplicateur_b,
                                                      *lMultiplicateur_bPrecedent,
                                                      *lDeriveeMultiplicateur_bPrecedent,
                                                      *lMultiplicateurLambda_l,
                                                      *lMultiplicateurLambda_lPrecedent,
                                                      *lIndicateurEtat,
                                                      *lPointeurVersLePasDeTempsDuSolveur,
                                                      lTheta);

        lChampTenseurS->asgnChampsAvecVieillissement(*lElasticite,
                                                     *lFonctionElasticiteViscosite,
                                                     *lMultiplicateur_b,
                                                     *lDeriveeMultiplicateur_b,
                                                     *lMultiplicateurLambda_l,
                                                     *lIndicateurEtat,
                                                     *lPointeurVersLePasDeTempsDuSolveur,
                                                     lTheta);


        lChampSigmaO4Courant  ->asgnChampsAvecVieillissement(*lU,
                                                             *lUPrec,
                                                             *lFonctionElasticiteViscosite,
                                                             *lChampSigmaO4Precedent,
                                                             *lMultiplicateur_b,
                                                             *lDeriveeMultiplicateur_b,
                                                             *lMultiplicateurLambda_l,
                                                             *lMultiplicateurLambda_lPrecedent,
                                                             *lMultiplicateur_bPrecedent,
                                                             *lDeriveeMultiplicateur_bPrecedent,
                                                             *lIndicateurEtat,
                                                             *lPointeurVersLePasDeTempsDuSolveur,
                                                             lTheta);

      }
      if(!lMateriauVieillissant) {

         lChampSigmaO4Precedent                  ->asgnChamp(*lChampSigmaO4Courant);

        lChampTenseurT ->asgnChampsSansVieillissement(*lElasticite,
                                                      *lFonctionElasticiteViscosite,
                                                      *lPointeurVersLePasDeTempsDuSolveur,
                                                      lTheta);

        lChampTenseurQ ->asgnChampsSansVieillissement(*lElasticite,
                                                      *lFonctionElasticiteViscosite,
                                                      *lPointeurVersLePasDeTempsDuSolveur,
                                                      lTheta);

        lChampTenseurS->asgnChampsSansVieillissement(*lElasticite,
                                                     *lFonctionElasticiteViscosite,
                                                     *lPointeurVersLePasDeTempsDuSolveur,
                                                     lTheta);


        lChampSigmaO4Courant  ->asgnChampsSansVieillissement(*lU,
                                                            *lUPrec,
                                                             *lFonctionElasticiteViscosite,
                                                             *lChampSigmaO4Precedent,
                                                             *lPointeurVersLePasDeTempsDuSolveur,
                                                             lTheta);
      }

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

    }

    if(lMsg) {
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].Elasticite                       = lElasticite;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].FonctionElasticiteViscosite      = lFonctionElasticiteViscosite;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].MateriauVieillissant             = lMateriauVieillissant;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].Multiplicateur_b                 = lMultiplicateur_b;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].Multiplicateur_bPrecedent        = lMultiplicateur_bPrecedent;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].DeriveeMultiplicateur_b          = lDeriveeMultiplicateur_b;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].DeriveeMultiplicateur_bPrecedent = lDeriveeMultiplicateur_bPrecedent;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].MultiplicateurLambda_l           = lMultiplicateurLambda_l;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].MultiplicateurLambda_lPrecedent  = lMultiplicateurLambda_lPrecedent;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].IndicateurEtat                   = lIndicateurEtat;

      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].TenseurS                         = lChampTenseurS;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].TenseurT                         = lChampTenseurT;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].TenseurQ                         = lChampTenseurQ;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].PseudoContraintesCourantes   = lChampSigmaO4Courant;
      lVectProprieteBrancheViscoElastique[lNumeroDeBranche].PseudoContraintesPrecedentes = lChampSigmaO4Precedent;
    }

  }

  //
  // 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é
  //
  // 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 "lNomFichierChamp".temps
  // (un temps par ligne avec le nombre de temps impose sur la premiere ligne)


  bool lTemps(false);

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

  if(lTemps && lMsg){

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

    if(lMsg && lNbTempsImpose != 0){
      // 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 du champ U:

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

  if(lMsg && lNombreEval>0) {
    std::cout<<lNombreEval<<" évaluations des variables d'état aux points: "<<std::endl;
    lCoordonnees.asgnLongueur(lNombreEval);
    lCoordFleche.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;
    }
  }

  //
  // Est-ce qu'on a impose un nom pour les sorties dans le fichier champ
  //
  if(lMsg) lMsg = lFichierChamp.reqChaineCar("NomSortie",lNomSortie,false);

  // *********************************************
  // * BLOC: Déclaration des termes de formulation
  // *       et conditions aux limites
  // *********************************************

  // On crée l'objet qui contiendra la liste des conditions limites
  ListeConditionsLimites lListeCondLim;

  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"
    lListeCondLim.asgnGestionFichierChamps(lFichierChamp);
    lListeCondLim.asgnChampGeoEtListeEntitesGeo(*lChampGeo, lListeEntite);
    // On assigne le schéma d'intégration aux conditions aux limites
    lListeCondLim.asgnSchemaIntg(lSchemaIntg);
  }

  // On déclare les termes de formulation

  GeorgesVisco_TFElastiqueViscoElastiqueMaxwellHPPEpsilonLin        lTFRigiditeComposanteElastiquePure;
  VectDyn<GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin> lVectTFRigiditeBrancheViscoElastique(lNombreBrancheViscoElastique);
  TFSourceFVVect3D                                     lTFChargeVolumique;

  if (lMsg && lERef_Inf) {
    lMsg = lTFRigiditeComposanteElastiquePure.asgnChamps(*lChampGeo,
                                                         *lU,
                                                         *lU,
                                                         *lUPrec,
                                                         *(lProprietePurementElastique.ContraintesPrecedentes),
                                                         *lE_Inf,
                                                         *(lProprietePurementElastique.IndicateurEtat),
                                                         lSchemaIntg);
    if(lMsg) lProbMec.ajouteTF(lTFRigiditeComposanteElastiquePure);
  }

  for(Entier lNumeroDeBranche=0;lNumeroDeBranche<lNombreBrancheViscoElastique && lMsg; ++lNumeroDeBranche){
    ProprieteViscoElastique& lBranche = lVectProprieteBrancheViscoElastique[lNumeroDeBranche];
    GeorgesVisco_TFCelluleViscoElastiqueMaxwellHPPEpsilonLin& lTFRigiditeBranche = lVectTFRigiditeBrancheViscoElastique[lNumeroDeBranche];
    lMsg = lTFRigiditeBranche.asgnChamps(*lChampGeo,
                                         *lU,
                                         *lU,
                                         *lUPrec,
                                         *(lBranche.PseudoContraintesPrecedentes),
                                         *(lBranche.TenseurS),
                                         *(lBranche.TenseurT),
                                         *(lBranche.TenseurQ),
                                         lSchemaIntg);
    if(lMsg) lProbMec.ajouteTF(lTFRigiditeBranche);

  }

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

  //On lit et initialise nos conditions aux limites
  if(lMsg){
    lMsg = lProbMec.lisEtInitialiseConditionsLimites(lNomCL, lListeCondLim);
  }

  // On demande le pas d'impression
  if(lMsg) lMsg = lFichierChamp.reqScalaireConstant("FrequenceImpression",lFrequenceImpression,false);

  ExportVU lExportVU;
  if(lMsg){
    lExportVU.asgnChampGeometrique(*lChampGeo);
    lProbMec.ajouteExportation(lExportVU, lFrequenceImpression, lNomSortie);
    lMsg = lExportVU.ajouteChamp(*lU, "U");
  }

  // Points d'interet dans le domaine:
  //pour U
  PPEvalueChampXYZ<ChampVect3DContinu> lEvalueU;
  VectDyn<Vecteur3D> lUAuxPointsDonnes(lNombreEval);

//    //pour NormeL2EcartFleche
//   PPEvalueChampXYZ<ChampScalaire> lEvalueNormeFleche;
//   VectDyn<DReel> lNormeFlecheAuxPointsDonnes(lNombreEval);

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

  if(lMsg) lFichierUz.open((lNomSortie + ".U").c_str());

  // On a donné dans le fichier champs des points ou on veut une évaluation:
  if(lNombreEval && lMsg){
    lEvalueU.asgnParametres(*lU,lCoordonnees,lUAuxPointsDonnes);
    lAfficheUz.afficheLeTemps(lProbMec,lVrai);
    lAfficheUz.asgnParametres("U",lUAuxPointsDonnes,lFichierUz);
    lFichierUz << std::setprecision(8);
    lFichierUz << std::scientific;
  }


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

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

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

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

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

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

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

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

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

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

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

    if(lUtilisePasAdapte && lMsg) {

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

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

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

      //
      // On utilise un pas de temps adapté
      //
      if(lMsg) lSolveurQuasiStatique->utilisePasDeTempsVariable(lLongueurPasDeTemps,lLongueurMaximalePasDeTemps,lTempsFinal,lToleranceAdaptation);
    }

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

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

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

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

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

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

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

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

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

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

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

  //
  if(lMsg){
    std::cout<<"  Resolution debut  "<<std::endl;
    lSolveurQuasiStatique->resoudre(lProbMec);
    std::cout<<"  Resolution sortie  "<<std::endl;
  }


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

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

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

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

  return 0;
}