// ************************************************************************
// --- Logiciel MEF++
// ---
// --- Copyright 1996-2006, 2011 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 code convient à un emploi quelconque. Celui-ci est distribué
// --- sans aucune garantie implicite ou explicite.
// ************************************************************************

// ********************************************************************
//  Nom du fichier: TFSigmaUEpsilonV3DIncrementalVieillissant.cc
//
// ********************************************************************

#include "TFSigmaUEpsilonV3DIncrementalVieillissant.h"

#include "BlocMElem.h"
#include "BlocVElem.h"
#include "ChampGeometrique.h"
#include "ChampTensorielO2Sym.h"
#include "ChampTensorielO4Sym.h"
#include "ChampTensorielO4SymMin.h"
#include "ChampVectoriel3D.h"
#include "FonctionsInterpolation.h"
#include "fonctionsTenseurs.h"
#include "InfoLocalesChamp.h"
#include "InfoLocalesInterpolation.h"
#include "InfoLocalesTransformation.h"
#include "Maillage.h"
#include "SchemaIntg.h"
#include "TenseurO2Sym.h"
#include "TenseurO4Sym.h"
#include <cmath_giref>

// ********************************************************************
// Description: Constructeur
//
//  Entrée: Aucune
//
//  Sortie: Aucune
//
// ********************************************************************
TFSigmaUEpsilonV3DIncrementalVieillissant::TFSigmaUEpsilonV3DIncrementalVieillissant()
  :TFGenerique(),
    aChampU(0),
    aChampV(0),
    aChampDeplacements(0),
    aChampE(0),
    aChampEPrecedent(0),
    aChampM(0),
    aChampMPrecedent(0),
    aChampMInitial(0),
    aChampBeta(0),
    aVectInfoTransformation(0),
    aVectInfoInterpolationNV(0),
    aVectInfoInterpolationNU(0),
    aVectInfoChampU(0),
    aVectInfoChampDeplacements(0),
    aVectInfoChampE(0),
    aVectInfoChampEPrecedent(0),
    aVectInfoChampM(0),
    aVectInfoChampMPrecedent(0),
    aVectInfoChampMInitial(0),
    aVectInfoChampBeta(0),
    aVectInfoChampBetaPrecedent(0)
{
  _GIREF_CONSTRUCTEUR_INVARIANTS;
  _GIREF_TEST_INVARIANTS;
}

// ********************************************************************
//
//  Description: Constructeur par copie
//
//  Entrée: pTF le terme de formulation à copier
//
//  Sortie: Aucune
//
// ********************************************************************
TFSigmaUEpsilonV3DIncrementalVieillissant::TFSigmaUEpsilonV3DIncrementalVieillissant(const TFSigmaUEpsilonV3DIncrementalVieillissant& pTF) : TFGenerique(pTF)

{
  _GIREF_CONSTRUCTEUR_INVARIANTS;
  _GIREF_APPEL_TEST_INVARIANTS(pTF);

  _GIREF_ASSERTION(false, "devrait jamais passer ici");

  _GIREF_TEST_INVARIANTS;
}

// ********************************************************************
//
//  Description: Destructeur
//
//  Entrée: Aucune
//
//  Sortie: Aucune
//
// ********************************************************************
TFSigmaUEpsilonV3DIncrementalVieillissant::~TFSigmaUEpsilonV3DIncrementalVieillissant()
{
  _GIREF_TEST_INVARIANTS;
  _GIREF_DESTRUCTEUR_INVARIANTS;
}

// ********************************************************************
//
//  Description: Assigne les champs pour les calculs dans ce terme de formulation
//
//  Entrée: ChampGeometrique&    pChampGeo:        Le champ géométrique pour le calcul
//          ChampVectoriel3D&    pChampU:          Le champ inconnu à calculer
//          ChampVectoriel3D&    pChampV:          Le champ de fonctions tests V utilisé dans le calcul
//          ChampTensorielO4Sym& pChampE:          Le champ fonction de la valeur courante du tenseur d'élasticité
//          ChampTensorielO4Sym& pChampEPrecedent: Le champ fonction de la valeur précédente du tenseur d'élasticité
//          SchemaIntg&          pSchemaIntg:      Le schéma d'intégration pour ce terme de formulation
//
//  Sortie: Message d'erreur si l'asgnCouplage ne fonctionne pas
//
// ********************************************************************
ERMsg TFSigmaUEpsilonV3DIncrementalVieillissant::asgnChamps(ChampGeometrique&       pChampGeo,
                                                            ChampVectoriel3D&       pChampU,
                                                            ChampVectoriel3D&       pChampV,
                                                            ChampVectoriel3D&       pChampDeplacements,
                                                            ChampTensorielO4Sym&    pChampE,
                                                            const SchemaIntg&       pSchemaIntg)
{
  _GIREF_TEST_INVARIANTS;

  ERMsg lMsg(ERMsg::OK);

  asgnChampGeometrique(pChampGeo);

  aChampU              = &pChampU;
  aChampV              = &pChampV;
  aChampDeplacements   = &pChampDeplacements;
  aChampE              = &pChampE;
  aChampEPrecedent     = 0;

  asgnSchemaIntg(pSchemaIntg);

  ajouteVectInfoTransformation(aVectInfoTransformation, pChampGeo, pSchemaIntg);

  ajouteVectInfoInterpolation (aVectInfoInterpolationNV, pChampV.reqFonctionsInterpolation(), aVectInfoTransformation, pChampGeo, pSchemaIntg);
  ajouteVectInfoInterpolation (aVectInfoInterpolationNU, pChampU.reqFonctionsInterpolation(), aVectInfoTransformation, pChampGeo, pSchemaIntg);

  ajouteVectInfoChamp         (aVectInfoChampU,     pChampU,   aVectInfoTransformation, pSchemaIntg);
  ajouteVectInfoChamp         (aVectInfoChampDeplacements,   pChampDeplacements,   aVectInfoTransformation, pSchemaIntg);
  ajouteVectInfoChamp         (aVectInfoChampE,              pChampE,              aVectInfoTransformation, pSchemaIntg);

  lMsg = asgnCouplage(pChampV,pChampU);

  //La matrice élémentaire est-elle symétrique?
  //Si oui, on ne calcule que la moitié inférieure que l'on recopie dans la moitié supérieure.
  aMatriceElementaireSymetrique = aChampV == aChampU;

  _GIREF_TEST_INVARIANTS;
  return lMsg;
}

// ********************************************************************
//
//  Description: Assigne les champs pour les calculs dans ce terme de formulation
//               dans le cas ou le tenseur d'élasticité depend du temps.
//
//  Entrée: ChampGeometrique&    pChampGeo:        Le champ géométrique pour le calcul
//          ChampVectoriel3D&    pChampU:          Le champ inconnu à calculer
//          ChampVectoriel3D&    pChampV:          Le champ de fonctions tests V utilisé dans le calcul
//          ChampTensorielO4Sym& pChampE:          Le champ fonction de la valeur courante du tenseur d'élasticité
//          ChampTensorielO4Sym& pChampEPrecedent: Le champ fonction de la valeur précédente du tenseur d'élasticité
//          SchemaIntg&          pSchemaIntg:      Le schéma d'intégration pour ce terme de formulation
//
//  Sortie: Message d'erreur si l'asgnCouplage ne fonctionne pas
//
// ********************************************************************
ERMsg TFSigmaUEpsilonV3DIncrementalVieillissant::asgnChamps(ChampGeometrique&       pChampGeo,
                                                            ChampVectoriel3D&       pChampU,
                                                            ChampVectoriel3D&       pChampV,
                                                            ChampVectoriel3D&       pChampDeplacements,
                                                            ChampTensorielO4Sym&    pChampE,
                                                            ChampTensorielO4Sym&    pChampEPrecedent,
                                                            const SchemaIntg&       pSchemaIntg)
{
  _GIREF_TEST_INVARIANTS;

  ERMsg lMsg(ERMsg::OK);

  asgnChampGeometrique(pChampGeo);

  aChampU              = &pChampU;
  aChampV              = &pChampV;
  aChampDeplacements   = &pChampDeplacements;
  aChampE              = &pChampE;
  aChampEPrecedent     = &pChampEPrecedent;

  asgnSchemaIntg(pSchemaIntg);

  ajouteVectInfoTransformation(aVectInfoTransformation, pChampGeo, pSchemaIntg);

  ajouteVectInfoInterpolation (aVectInfoInterpolationNV, pChampV.reqFonctionsInterpolation(), aVectInfoTransformation, pChampGeo, pSchemaIntg);
  ajouteVectInfoInterpolation (aVectInfoInterpolationNU, pChampU.reqFonctionsInterpolation(), aVectInfoTransformation, pChampGeo, pSchemaIntg);

  ajouteVectInfoChamp         (aVectInfoChampU,     pChampU,   aVectInfoTransformation, pSchemaIntg);
  ajouteVectInfoChamp         (aVectInfoChampDeplacements,   pChampDeplacements,   aVectInfoTransformation, pSchemaIntg);
  ajouteVectInfoChamp         (aVectInfoChampE,              pChampE,              aVectInfoTransformation, pSchemaIntg);
  ajouteVectInfoChamp         (aVectInfoChampEPrecedent,     pChampEPrecedent,     aVectInfoTransformation, pSchemaIntg);

  lMsg = asgnCouplage(pChampV,pChampU);

  //La matrice élémentaire est-elle symétrique?
  //Si oui, on ne calcule que la moitié inférieure que l'on recopie dans la moitié supérieure.
  aMatriceElementaireSymetrique = aChampV == aChampU;

  _GIREF_TEST_INVARIANTS;
  return lMsg;
}

//********************************************************************
//
//  Description: Méthode d'assignation d'un champ scalaire de couplage
//
//  Entrée: pChampM:         le champ faisant le couplage
//          pChampMPrecedent la valeur du champ au pas precedent
//          pChampMInitial   la valeur initiale du champ
//          pChampDilatation le tenseur de dilatation
//
//  Sortie: Aucune
//
//********************************************************************
void TFSigmaUEpsilonV3DIncrementalVieillissant::asgnCouplageFctConnue(ChampScalaire&       pChampM,
                                                                      ChampScalaire&       pChampMPrecedent,
                                                                      ChampScalaire&       pChampMInitial,
                                                                      ChampTensorielO2Sym& pChampBeta)
{
  _GIREF_TEST_INVARIANTS;
  _GIREF_PRECONDITION(0!= &pChampM,           "Pointeur nul!");
  _GIREF_PRECONDITION(0!= &pChampMInitial,    "Pointeur nul!");
  _GIREF_PRECONDITION(0!= &pChampBeta,        "Pointeur nul!");

  aChampM             = &pChampM;
  aChampMPrecedent    = &pChampMPrecedent;
  aChampMInitial      = &pChampMInitial;
  aChampBeta          = &pChampBeta;
  aChampBetaPrecedent = 0;

  ajouteVectInfoChamp(aVectInfoChampM,          *aChampM,          aVectInfoTransformation,  reqSchemaIntg());
  ajouteVectInfoChamp(aVectInfoChampMPrecedent, *aChampMPrecedent, aVectInfoTransformation,  reqSchemaIntg());
  ajouteVectInfoChamp(aVectInfoChampMInitial,   *aChampMInitial,   aVectInfoTransformation,  reqSchemaIntg());
  ajouteVectInfoChamp(aVectInfoChampBeta,       *aChampBeta,       aVectInfoTransformation,  reqSchemaIntg());

  _GIREF_TEST_INVARIANTS;
}

//********************************************************************
//
//  Description: Méthode d'assignation d'un champ scalaire de couplage
//               dans le cas où l'on a dependance en temps du tenseur de
//               dilatation
//
//  Entrée: pChampM:             le champ faisant le couplage
//          pChampMPrecedent     la valeur du champ au pas precedent
//          pChampMInitial       la valeur initiale du champ
//          pChampDilatation     le tenseur de dilatation
//          pChampDilatationPrec le tenseur de dilatation au pas precedent
//
//  Sortie: Aucune
//
//********************************************************************
void TFSigmaUEpsilonV3DIncrementalVieillissant::asgnCouplageFctConnue(ChampScalaire&       pChampM,
                                                                      ChampScalaire&       pChampMPrecedent,
                                                                      ChampScalaire&       pChampMInitial,
                                                                      ChampTensorielO2Sym& pChampBeta,
                                                                      ChampTensorielO2Sym& pChampBetaPrecedent)
{
  _GIREF_TEST_INVARIANTS;
  _GIREF_PRECONDITION(0!= &pChampM,        "Pointeur nul!");
  _GIREF_PRECONDITION(0!= &pChampMInitial, "Pointeur nul!");
  _GIREF_PRECONDITION(0!= &pChampBeta,     "Pointeur nul!");

  aChampM             = &pChampM;
  aChampMPrecedent    = &pChampMPrecedent;
  aChampMInitial      = &pChampMInitial;
  aChampBeta          = &pChampBeta;
  aChampBetaPrecedent = &pChampBetaPrecedent;

  ajouteVectInfoChamp(aVectInfoChampM,             *aChampM,             aVectInfoTransformation,  reqSchemaIntg());
  ajouteVectInfoChamp(aVectInfoChampMPrecedent,    *aChampMPrecedent,    aVectInfoTransformation,  reqSchemaIntg());
  ajouteVectInfoChamp(aVectInfoChampMInitial,      *aChampMInitial,      aVectInfoTransformation,  reqSchemaIntg());
  ajouteVectInfoChamp(aVectInfoChampBeta,          *aChampBeta,          aVectInfoTransformation,  reqSchemaIntg());
  ajouteVectInfoChamp(aVectInfoChampBetaPrecedent, *aChampBetaPrecedent, aVectInfoTransformation,  reqSchemaIntg());

  _GIREF_TEST_INVARIANTS;
}

// ********************************************************************
//
//  Description: Ajoute aux champs les évaluations demandées pour le calcul
//
//  Entrée: bool pAssembleMatrice: Booléen qui nous indique si on assemble la matrice
//          bool pAssembleResidu : Booléen qui nous indique si on assemble le résidu
//
//  Sortie: Aucune
//
// ********************************************************************
void TFSigmaUEpsilonV3DIncrementalVieillissant::ajouteEvaluationsAFaire(bool pAssembleMatrice,bool pAssembleResidu)
{
  _GIREF_TEST_INVARIANTS;
  _GIREF_PRECONDITION(0 != aChampU, "Les champs ne sont pas assignés");

  Entier lEvaluationsChampGeo    = EvalsChampLocal::Aucune;
  Entier lEvaluationsGradFctBase = EvalsChampLocal::Aucune;
  Entier lEvaluationsGradChamp   = EvalsChampLocal::Aucune;
  Entier lEvaluationsEvalChamp   = EvalsChampLocal::Aucune;

  lEvaluationsChampGeo    |= EvalsChampLocal::EvalueDetJxPoid;
  lEvaluationsGradFctBase |= EvalsChampLocal::EvalueDeriveesPremieresFctBase;
  lEvaluationsGradChamp   |= EvalsChampLocal::EvalueDeriveesPremieresChamp;
  lEvaluationsEvalChamp   |= EvalsChampLocal::EvalueChamp;

  aChampV         -> ajouteEvaluationsAFaire(lEvaluationsGradFctBase,    reqSchemaIntg());

  reqChampGeometrique().ajouteEvaluationsTransformationAFaire(lEvaluationsChampGeo,
                                                              reqSchemaIntg());

  // Si on assemble la matrice, les évaluations suivantes seront nécessaires.
  if (pAssembleMatrice) {
    aChampU         -> ajouteEvaluationsAFaire(lEvaluationsGradFctBase,    reqSchemaIntg());
  }

  if(pAssembleResidu){
    aChampU           ->ajouteEvaluationsAFaire(lEvaluationsGradChamp, reqSchemaIntg());
    aChampDeplacements->ajouteEvaluationsAFaire(lEvaluationsGradChamp, reqSchemaIntg());
    aChampE           ->ajouteEvaluationsAFaire(lEvaluationsEvalChamp, reqSchemaIntg());
    if(aChampEPrecedent) aChampEPrecedent  ->ajouteEvaluationsAFaire(lEvaluationsEvalChamp, reqSchemaIntg());

    if (aChampM) {
      aChampM         -> ajouteEvaluationsAFaire(lEvaluationsEvalChamp, reqSchemaIntg());
      aChampMPrecedent-> ajouteEvaluationsAFaire(lEvaluationsEvalChamp, reqSchemaIntg());
      aChampMInitial  -> ajouteEvaluationsAFaire(lEvaluationsEvalChamp, reqSchemaIntg());
      aChampBeta      -> ajouteEvaluationsAFaire(lEvaluationsEvalChamp, reqSchemaIntg());
      if(aChampBetaPrecedent) aChampBetaPrecedent -> ajouteEvaluationsAFaire(lEvaluationsEvalChamp, reqSchemaIntg());
    }
  }

  _GIREF_TEST_INVARIANTS;
}

// ********************************************************************
//
//  Description: Calcul élémentaire et assemblage pour un élément quelconque.
//
//  Entrée: BlocMElem& pBlocMElemRigidite: La matrice rigidité élémentaire
//          BlocVElem& pBlocVElemResidu: Le vecteur résidu élémentaire.
//
//  Sortie: Aucun
//
// ********************************************************************
void TFSigmaUEpsilonV3DIncrementalVieillissant::evalueIntegrale(BlocMElem& pBlocMElem,
                                                                BlocVElem& pBlocVElem)
{
  _GIREF_TEST_INVARIANTS;

  const bool lAssembleMatrice  = reqAssembleMatrice();
  const bool lAssembleResidu  = reqAssembleResidu();

  // Boucle sur les points d'intégration
  const Entier lNbPoints = aVectInfoTransformation->reqLongueur();

  for (Entier k=0; k<lNbPoints; ++k) {

    const DReel&          lDetJxPoid  = (*aVectInfoTransformation)[k].reqDetJxPoid();

    _GIREF_ASSERTION(lDetJxPoid != 0.0, "Le Poids*Jacobien de l'élément est 0.0");

    TenseurO4Sym lE = (*aVectInfoChampE)[k].reqValeurChamp();

    lE *= lDetJxPoid;

    //
    // On commence par calculer le tenseur résultant contribuant au résidu:
    //

    TenseurO2Sym lContraintesPrecedentes(0.);

    if (lAssembleResidu) {

      TenseurO4Sym lDeltaE(0);
      if(aChampEPrecedent) lDeltaE = lE  - ((*aVectInfoChampEPrecedent)[k].reqValeurChamp())*lDetJxPoid;

      DReel lMMoinsM0(0);
      DReel lMMoinsMP(0);

      if (aChampBeta) {
        lMMoinsM0 = (*aVectInfoChampMPrecedent)[k].reqValeurChamp() - (*aVectInfoChampMInitial)[k].reqValeurChamp();
        lMMoinsMP = (*aVectInfoChampM)[k].reqValeurChamp() - (*aVectInfoChampMPrecedent)[k].reqValeurChamp();
      }

      if(aChampEPrecedent){

        const InfoLocalesChamp<Vecteur3D>& lLocalU   = (*aVectInfoChampDeplacements)[k];
        const Vecteur3D&                   ldUdX     = lLocalU.reqValeurDeriveeXChamp();
        const Vecteur3D&                   ldUdY     = lLocalU.reqValeurDeriveeYChamp();
        const Vecteur3D&                   ldUdZ     = lLocalU.reqValeurDeriveeZChamp();

        DReel lEpsilonXX = ldUdX.reqX();
        DReel lEpsilonYY = ldUdY.reqY();
        DReel lEpsilonZZ = ldUdZ.reqZ();

        DReel lGammaXY = ldUdX.reqY() + ldUdY.reqX();
        DReel lGammaYZ = ldUdY.reqZ() + ldUdZ.reqY();
        DReel lGammaXZ = ldUdX.reqZ() + ldUdZ.reqX();

        // e(U)
        TenseurO2Sym lTenseurResidu(lEpsilonXX,lEpsilonYY,lEpsilonZZ,0.5*lGammaXY,0.5*lGammaYZ,0.5*lGammaXZ);

        // e(U)-B(M-M0)
        if (aChampBeta) {
          TenseurO2Sym lBeta(0);
          lBeta = (*aVectInfoChampBeta)[k].reqValeurChamp();
          lTenseurResidu -= lBeta*lMMoinsM0;
        }

        // dE:(e(U)-B(M-M0))
        lContraintesPrecedentes = lDeltaE.contractionDoubleIndicesDeGauche(lTenseurResidu);
      }

      // dE:(e(U)-B(M-M0)) - E:dB(M-M0)
      if (aChampBeta && aChampBetaPrecedent) {
        TenseurO2Sym lDeltaBeta((*aVectInfoChampBeta)[k].reqValeurChamp());
        lDeltaBeta -= (*aVectInfoChampBetaPrecedent)[k].reqValeurChamp();
        lContraintesPrecedentes -= (*aVectInfoChampE)[k].reqValeurChamp().contractionDoubleIndicesDeGauche(lDeltaBeta*lMMoinsM0);
      }

      // dE:(e(U)-B(M-M0)) - E:dB(M-M0) - E:Beta(M-MP)
      if (aChampBeta) {
        TenseurO2Sym lBeta((*aVectInfoChampBeta)[k].reqValeurChamp());
        lContraintesPrecedentes -= (*aVectInfoChampE)[k].reqValeurChamp().contractionDoubleIndicesDeGauche(lBeta*lMMoinsMP);
      }

      // Puisque l'increment est l'inconnue et qu'on est en correction on a aussi:

      const InfoLocalesChamp<Vecteur3D>& lLocalDeltaU   = (*aVectInfoChampU)[k];
      const Vecteur3D&                   ldDeltaUdX     = lLocalDeltaU.reqValeurDeriveeXChamp();
      const Vecteur3D&                   ldDeltaUdY     = lLocalDeltaU.reqValeurDeriveeYChamp();
      const Vecteur3D&                   ldDeltaUdZ     = lLocalDeltaU.reqValeurDeriveeZChamp();

      DReel lEpsilonDuXX = ldDeltaUdX.reqX();
      DReel lEpsilonDuYY = ldDeltaUdY.reqY();
      DReel lEpsilonDuZZ = ldDeltaUdZ.reqZ();

      DReel lGammaDuXY = ldDeltaUdX.reqY() + ldDeltaUdY.reqX();
      DReel lGammaDuYZ = ldDeltaUdY.reqZ() + ldDeltaUdZ.reqY();
      DReel lGammaDuXZ = ldDeltaUdX.reqZ() + ldDeltaUdZ.reqX();

      // e(deltaU)
      TenseurO2Sym lTenseurDu(lEpsilonDuXX,lEpsilonDuYY,lEpsilonDuZZ,0.5*lGammaDuXY,0.5*lGammaDuYZ,0.5*lGammaDuXZ);

      lContraintesPrecedentes += (*aVectInfoChampE)[k].reqValeurChamp().contractionDoubleIndicesDeGauche(lTenseurDu);
    }

    //
    // Fin calcul du tenseur pour le résidu
    //

    // On récupère une référence aux informations locales nécessaire dans tout les cas
    const InfoLocalesInterpolation&    lLocalNV = (*aVectInfoInterpolationNV)[k];

    const VectDyn<DReel>& ldNVdX      = lLocalNV.reqDeriveeXFctBase();
    const VectDyn<DReel>& ldNVdY      = lLocalNV.reqDeriveeYFctBase();
    const VectDyn<DReel>& ldNVdZ      = lLocalNV.reqDeriveeZFctBase();

    const InfoLocalesInterpolation&    lLocalNU  = (*aVectInfoInterpolationNU)[k];
    const VectDyn<DReel>&              ldNUdX    = lLocalNU.reqDeriveeXFctBase();
    const VectDyn<DReel>&              ldNUdY    = lLocalNU.reqDeriveeYFctBase();
    const VectDyn<DReel>&              ldNUdZ    = lLocalNU.reqDeriveeZFctBase();

    const Entier          lNbFctBaseV = ldNVdX.reqLongueur();
    const Entier          lNbDDLsU    = ldNUdX.reqLongueur();

    // Optimisation de l'assemblage de la matrice
    const DReel lE0000 = lE.reqComposante0000();
    const DReel lE1111 = lE.reqComposante1111();
    const DReel lE2222 = lE.reqComposante2222();
    const DReel lE1000 = lE.reqComposante1000();
    const DReel lE2000 = lE.reqComposante2000();
    const DReel lE1010 = lE.reqComposante1010();
    const DReel lE2010 = lE.reqComposante2010();
    const DReel lE2020 = lE.reqComposante2020();
    const DReel lE1100 = lE.reqComposante1100();
    const DReel lE1110 = lE.reqComposante1110();
    const DReel lE2011 = lE.reqComposante2011();
    const DReel lE2100 = lE.reqComposante2100();
    const DReel lE2110 = lE.reqComposante2110();
    const DReel lE2120 = lE.reqComposante2120();
    const DReel lE2200 = lE.reqComposante2200();
    const DReel lE2210 = lE.reqComposante2210();
    const DReel lE2220 = lE.reqComposante2220();
    const DReel lE2111 = lE.reqComposante2111();
    const DReel lE2121 = lE.reqComposante2121();
    const DReel lE2211 = lE.reqComposante2211();
    const DReel lE2221 = lE.reqComposante2221();

    for (Entier i=0; i<lNbFctBaseV; ++i) {

      const DReel&   ldNVdXi =  ldNVdX[i];
      const DReel&   ldNVdYi =  ldNVdY[i];
      const DReel&   ldNVdZi =  ldNVdZ[i];

      const Entier i0 = 3*i;
      const Entier i1 = i0+1;
      const Entier i2 = i1+1;

      if(lAssembleResidu){

        aGamma_phi_i_0.asgnComposantes( ldNVdXi, 0      , 0      , 0.5*ldNVdYi, 0          , 0.5*ldNVdZi );
        aGamma_phi_i_1.asgnComposantes( 0      , ldNVdYi, 0      , 0.5*ldNVdXi, 0.5*ldNVdZi, 0           );
        aGamma_phi_i_2.asgnComposantes( 0      , 0      , ldNVdZi, 0          , 0.5*ldNVdYi, 0.5*ldNVdXi );

        pBlocVElem[i0] -= lContraintesPrecedentes.contractionDouble(aGamma_phi_i_0);
        pBlocVElem[i1] -= lContraintesPrecedentes.contractionDouble(aGamma_phi_i_1);
        pBlocVElem[i2] -= lContraintesPrecedentes.contractionDouble(aGamma_phi_i_2);
      }


      if(lAssembleMatrice) {

        const DReel lEpsilonV_i00 = lE0000*ldNVdXi + lE1000*ldNVdYi + lE2000*ldNVdZi;
        const DReel lEpsilonV_i01 = lE1100*ldNVdXi + lE1110*ldNVdYi + lE2011*ldNVdZi;
        const DReel lEpsilonV_i02 = lE2200*ldNVdXi + lE2210*ldNVdYi + lE2220*ldNVdZi;
        const DReel lEpsilonV_i03 = lE1000*ldNVdXi + lE1010*ldNVdYi + lE2010*ldNVdZi;
        const DReel lEpsilonV_i04 = lE2100*ldNVdXi + lE2110*ldNVdYi + lE2120*ldNVdZi;
        const DReel lEpsilonV_i05 = lE2000*ldNVdXi + lE2010*ldNVdYi + lE2020*ldNVdZi;

        const DReel lEpsilonV_i10 = lE1100*ldNVdYi + lE1000*ldNVdXi + lE2100*ldNVdZi;
        const DReel lEpsilonV_i11 = lE1111*ldNVdYi + lE1110*ldNVdXi + lE2111*ldNVdZi;
        const DReel lEpsilonV_i12 = lE2211*ldNVdYi + lE2210*ldNVdXi + lE2221*ldNVdZi;
        const DReel lEpsilonV_i13 = lE1110*ldNVdYi + lE1010*ldNVdXi + lE2110*ldNVdZi;
        const DReel lEpsilonV_i14 = lE2111*ldNVdYi + lE2110*ldNVdXi + lE2121*ldNVdZi;
        const DReel lEpsilonV_i15 = lE2011*ldNVdYi + lE2010*ldNVdXi + lE2120*ldNVdZi;

        const DReel lEpsilonV_i20 = lE2200*ldNVdZi + lE2100*ldNVdYi + lE2000*ldNVdXi;
        const DReel lEpsilonV_i21 = lE2211*ldNVdZi + lE2111*ldNVdYi + lE2011*ldNVdXi;
        const DReel lEpsilonV_i22 = lE2222*ldNVdZi + lE2221*ldNVdYi + lE2220*ldNVdXi;
        const DReel lEpsilonV_i23 = lE2210*ldNVdZi + lE2110*ldNVdYi + lE2010*ldNVdXi;
        const DReel lEpsilonV_i24 = lE2221*ldNVdZi + lE2121*ldNVdYi + lE2120*ldNVdXi;
        const DReel lEpsilonV_i25 = lE2220*ldNVdZi + lE2120*ldNVdYi + lE2020*ldNVdXi;

        for(Entier j=0; j<lNbDDLsU && (!aMatriceElementaireSymetrique || j <=i); ++j) {

          const DReel& ldNUdXj = ldNUdX[j];
          const DReel& ldNUdYj = ldNUdY[j];
          const DReel& ldNUdZj = ldNUdZ[j];

          const Entier j0 = 3*j;
          const Entier j1 = j0+1;
          const Entier j2 = j1 +1;

          const DReel lContribution_00 = ldNUdXj*lEpsilonV_i00 + ldNUdYj*lEpsilonV_i03 + ldNUdZj*lEpsilonV_i05;
          const DReel lContribution_01 = ldNUdYj*lEpsilonV_i01 + ldNUdXj*lEpsilonV_i03 + ldNUdZj*lEpsilonV_i04;
          const DReel lContribution_02 = ldNUdZj*lEpsilonV_i02 + ldNUdYj*lEpsilonV_i04 + ldNUdXj*lEpsilonV_i05;

          const DReel lContribution_10 = ldNUdXj*lEpsilonV_i10 + ldNUdYj*lEpsilonV_i13 + ldNUdZj*lEpsilonV_i15;
          const DReel lContribution_11 = ldNUdYj*lEpsilonV_i11 + ldNUdXj*lEpsilonV_i13 + ldNUdZj*lEpsilonV_i14;
          const DReel lContribution_12 = ldNUdZj*lEpsilonV_i12 + ldNUdYj*lEpsilonV_i14 + ldNUdXj*lEpsilonV_i15;

          const DReel lContribution_20 = ldNUdXj*lEpsilonV_i20 + ldNUdYj*lEpsilonV_i23 + ldNUdZj*lEpsilonV_i25;
          const DReel lContribution_21 = ldNUdYj*lEpsilonV_i21 + ldNUdXj*lEpsilonV_i23 + ldNUdZj*lEpsilonV_i24;
          const DReel lContribution_22 = ldNUdZj*lEpsilonV_i22 + ldNUdYj*lEpsilonV_i24 + ldNUdXj*lEpsilonV_i25;

          pBlocMElem(i0, j0) += lContribution_00;
          pBlocMElem(i0, j1) += lContribution_01;
          pBlocMElem(i0, j2) += lContribution_02;

          pBlocMElem(i1, j0) += lContribution_10;
          pBlocMElem(i1, j1) += lContribution_11;
          pBlocMElem(i1, j2) += lContribution_12;

          pBlocMElem(i2, j0) += lContribution_20;
          pBlocMElem(i2, j1) += lContribution_21;
          pBlocMElem(i2, j2) += lContribution_22;

        } // Fin de la boucle for(Entier j=0; j<lNbDDLsU; ++j)
      }  // Fin if (lAssembleMatrice)
    } // Fin de la boucle for(Entier i=0; j<lNbDDLsV; ++i)
  } // Fin de la boucle des points de Gauss

  //Si la matrice est symétrique, on fait la copie dans la moitié supérieure.
  if (lAssembleMatrice && aMatriceElementaireSymetrique) {
    for (Entier i = 0; i < pBlocMElem.reqNbLignes(); ++i) {
      for (Entier j = 0; j < i  ; ++j) {
        pBlocMElem(j,i) =  pBlocMElem(i,j);
      }
    }
  }

  _GIREF_TEST_INVARIANTS;
}

// ********************************************************************
//
//  Description: Méthode qui teste les invariants de la classe
//
//  Entrée: Aucune
//
//  Sortie: Arrêt du programme en cas d'erreur
//
// ********************************************************************
#ifdef GIREF_MODE_DEBUG
void TFSigmaUEpsilonV3DIncrementalVieillissant::testInvariants() const
{
  _GIREF_TEST_OBJET_INVARIANTS;
  TFGenerique::testInvariants();
}
#endif // #ifdef GIREF_MODE_DEBUG