// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL  : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.

//#include "Debug.h"
#include "ElemThermi.h"
#include <iomanip>
#include "ConstMath.h"
#include "Util.h"
#include "Coordonnee2.h"
#include "Coordonnee3.h"
#include "TypeConsTens.h"
#include "TypeQuelconqueParticulier.h"

        
// retourne la liste de tous les  types de ddl interne actuellement utilisés
// par l'élément (actif ou non), sont exclu de cette liste les ddl des noeuds
// reliés à l'élément (ddl implique grandeur uniquement scalaire !)
List_io <Ddl_enum_etendu> ElemThermi::Les_type_de_ddl_internes(bool) const
   { List_io <Ddl_enum_etendu> ret;
     switch (ParaGlob::Dimension())
     { case 1 :
        {ret.push_back(TEMP); ret.push_back(FLUXD1);
         ret.push_back(GRADT1); ret.push_back(DGRADT1);
//         ret.push_back(string("Green-Lagrange11"));ret.push_back(string("Almansi11"));ret.push_back(string("logarithmique11"));
         break;
         }
       case 2 :
        {ret.push_back(TEMP); ret.push_back(FLUXD1);ret.push_back(FLUXD2);
         ret.push_back(GRADT1); ret.push_back(GRADT2);
         ret.push_back(DGRADT1); ret.push_back(DGRADT2);
//         ret.push_back(EPS22); ret.push_back(EPS12);
         break;
         }
       case 3 :
        {ret.push_back(TEMP);
         ret.push_back(FLUXD1);ret.push_back(FLUXD2);ret.push_back(FLUXD3);
         ret.push_back(GRADT1); ret.push_back(GRADT2); ret.push_back(GRADT3);
         ret.push_back(DGRADT1); ret.push_back(DGRADT2);
//         ret.push_back(string("Green-Lagrange11")); ret.push_back(string("Green-Lagrange22"));
         break;
         }
      };
    // cas des invariants
//    tab_Dee(107).nom = "norme_gradT"     ;tab_Dee(107).enu = GRADT1;tab_Dee(107).posi_nom = nbenumddl + 107;
//    tab_Dee(108).nom = "norme_DgradT"     ;tab_Dee(108).enu = DGRADT1;tab_Dee(108).posi_nom = nbenumddl + 108;
//    tab_Dee(109).nom = "norme_dens_flux"     ;tab_Dee(109).enu = FLUXD1;tab_Dee(109).posi_nom = nbenumddl + 109;
    
    ret.push_back(string("norme_gradT"));
    ret.push_back(string("norme_DgradT"));
    ret.push_back(string("norme_dens_flux"));
    // cas de l'existence de l'erreur
    if (fluxErreur != NULL)  ret.push_back(ERREUR);
//    // cas des énergies locales
//    ret.push_back(string("energie_elastique"));
//    ret.push_back(string("dissipation_plastique"));
//    ret.push_back(string("dissipation_visqueuse"));
    // cas des coordonnées du points
    switch (ParaGlob::Dimension())
      {case 3: ret.push_back(X3);
       case 2: ret.push_back(X2);
       case 1: ret.push_back(X1);
       };
    // retour
    return ret;   
   };

// retourne la liste abondée de tous les  données particulières interne actuellement utilisés
// par l'élément (actif ou non), sont exclu de cette liste les données particulières des noeuds
// reliés à l'élément
List_io <TypeQuelconque> ElemThermi::Les_types_particuliers_internes(bool absolue) const
  { List_io <TypeQuelconque> ret;
    // on s'intéresse au information stockée au niveau des points d'intégration
    // 1° récupéreration les grandeurs 
    //    a) cas des grandeurs liées à la déformation
    def->ListeGrandeurs_particulieres(absolue,ret);
    //    b) cas des grandeurs liées à la loi de comportement mécanique 
    loiTP->ListeGrandeurs_particulieres(absolue,ret);
    //    c) cas des grandeurs liées à la loi de comportement thermo-physique 
    if (loiComp != NULL) loiComp->ListeGrandeurs_particulieres(absolue,ret);
    //    d) cas des infos stockés aux points d'intégration, donc par l'élément
    Grandeur_scalaire_double grand_courant; // def d'une grandeur courante 
    // def d'un type quelconque représentatif à chaque grandeur
    // a priori ces grandeurs sont défini aux points d'intégration 
    //    e) cas des infos stockés à l'élément
    //  $$$ cas du volume total : 
    TypeQuelconque typQ3(VOLUME_ELEMENT,SIG11,grand_courant);
    ret.push_back(typQ3);

    //    f) cas des infos calculable à la demande sur l'élément
    //  $$$ cas des volumes entre les plans et l'élément dans le cas ou l'élément des de surface
    //      et que le calcul est 3D :
    if ((ParaGlob::Dimension()==3)&&(Type_geom_generique(id_geom)==SURFACE))
    {Grandeur_coordonnee grandCoordonnee(3); // un type courant 
       TypeQuelconque typQ4(VOL_ELEM_AVEC_PLAN_REF,SIG11,grandCoordonnee);
       ret.push_back(typQ4);
     };
    // liberation des tenseurs intermediaires
    LibereTenseur(); 
    return ret;
    }; 
     
// retourne la liste des grandeurs évoluées interne actuellement utilisés
// par l'élément, c'est-à-dire directement sous forme de vecteur, tenseurs ....		 
List_io <TypeQuelconque> ElemThermi::Les_type_evolues_internes(bool absolue) const
  { List_io <TypeQuelconque> liTQ;
    // def des grandeurs courantes de type coordonnee
   {Coordonnee coor(ParaGlob::Dimension()); // un type coordonnee typique
    // maintenant on définit une grandeur typique de type Grandeur_coordonnee
    Grandeur_coordonnee gt(coor); 
    // def d'un type quelconque représentatif  pour chaque grandeur
    {TypeQuelconque typQ(FLUXD,FLUXD1,gt);liTQ.push_back(typQ);};
    {TypeQuelconque typQ(GRADT,GRADT1,gt);liTQ.push_back(typQ);};
    {TypeQuelconque typQ(DGRADT,DGRADT1,gt);liTQ.push_back(typQ);};
    // pour la position il s'agit de la position du point d'intégration
    {TypeQuelconque typQ(POSITION_GEOMETRIQUE,FLUXD1,gt);liTQ.push_back(typQ);};
    };
   
    // retour
    return liTQ;   
   };
        
// récupération de grandeurs particulières au numéro d'ordre  = iteg  
// celles-ci peuvent être quelconques
// en retour liTQ est modifié et contiend les infos sur les grandeurs particulières
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
void ElemThermi::Grandeur_particuliere (bool absolue,List_io<TypeQuelconque>& liTQ,int iteg)
   { // on s'intéresse au information stockée au niveau des points d'intégration
     // pour tenir compte du fait qu'il peut y avoir plusieurs fois la même information, ces dernières sont systématiquement
     // stockées dans des tableaux, et on associe à ce tableau un indice qui indique aux méthodes appelées, à quel endroit
     // elles doivent écrire les infos. les indices sont stockés dans une liste de même dim que liTQ
     list <int> decal(liTQ.size(),0.);  // on initi tous les membres à 0

     // 1° on commence par récupérer les grandeurs 
     //    a) cas des grandeurs liées à la déformation
     // met à jour les données spécifiques du point considéré
     def->Mise_a_jour_data_specif(tabSaveDefDon(iteg));		
     def->Grandeur_particuliere(absolue,liTQ);
     //    b) cas des grandeurs liées à la loi de comportement thermo-physique 
     loiComp->Grandeur_particuliere(absolue,liTQ,tabSaveDon(iteg),decal);
// pour le débug ----
//cout << "\n debug ElemThermi::Grandeur_particuliere  " <<  liTQ << endl; // débug
// fin pour le débug ----     
     //    c) cas des grandeurs liées à la loi de comportement thermo-physique 
     if (loiTP != NULL) loiTP->Grandeur_particuliere(absolue,liTQ,tabSaveTP(iteg),decal);
     // 2° on change éventuellement de repère (local en global) si nécessaire
     //    a) définition des grandeurs qui sont indépendante de la boucle qui suit
     def->ChangeNumInteg(iteg);  // on change le numéro de point d'intégration courant
     // dimension des vecteurs de base de la métrique
     int dim = met->Dim_vec_base();
     // calcul des matrices de passages
     // calcul des matrices de passages
     Mat_pleine* gamma0=NULL; Mat_pleine* gammat=NULL; Mat_pleine* gammatdt=NULL;
     Mat_pleine* beta0=NULL; Mat_pleine* betat=NULL; Mat_pleine* betatdt=NULL;
     List_io<TypeQuelconque>::iterator il,ilfin = liTQ.end();
     for (il=liTQ.begin();il!=ilfin;il++)
      { TypeQuelconque& tipParticu = (*il); // pour simplifier
        // n'intervient que si la grandeur est active: 
        if ((*il).Activite())
        { // on regarde si c'est une grandeur locale ou pas 
          if (tipParticu.TypeExpressionGrandeur()==0)
            // cas d'une grandeur exprimée dans le repère locale, changement de repère
            { // calcul des matrices de passages, au premier passage
              if (gamma0 == NULL)
               { gamma0 = new Mat_pleine(dim,dim); gammat = new Mat_pleine(dim,dim);gammatdt= new Mat_pleine(dim,dim);
                 const Met_abstraite::Info0_t_tdt& ex = def->Remont0_t_tdt(absolue,*gamma0,*gammat,*gammatdt);
                 beta0 = new Mat_pleine(gamma0->Inverse().Transpose());
                 betat = new Mat_pleine(gammat->Inverse().Transpose());
                 betatdt = new Mat_pleine(gammatdt->Inverse().Transpose());
               };
              tipParticu.Change_repere(*beta0,*betat,*betatdt,*gamma0,*gammat,*gammatdt);
            };
          // on utilise la boucle également pour mettre à jour les grandeurs particulières liés à l'élément  
            
          // 1) -----cas du module de compressibilité ou de cisaillement totale
          switch (tipParticu.EnuTypeQuelconque().EnumTQ())
            {//    e) cas des infos stockés à l'élément
             case  VOLUME_ELEMENT:
              { *((Grandeur_scalaire_double*) tipParticu.Grandeur_pointee())=volume;break;}
             //    d) cas des infos calculable à la demande sur l'élément
             case VOL_ELEM_AVEC_PLAN_REF:
              { *((Grandeur_coordonnee*) tipParticu.Grandeur_pointee())=VolumePlan();break;}
             default: break; // sinon rien
            };
        };

       };
      // suppression des bases de passages
      if (gamma0 != NULL)
       {delete gamma0; delete gammat; delete gammatdt; delete beta0; delete betat; delete betatdt;};
      // liberation des tenseurs intermediaires
      LibereTenseur();   
    };                             	
// retourne la liste de toutes les  grandeurs quelconques relatives aux faces de
// l'élément (actif ou non),
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
List_io <TypeQuelconque> ElemThermi::Les_type_quelconque_de_face(bool absolue) const
   {  cout << "\n **** méthode  non implantée **** "
           << "\n List_io <TypeQuelconque> ElemThermi::Les_type_quelconque_de_face(.... "
           << flush;
      Sortie(1);
   };
   
// retourne la liste de toutes les  grandeurs quelconques relatives aux arêtes de
// l'élément (actif ou non),
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
List_io <TypeQuelconque> ElemThermi::Les_type_quelconque_de_arete(bool absolue) const
{  cout << "\n **** méthode  non implantée **** "
        << "\n List_io <TypeQuelconque> ElemThermi::Les_type_quelconque_de_arete(.... "
        << flush;
   Sortie(1);
};
                         
// récupération de grandeurs particulières pour une face au numéro d'ordre  = iteg
// celles-ci peuvent être quelconques
// en retour liTQ est modifié et contiend les infos sur les grandeurs particulières
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
void ElemThermi::Grandeur_particuliere_face(bool absolue,List_io<TypeQuelconque>& liTQ,int face, int iteg)
{  cout << "\n **** méthode  non implantée **** "
        << "\n List_io <TypeQuelconque> ElemThermi::Grandeur_particuliere_face(.... "
        << flush;
   Sortie(1);
};

// récupération de grandeurs particulières pour une arête au numéro d'ordre  = iteg
// celles-ci peuvent être quelconques
// en retour liTQ est modifié et contiend les infos sur les grandeurs particulières
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
void ElemThermi::Grandeur_particuliere_arete(bool absolue,List_io<TypeQuelconque>& liTQ,int arete, int iteg)
{  cout << "\n **** méthode  non implantée **** "
        << "\n List_io <TypeQuelconque> ElemThermi::Grandeur_particuliere_arete(.... "
        << flush;
   Sortie(1);
};

// --- transfert des grandeurs des points d'intégration aux noeuds 
// transfert de ddl des points d'intégrations (de tous) aux noeuds d'un éléments (on ajoute aux noeuds, on ne remplace pas)
// les ddl doivent déjà exister aux noeuds sinon erreur
// il doit s'agir du même type de répartition de pt d'integ pour toutes les grandeurs
// tab_val(i)(j) : valeur associée au i ième pt d'integ et au j ième ddl_enum_etendu
void ElemThermi::TransfertAjoutAuNoeuds(const List_io < Ddl_enum_etendu >& lietendu
                                     ,const Tableau <Tableau <double> > & tab_val,int cas) 
{List_io < Ddl_enum_etendu >::const_iterator ili,ilifin = lietendu.end();
 // on récupère l'élément géométrique associé au ddl pur du premier ddl étendu,
 ili = lietendu.begin();
 Enum_ddl enu = (*ili).Enum();
 const ElemGeomC0& elemgeom = ElementGeometrie(enu);
 // on récupère le nombre de noeud de l'élément géométrique (qui peut-être différent de celui de l'élément
 // par ex: element sfe, donc on utilise les noeuds de l'élément géom. On considère que ce sont les premiers
 // du tableau des noeuds
 int nbn = elemgeom.Nbne();
 
 int indice = 1;
 for (ili = lietendu.begin();ili!=ilifin;ili++,indice++)
  { // récupération des tableaux d'indiçage pour le calcul aux noeuds
    const ElemGeomC0::ConteneurExtrapolation & contExtra = elemgeom.ExtrapolationNoeud(cas);
    // calcul des valeurs aux noeuds
    // val_au_noeud(i) = somme_(de j=indir(i)(1) à indir(i)(taille(indice(i)) )) {tab(i)(j) * val_pt_integ(j) }
	   for (int ine=1;ine<=nbn;ine++)
     { Ddl_etendu& de = tab_noeud(ine)->ModifDdl_etendu(*ili);
       Tableau <int>& indir = contExtra.indir(ine); // pour simplifier
       int nbindir = indir.Taille();
       for (int j=1; j<= nbindir; j++)
           de.Valeur() += contExtra.tab(ine)(indir(j)) * tab_val(indir(j))(indice);
     };
  };
};

// transfert de type quelconque des points d'intégrations (de tous) aux noeuds d'un éléments (on ajoute aux noeuds, 
// on ne remplace pas). Les types quelconques doivent déjà exister.
// un tableau tab_liQ  correspondent aux grandeurs quelconque pour tous les pt integ. tab_liQ(i) est associé au pt d'integ i
// Toutes les listes sont identiques au niveau des descripteurs (types...) ce sont uniquement les valeurs numériques
// c-a-dire les valeurs associées à TypeQuelconque::Grandeur qui sont différentes (elles sont associées à chaque pt d'integ)
// liQ_travail: est une liste de travail qui sera utilisée dans le transfert                        
void ElemThermi::TransfertAjoutAuNoeuds(const Tableau <List_io < TypeQuelconque > >& tab_liQ
                                     ,List_io < TypeQuelconque > & liQ_travail,int cas)	  
{int nb_pt_integ = lesPtIntegThermiInterne->NbPti();
// ---- modif transitoire pour la sortie des coques !!!!
 // on s'occupe uniquement de la première couche d'élément
 if ((this->PoutrePlaqueCoque() == COQUE)||(this->PoutrePlaqueCoque() == PLAQUE))
   nb_pt_integ = this->NbPtIntegSurface(FLUXD1); // SIG11 arbitraire car c'est les pt d'integ classique


//------- fin modif transitoire --------
 if (tab_liQ.Taille() == 0)
    return; // cas où on n'a rien à faire
 #ifdef MISE_AU_POINT
    if (nb_pt_integ != tab_liQ.Taille())
      { cout << "\n erreur de dimension, la taille de tab_liQ= " << tab_liQ.Taille() << " n'est pas egale "
             << " au nombre de pt_d'integ de l'element : " << nb_pt_integ ;
        cout << "\n ElemThermi::TransfertAjoutAuNoeuds(... " << endl;
        Sortie(1);
      }  
 #endif
//---- pour le debug
//cout << " \n debug ElemThermi::TransfertAjoutAuNoeuds " << tab_liQ << endl;
//--- fin pour le debug 
 // sinon on a au moins un élément, on récupère la 1ère liste qui nous servira pour les descripteurs
 List_io <TypeQuelconque > & liQ_ref = tab_liQ(1); // pour simplifier
 List_io < TypeQuelconque >::const_iterator ili_ref,ilifin_ref=liQ_ref.end();
 // la liste de travail
 List_io < TypeQuelconque >::iterator ili_travail = liQ_travail.begin();
 // on définit un tableau d'itérator de dimension = le nb de pt d'integ
 static Tableau < List_io < TypeQuelconque >::const_iterator > tab_iterator;
 // on le redimentionne si besoin est
 tab_iterator.Change_taille(nb_pt_integ);
 // on définit ses éléments
 for (int ia=1;ia<=nb_pt_integ;ia++)
   tab_iterator(ia) = tab_liQ(ia).begin();
 // on boucle sur les éléments de la liste de référence
 for (ili_ref = liQ_ref.begin();ili_ref != ilifin_ref;ili_ref++,ili_travail++)
   { // on récupère l'élément géométrique associé au ddl pur dont les valeurs sont calculées aux mêmes
     // endroits que la grandeur quelconque
     Enum_ddl enu = (*ili_ref). Enum();
     const ElemGeomC0& elemgeom = ElementGeometrie(enu);
     // on récupère le nombre de noeud de l'élément géométrique (qui peut-être différent de celui de l'élément
     // par ex: element sfe, donc on utilise les noeuds de l'élément géom. On considère que ce sont les premiers
     // du tableau des noeuds
     int nbn = elemgeom.Nbne();
     // récupération des tableaux d'indiçage pour le calcul aux noeuds
     const ElemGeomC0::ConteneurExtrapolation & contExtra = elemgeom.ExtrapolationNoeud(cas);
     // on boucle sur les noeuds, ceci est due à la technique utilisée pour accumuler dans les conteneurs du noeud
     for (int ine=1;ine<=nbn;ine++)
        {TypeQuelconque& tq = tab_noeud(ine)->ModifGrandeur_quelconque((*ili_ref).EnuTypeQuelconque());
         TypeQuelconque::Grandeur & grandeur_au_noeud = *(tq.Grandeur_pointee());
	     // on récupère une grandeur du même type que les grandeurs à manipuler, ceci pour accumuler
	     TypeQuelconque::Grandeur & stockage_inter = *((*ili_travail).Grandeur_pointee());
         // calcul des valeurs aux noeuds
         // val_au_noeud(i) = somme_(de j=indir(i)(1) à indir(i)(taille(indice(i)) )) {tab(i)(j) * val_pt_integ(j) }
         Tableau <int>& indir = contExtra.indir(ine); // pour simplifier
         int nbindir = indir.Taille();
         for (int j=1; j<= nbindir; j++)
       	  { stockage_inter *= 0.; // met à 0 la grandeur inter
       	    // tab_iterator(j) pointe sur la grandeur quelconque du point d'intégration voulu 
       	    stockage_inter = *((*(tab_iterator(indir(j)))).Const_Grandeur_pointee()); // récup de la grandeur pointée
       	    stockage_inter *= contExtra.tab(ine)(indir(j)); // * par le coeff
       	    grandeur_au_noeud += stockage_inter; // stockage dans le noeud	
       	  };
       };
     // on incrémente les itérators de points d'intégration
     for (int ib=1;ib<=nb_pt_integ;ib++)
         (tab_iterator(ib))++;
       	};
};
	     
 
// accumulation aux noeuds de grandeurs venant des éléments vers leurs noeuds (exemple la pression appliquée)
// autres que celles aux pti classiques, mais directements disponibles
// le contenu du conteneur stockées dans liQ est utilisé en variable intermédiaire
void ElemThermi::Accumul_aux_noeuds(const List_io < Ddl_enum_etendu >& lietendu
                                    ,List_io < TypeQuelconque > & liQ,int cas)
{// ... cas des ddl étendu ...
 {List_io < Ddl_enum_etendu >::const_iterator ili,ilifin = lietendu.end();
  ili = lietendu.begin();
//  for (ili = lietendu.begin();ili!=ilifin;ili++)
//   {// on traite en fonctions des cas particuliers
////debug
////cout << "\n debug ElemMeca::Accumul_aux_noeuds"
////     << " *ili = " << *ili << endl;
//// fin debug
//   
//    switch ((*ili).Position()-NbEnum_ddl())
//     { case 94 : // cas des pressions externes 
//        {if (lesChargeExtSurEle != NULL)
//          if (lesChargeExtSurEle->LesPressionsExternes() != NULL) // cas où il existe des pressions sauvegardées
//           { Tableau <Tableau <Pression_appliquee> >& lesPressionsExternes = *lesChargeExtSurEle->LesPressionsExternes();
//             int nb_face = lesPressionsExternes.Taille();
//             for (int n_face=1;n_face<=nb_face;n_face++) // on boucle sur les faces
//              if (lesPressionsExternes(n_face).Taille() != 0)
//              { Tableau <Pression_appliquee>& tab_press_appliquee = (lesPressionsExternes(n_face)); // pour simplifier
//                int t_tail = tab_press_appliquee.Taille();
//                if (t_tail != 0) // cas où on a une face chargée
//                  { // on récupère l'élément frontière associée (on considère qu'il y en a forcément 1)
//                    int num_face = n_face; // init par défaut
//                    if(ParaGlob::AxiSymetrie())  // dans le cas d'éléments axysymétriques, les faces
//                      {num_face += posi_tab_front_lin;}; // sont en faites les frontières lignes
//                    const ElemGeomC0 & elemgeom = tabb(num_face)->ElementGeometrique();
//                    int nbn = elemgeom.Nbne();
//                    const ElemGeomC0::ConteneurExtrapolation & contExtra = elemgeom.ExtrapolationNoeud(cas);
//                    Tableau <Noeud *>& tab_noe = tabb(num_face)->TabNoeud(); // récup du tableau de noeuds de la frontière
//                    for (int ine=1;ine<=nbn;ine++)
//                      { Ddl_etendu& de = tab_noe(ine)->ModifDdl_etendu(*ili);
//                        Tableau <int>& indir = contExtra.indir(ine); // pour simplifier
//                        int nbindir = indir.Taille();
//                        for (int j=1; j<= nbindir; j++)
//{                           de.Valeur() += contExtra.tab(ine)(indir(j)) * tab_press_appliquee(indir(j)).press;
////debug
////cout << "\n debug ElemMeca::Accumul_aux_noeuds"
////     << " de.Valeur() = " << de.Valeur() << endl;
//// fin debug
//}                      };
//                  };
//              };
//           };
//          break;
//        }
//       default :
//         break;
//        // sinon on ne fait rien, cela signifie que ce n'est pas un ddl géré par l'élément
//     };
//   };
 };

 // ... cas des grandeurs quelconque ...
 // on traite d'abord le cas où on n'a rien à faire
 if (liQ.size() == 0)
    return; 
//---- pour le debug
//cout << " \n debug ElemMeca::TransfertAjoutAuNoeuds " << tab_liQ << endl;
//--- fin pour le debug 
 {List_io < TypeQuelconque >::const_iterator ili,ilifin=liQ.end();
  
  // on boucle sur les grandeurs de la liste 
  for (ili = liQ.begin();ili != ilifin;ili++)
   { // on récupère l'énuméré
     TypeQuelconque_enum_etendu enutq = (*ili).EnuTypeQuelconque();
     switch (enutq.EnumTQ())
      { //        case VECT_REAC:
//         // là on ne fait rien car la grandeur est directement à construire via les noeuds
//         // et ceci est fait dans la routine:
//         // lesMaillages->PassageInterneDansNoeud_composantes_vers_vectorielles()
//         break;

       default:
     // VECT_PRESSION pour l'instant il n'y a aucune grandeur a priori qui est disponible directement
     // donc message d'erreur
        cout << "\n cas de grandeur actuellement non traité "
             << "\n pas de grandeur " << ((*ili).EnuTypeQuelconque().NomPlein()) << "  dans l'element "
             << "\n ou cas non implante pour l'instant"
             << "\n ElemMeca::Accumul_aux_noeuds(....";
      };
   };
  };
};

	     
// procedure permettant de completer l'element, en ce qui concerne les variables gérés
// par ElemThermi, apres sa creation avec les donnees du bloc transmis
// peut etre appeler plusieurs fois
Element* ElemThermi::Complete_ElemThermi(BlocGen & bloc,LesFonctions_nD*  lesFonctionsnD)
    { // cas de la masse volumique
      if (bloc.Nom(1) == "masse_volumique")
       { masse_volumique = bloc.Val(1);
         return this;
        }
      // cas de données thermique   
      else if (bloc.Nom(1) == "dilatation_thermique")
       { if (bloc.Val(1) == 1.) {dilatation=true;} else {dilatation = false;};
         return this;
        }
      // cas de la définition d'une intégrale sur le volume
      // on prévoit une intégrale via les pti, on utilise donc le ddl EPS11 pour la position
      else if (bloc.Nom(1) == "integrale_sur_volume_")
       { if (bloc.Nom(2) == "un_ddl_etendu_")
           { // il s'agit d'une intégrale d'un ddl étendu de nom: bloc.Nom(3)
             Ddl_etendu ddl(Ddl_enum_etendu(bloc.Nom(3))) ;// on crée un ddl ad hoc
             Grandeur_Ddl_etendu grand_courant(ddl,bloc.Nom(4));// le conteneur ad hoc
             TypeQuelconque typQ(INTEG_SUR_VOLUME,EPS11,grand_courant);
             if (integ_vol_typeQuel == NULL)
                {integ_vol_typeQuel = new Tableau <TypeQuelconque>(1);
                 (*integ_vol_typeQuel)(1) = typQ;
                 index_Integ_vol_typeQuel = new Tableau <int>(1);
                 (*index_Integ_vol_typeQuel)(1) = (int) bloc.Val(1);
                }
             else {int taille = integ_vol_typeQuel->Taille()+1;
                   integ_vol_typeQuel->Change_taille(taille);
                   (*integ_vol_typeQuel)(taille) = typQ;
                   index_Integ_vol_typeQuel->Change_taille(taille);
                   (*index_Integ_vol_typeQuel)(taille) = (int) bloc.Val(1);
                  };
           }
         else if (bloc.Nom(2) == "une_fonction_nD_")
           { // il s'agit d'une intégrale d'une fonction nD de nom: bloc.Nom(3)
             // on récupère le pointeur de fonction correspondant:
             Fonction_nD * fct = lesFonctionsnD->Trouve(bloc.Nom(3));
             Grandeur_Vecteur_Nommer grand_courant(bloc.Nom(3),1,fct);// par défaut une taille de 1
             TypeQuelconque typQ(INTEG_SUR_VOLUME,EPS11,grand_courant);
             if (integ_vol_typeQuel == NULL)
               { integ_vol_typeQuel = new Tableau <TypeQuelconque>(1);
                 (*integ_vol_typeQuel)(1) = typQ;
                 index_Integ_vol_typeQuel = new Tableau <int>(1);
                 (*index_Integ_vol_typeQuel)(1) = (int) bloc.Val(1);
               }
             else {int taille = integ_vol_typeQuel->Taille()+1;
                   integ_vol_typeQuel->Change_taille(taille);
                   (*integ_vol_typeQuel)(taille) = typQ;
                   index_Integ_vol_typeQuel->Change_taille(taille);
                   (*index_Integ_vol_typeQuel)(taille) = (int) bloc.Val(1);
                  };
           };
         return this;
       }
     
      // cas de la définition d'une intégrale sur le volume et le temps
      else if (bloc.Nom(1) == "integrale_sur_vol_et_temps_")
        { if (bloc.Nom(2) == "un_ddl_etendu_")
           { // il s'agit d'une intégrale d'un ddl étendu de nom: bloc.Nom(3)
             Ddl_etendu ddl(Ddl_enum_etendu(bloc.Nom(3))) ;// on crée un ddl ad hoc
             Grandeur_Ddl_etendu grand_courant(ddl,bloc.Nom(4));// le conteneur ad hoc
             TypeQuelconque typQ(INTEG_SUR_VOLUME_ET_TEMPS,EPS11,grand_courant);
             if (integ_vol_t_typeQuel == NULL)
                {integ_vol_t_typeQuel = new Tableau <TypeQuelconque>(1);
                 (*integ_vol_t_typeQuel)(1) = typQ;
                 index_Integ_vol_t_typeQuel = new Tableau <int>(1);
                 (*index_Integ_vol_t_typeQuel)(1) = (int) bloc.Val(1);
                }
             else {int taille = integ_vol_t_typeQuel->Taille()+1;
                   integ_vol_t_typeQuel->Change_taille(taille);
                   (*integ_vol_t_typeQuel)(taille) = typQ;
                   index_Integ_vol_t_typeQuel->Change_taille(taille);
                   (*index_Integ_vol_t_typeQuel)(taille) = (int) bloc.Val(1);
                  };
           }
         else if (bloc.Nom(2) == "une_fonction_nD_")
           { // il s'agit d'une intégrale d'une fonction nD de nom: bloc.Nom(3)
             // on récupère le pointeur de fonction correspondant:
             Fonction_nD * fct = lesFonctionsnD->Trouve(bloc.Nom(3));
             Grandeur_Vecteur_Nommer grand_courant(bloc.Nom(3),1,fct);// par défaut une taille de 1
             TypeQuelconque typQ(INTEG_SUR_VOLUME_ET_TEMPS,EPS11,grand_courant);
             if (integ_vol_t_typeQuel == NULL)
                {integ_vol_t_typeQuel = new Tableau <TypeQuelconque>(1);
                 (*integ_vol_t_typeQuel)(1) = typQ;
                 index_Integ_vol_t_typeQuel = new Tableau <int>(1);
                 (*index_Integ_vol_t_typeQuel)(1) = (int) bloc.Val(1);
                }
             else {int taille = integ_vol_t_typeQuel->Taille()+1;
                   integ_vol_t_typeQuel->Change_taille(taille);
                   (*integ_vol_t_typeQuel)(taille) = typQ;
                   index_Integ_vol_t_typeQuel->Change_taille(taille);
                   (*index_Integ_vol_t_typeQuel)(taille) = (int) bloc.Val(1);
                  };
           };
         return this;
       }
      // sinon rien
       else  
        { return NULL;};
    };       

    
// récupération de la base locales au noeud noe, pour le temps: temps 
const BaseB & ElemThermi::Gib_elemeca(Enum_dure temps, const Noeud * noe)
   { int nb_noeud = tab_noeud.Taille(); // nombre de noeud de l'élément local
     // on fait une rotation sur les noeuds de l'élément 
     int num=1;
     for (num=1;num<= nb_noeud;num++)
       { if (noe == tab_noeud(num)) break; }
     // récupération des coordonnées locales du noeud
     ElemGeomC0& elem_geom = ElementGeometrique(); // recup de l'élément géométrique correspondant
     // récup de la référence des coordonnées locales du noeud num
     const Coordonnee & coo = elem_geom.PtelemRef()(num);
     // calcul de la dérivée des fonctions au point
     const Mat_pleine& ddphi = elem_geom.Dphi_point(coo);
     const Vecteur& pphi = elem_geom.Phi_point(coo);
     // en fonction du temps on calcul la base naturelle
     switch (temps)
      { case  TEMPS_0: return (met->BaseNat_0(tab_noeud,ddphi,pphi)); break;
        case  TEMPS_t: return (met->BaseNat_t(tab_noeud,ddphi,pphi)); break;
        case  TEMPS_tdt: return (met->BaseNat_tdt(tab_noeud,ddphi,pphi)); break;
       };
     // pour faire taire le compilateur
     cout << "\n erreur, on ne devrait jamais passer ici"
          << "\n ElemThermi::Gib_elemeca(...";
     Sortie(1);
     BaseB *toto = new BaseB(0);
     return (*toto);       
   };
   
// calcul du mini du module d'young équivalent pour tous les points d'intégrations
double ElemThermi::Calcul_mini_E_qui(Enum_dure temps)
 { double mini = 0.; // init
   int ni;
   for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
     { def->ChangeNumInteg(ni);def->Mise_a_jour_data_specif(tabSaveDefDon(ni));
       double E_equi = loiComp->Module_young_equivalent(temps,*def,tabSaveDon(ni));
       if (E_equi < mini) mini = E_equi;
      }
   return mini;   
   };    
   
// calcul du  maxi du module d'young équivalent pour tous les points d'intégrations
double  ElemThermi::Calcul_maxi_E_qui(Enum_dure temps)
 { double maxi = 0.; // init
   int ni;
   for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
     { def->ChangeNumInteg(ni);def->Mise_a_jour_data_specif(tabSaveDefDon(ni));
       double E_equi = loiComp->Module_young_equivalent(temps,*def,tabSaveDon(ni));
       if (E_equi > maxi) maxi = E_equi;
      }
   return maxi;   
   };    

	    
        // cas d'un chargement aero-hydrodynamique, sur les frontières de l'élément 
        // Il y a trois forces: une suivant la direction de la vitesse: de type traînée aerodynamique
        // Fn = poids_volu * fn(V) * S * (normale*u) * u, u étant le vecteur directeur de V (donc unitaire)
        // une suivant la direction normale à la vitesse de type portance
        // Ft = poids_volu * ft(V) * S * (normale*u) * w, w unitaire, normal à V, et dans le plan n et V
        // une suivant la vitesse tangente de type frottement visqueux
        // T = to(Vt) * S * ut, Vt étant la vitesse tangentielle et ut étant le vecteur directeur de Vt
        // retourne  le second membre résultant
        // calcul à l'instant tdt ou t en fonction de la variable atdt
        // coef_mul: est un coefficient multiplicateur global (de tout)
Vecteur& ElemThermi::SM_charge_hydrodyn_E (const double& poidvol,const Tableau <Vecteur>& taphi,int nbne
                                         ,Courbe1D* frot_fluid,const Vecteur& poids
                                         ,Courbe1D* coef_aero_n,int numfront,const double& coef_mul
                                         ,Courbe1D* coef_aero_t,const ParaAlgoControle & ,bool atdt)
   {  
//****************** attention voir w non normé pour les dimensions 2 et 1 ****************   
//****************** attention voir w non normé pour les dimensions 2 et 1 ****************   
//****************** attention voir w non normé pour les dimensions 2 et 1 ****************   
//****************** attention voir w non normé pour les dimensions 2 et 1 ****************   
//****************** attention voir w non normé pour les dimensions 2 et 1 ****************   




    // *** dans l'ensemble du sp on utilise des vecteurs et on gére directement la variance
    // *** sinon c'est intractable et long
    // remarque: l'initialisation du résidu n'est pas fait ici, mais dans les routines appelants
    Vecteur* SM_r;
    int dime = ParaGlob::Dimension();
    bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
    switch (dime)
     {case 3: // cas de la dimension 3
      { // définition du vecteur de retour 
        Vecteur& SM = *((*res_extS)(numfront));SM_r=&SM;
        // dans le cas ou le poid volumique est nul, on retourne sans calcul
        if (poidvol < ConstMath::trespetit) return SM;
        int ni; // compteur globale de point d'integration
        Deformation & defS = *defSurf(numfront); // pour simplifier l'écriture
        for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
         { // -- calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
           // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
           Met_abstraite::Expli ex;
           if (atdt) {ex = (defS.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t(); }
           else      {ex = defS.Cal_explicit_t(premier_calcul); };
           // détermination de la normale à la surface qui  vaut le produit vectoriel des 2 premiers vecteurs
                   // conversion explicite en coordonnées sans variance
           Coordonnee normale = Util::ProdVec_coorBN( (*ex.giB_t)(1), (*ex.giB_t)(2));
 //          cout << "\n normale avant normalisation";
 //          normale.Affiche();
           normale.Normer(); // on norme la normale
  //         cout << "\n normale après normalisation";
  //         normale.Affiche(); Sortie(1);

           // récupération de la vitesse
           const Coordonnee* Vp;
           if (atdt) {Vp = &(defS.VitesseM_tdt());} 
           else      {Vp = &(defS.VitesseM_t());};
           const Coordonnee& V = *Vp; double nV = V.Norme(); 
           // dans le cas où la vitesse est trop faible, le résidu est nul
           if (nV < ConstMath::trespetit) return SM;           
           Coordonnee Vt=V-(V*normale)*normale;// calcul de la vitesse tangente
           double nVt = Vt.Norme(); // norme de la vitesse tangente
           if (nVt < ConstMath::trespetit) { nVt=0.;Vt.Zero();} else {Vt /= nVt;}; // on norme la vitesse tangente
           Coordonnee u= V / nV; // on norme la vitesse, car ensuite on n'a besoin que de la direction unitaire de la vitesse
           double cos_directeur = u * normale; // permet de calculer surface normale à la vitesse
           Coordonnee w(3);
           if (coef_aero_t != NULL) 
             {w = Util::ProdVec_coor(u,Util::ProdVec_coor( u, normale));
              // si w est non nulle on le norme
              double norme_w = w.Norme();
              if (norme_w > ConstMath::trespetit) {w /= norme_w;};              
              };
           // calcul du résidu
           int ix=1; int dimf = 3;   // ne fonctionne qu'en dim 3  
           // cas de T: T = to(Vt) * S * ut, Vt étant la vitesse tangentielle et ut étant le vecteur directeur de Vt
           double coef_frot_fluid=0.;
           if (frot_fluid != NULL) { coef_frot_fluid += frot_fluid->Valeur(nVt) * poids(ni) * (*ex.jacobien_t);}; 
           // cas de Fn = poids_volu * fn(V) * S * (normale*u) * u, u unitaire suivant V
           double coef_n=0.;
           if (coef_aero_n != NULL) 
              { coef_n += poidvol * coef_aero_n->Valeur(nV) *  poids(ni) * (*ex.jacobien_t) * cos_directeur;}; 
           // cas de Ft = poids_volu * ft(V) * S * (normale*u) * w, w unitaire suivant la normal à V
           double coef_t=0.;
           if (coef_aero_t != NULL) 
              { coef_t += poidvol * coef_aero_t->Valeur(nV) *  poids(ni) * (*ex.jacobien_t) * cos_directeur;};
           for (int ne =1; ne<= nbne; ne++) 
             for (int i=1;i<= dimf;i++,ix++)
               { SM(ix) += taphi(ni)(ne) * coef_mul * (coef_frot_fluid * Vt(i) + coef_n * u(i) + coef_t * w(i));};
          } // fin boucle sur les points d'intégrations
        break;   
      }; // fin du cas de la dimension 3
      case 2: // cas de la dimension 2
      { // définition du vecteur de retour 
        Vecteur& SM = *((*res_extA)(numfront));SM_r=&SM;
        // dans le cas ou le poid volumique est nul, on retourne sans calcul
        if (poidvol < ConstMath::trespetit) return SM;
        int ni; // compteur globale de point d'integration
        Deformation & defA = *defArete(numfront); // pour simplifier l'écriture
        for (defA.PremierPtInteg(), ni = 1;defA.DernierPtInteg();defA.NevezPtInteg(),ni++)
         { // -- calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
           // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
           Met_abstraite::Expli ex;
           if (atdt) {ex = (defA.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t(); }
           else      {ex = defA.Cal_explicit_t(premier_calcul); };
           // détermination de la normale à la surface qui  vaut le vecteur perpendiculaire au g1
           // dans le sens directe 
           Coordonnee normale;  normale(1) = - (*ex.giB_t)(1)(2); normale(2) = (*ex.giB_t)(1)(1);
           normale.Normer(); // on norme la normale
           // récupération de la vitesse
           const Coordonnee* Vp;
           if (atdt) {Vp = &(defA.VitesseM_tdt());} 
           else      {Vp = &(defA.VitesseM_t());};
           const Coordonnee& V = *Vp; double nV = V.Norme();
           // dans le cas où la vitesse est trop faible, le résidu est nul
           if (nV < ConstMath::trespetit) return SM;           
           Coordonnee Vt=V-(V*normale)*normale;// calcul de la vitesse tangente
           double nVt = Vt.Norme(); // norme de la vitesse tangente
           if (nVt < ConstMath::trespetit) { nVt=0.;Vt.Zero();} else {Vt /= nVt;}; // on norme la vitesse tangente
           Coordonnee u= V / nV; // on norme la vitesse, car ensuite on n'a besoin que de la direction unitaire de la vitesse
           double cos_directeur = u * normale; // permet de calculer surface normale à la vitesse
           Coordonnee w(dime); // comme on est en 2D w ce calcul directement avec les composantes de u
           if (coef_aero_t != NULL) {w(1) = -u(2); w(2)=u(1);};
           // calcul du résidu
           int ix=1;      
           // cas de T: T = to(Vt) * S * ut, Vt étant la vitesse tangentielle et ut étant le vecteur directeur de Vt
           double coef_frot_fluid=0.;
           if (frot_fluid != NULL) { coef_frot_fluid += frot_fluid->Valeur(nVt) * poids(ni) * (*ex.jacobien_t);}; 
           // cas de Fn = poids_volu * fn(V) * S * (normale*u) * u, u unitaire suivant V
           double coef_n=0.;
           if (coef_aero_n != NULL) 
              { coef_n += poidvol * coef_aero_n->Valeur(nV) *  poids(ni) * (*ex.jacobien_t) * cos_directeur;}; 
           // cas de Ft = poids_volu * ft(V) * S * (normale*u) * w, w unitaire suivant la normal à V
           double coef_t=0.;
           if (coef_aero_t != NULL) 
              { coef_t += poidvol * coef_aero_t->Valeur(nV) *  poids(ni) * (*ex.jacobien_t) * cos_directeur;};
           for (int ne =1; ne<= nbne; ne++) 
             for (int i=1;i<= dime;i++,ix++)
               { SM(ix) += taphi(ni)(ne) * coef_mul * (coef_frot_fluid * Vt(i) + coef_n * u(i) + coef_t * w(i));};
          } // fin boucle sur les points d'intégrations
         } // fin boucle sur les points d'intégrations
        break;   
      case 1: // cas de la dimension 1
      { // définition du vecteur de retour 
        Vecteur& SM = *((*res_extN)(numfront));SM_r=&SM;
        // dans le cas ou le poid volumique est nul, on retourne sans calcul
        if (poidvol < ConstMath::trespetit) return SM;
        // il n'y a pas plusieurs point d'intégration
        // en 1D il ne peut y avoir que deux extrémités, la première à une normale = -1 , la seconde 1
        Coordonnee normale(1); if (numfront == 1) {normale(1)=-1.;} else {normale(1)=1.;};
        // récupération de la vitesse
        const Coordonnee* Vp;
        Noeud* noe = (tabb(numfront)->TabNoeud())(1); // on récupére tout d'abord le noeud fontière
        if (atdt) {Vp = &(def->VitesseM_tdt(noe));}   // on utilise la déformation générale, qui convient pour un noeud
        else      {Vp = &(def->VitesseM_t(noe));};
        const Coordonnee& V = *Vp; double nV = V.Norme();
        // dans le cas où la vitesse est trop faible, le résidu est nul
        if (nV < ConstMath::trespetit) return SM;
        // dans le cas 1D il n'y a qu'un effort de trainée (car pas de vitesse tangente)           
        Coordonnee u= V / nV; // on norme la vitesse, car ensuite on n'a besoin que de la direction unitaire de la vitesse
           double cos_directeur = u * normale; // permet de calculer surface normale à la vitesse
        // calcul du résidu
        int ix=1; int dimf = dime;     
        // cas de Fn = poids_volu * fn(V) * S * (normale*u) * u, u unitaire suivant V
        double coef_n=0.; 
        if (coef_aero_n != NULL) // le poids(1) contiend la section, et il n'y a pas de jacobien, cos directeur = 1
              { coef_n += poidvol * coef_aero_n->Valeur(nV) *  poids(1) ;}; 
        for (int ne =1; ne<= nbne; ne++)  // il n'y qu'un point d'intégration le noeud lui-même
          for (int i=1;i<= dimf;i++,ix++)
               { SM(ix) += taphi(1)(ne) * coef_mul * ( coef_n * u(i) );}; 
        break;   
      }; // fin du cas de la dimension 1 */
      
    }; //-- fin du switch sur la dimension  
    // retour du second membre  
    return *SM_r;  
  };
    
       // idem SM_charge_hydrodyn_E mais  -> implicite, 
	   // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur  
	   // retourne le second membre et la matrice de raideur correspondant
        // coef_mul: est un coefficient multiplicateur global (de tout)
Element::ResRaid ElemThermi::SM_charge_hydrodyn_I (const double& poidvol,const Tableau <Vecteur>& taphi,int nbne
                                         ,Courbe1D* frot_fluid,const Vecteur& poids,DdlElement & ddls
                                         ,Courbe1D* coef_aero_n,int numfront,const double& coef_mul
                                         ,Courbe1D* coef_aero_t,const ParaAlgoControle & pa)
   { 
   
   
//****************** attention voir w non normé pour les dimensions 2 et 1 ****************   
//****************** attention voir w non normé pour les dimensions 2 et 1 ****************   
//****************** attention voir w non normé pour les dimensions 2 et 1 ****************   
//****************** attention voir w non normé pour les dimensions 2 et 1 ****************   
//****************** attention voir w non normé pour les dimensions 2 et 1 ****************   
   
    
    // *** dans l'ensemble du sp on utilise des vecteurs et on gére directement la variance
    // *** sinon c'est intractable et long
    // remarque: l'initialisation du résidu n'est pas fait ici, mais dans les routines appelants
    Element::ResRaid el;               
    int nbddl = ddls.NbDdl();
    bool avec_raid = pa.Var_charge_externe();  // controle
    int dime = ParaGlob::Dimension();
    bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
    switch (dime)
     {case 3: // cas de la dimension 3
      { // dans le cas ou le poid volumique est nul, on retourne sans calcul
        if (poidvol < ConstMath::trespetit)
        {  el.res = (*res_extS)(numfront); el.raid = (*raid_extS)(numfront);       
           return el;
         };
        Vecteur& SM = (*((*res_extS)(numfront)));  //  pour simplifier
        Mat_pleine & KM = (*((*raid_extS)(numfront))); //  "      "       "
        el.res = (*res_extS)(numfront); el.raid = (*raid_extS)(numfront); // pour le retour
        int ni; // compteur globale de point d'integration
        Deformation & defS = *defSurf(numfront); // pour simplifier l'écriture
        for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
         { // -- calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
           // que pour un calcul primaire en implicit mais pour un calcul autre que mécanique
           const Met_abstraite::Impli& ex = defS.Cal_implicit(premier_calcul); 
           // détermination de la normale à la surface qui  vaut le produit vectoriel des 2 premiers vecteurs
           Coordonnee normale = Util::ProdVec_coorBN( (*ex.giB_tdt)(1), (*ex.giB_tdt)(2));
           normale.Normer(); // on norme la normale
           // récupération de la vitesse
           const Coordonnee V = defS.VitesseM_tdt(); double nV = V.Norme(); 
           // dans le cas où la vitesse est trop faible, le résidu est nul
           if (nV < ConstMath::trespetit) return el;           
           Coordonnee Vt=V-(V*normale)*normale;// calcul de la vitesse tangente
           double nVt = Vt.Norme(); // norme de la vitesse tangente
           if (nVt < ConstMath::trespetit) { nVt=0.;Vt.Zero();} else {Vt /= nVt;}; // on norme la vitesse tangente
           Coordonnee u= V / nV; // on norme la vitesse, car ensuite on n'a besoin que de la direction unitaire de la vitesse
           double cos_directeur = u * normale; // permet de calculer surface normale à la vitesse
           Coordonnee w(3);double norme_w=0.;
			  Coordonnee vecNul(3); // un vecteur nul pour l'initialisation
           if (coef_aero_t != NULL) 
             {w = Util::ProdVec_coor(u,Util::ProdVec_coor( u, normale));
              // si w est non nulle on le norme
              norme_w = w.Norme();
              if (norme_w > ConstMath::trespetit) {w /= norme_w;};              
              };
           // calcul éventuel des variations
           static Tableau <Coordonnee> D_Vt,D_u,D_w,D_w_non_normer;
           static Tableau <double> D_coef_frot_fluid,D_coef_n,D_coef_t;
           if (avec_raid)
            { // calcul de la variation de la normale
              // 1) variation du produit vectoriel
              Tableau <Coordonnee > D_pasnormale = 
                  Util::VarProdVect_coor((*ex.giB_tdt)(1).Coor(),(*ex.giB_tdt)(2).Coor(),(*ex.d_giB_tdt));
              // 2) de la normale        
              Tableau <Coordonnee> D_normale = Util::VarUnVect_coor(normale,D_pasnormale,normale.Norme());
              // récupération de la variation de vitesse
              bool ddl_vitesse;
              const Tableau <Coordonnee> D_V = defS.Der_VitesseM_tdt(ddl_vitesse);
              D_Vt.Change_taille(nbddl); // calcul de la variation de la vitesse tangente
              for (int ii=1;ii<=nbddl;ii++)
               {D_Vt(ii)=D_V(ii) - (D_V(ii) * normale + V * D_normale(ii)) * normale
                                 - (V*normale) * D_normale(ii);};
              // variation de u
              D_u.Change_taille(nbddl); D_u = Util::VarUnVect_coor( V, D_V,nV);
              // variation de la norme de V
              Tableau <double> D_nV = Util::VarNorme(D_V,V,nV);
              // variation de w s'il a une norme différente de 0
              D_w.Change_taille(nbddl,vecNul); // !! init obligatoire car peut ne pas rentrer dans les if qui suivent
              if (norme_w > ConstMath::trespetit)
               {if (coef_aero_t != NULL)
                 { D_w_non_normer.Change_taille(nbddl);
					    // tout d'abord la variation du produit vectoriel 
                   Coordonnee k = Util::ProdVec_coor( u, normale);              
                   Tableau <Coordonnee> D_k = Util::VarProdVect_coor(u,normale,D_u,D_normale);
                   D_w_non_normer=Util::VarProdVect_coor(u,k,D_u,D_k);
                   D_w = Util::VarUnVect_coor(w,D_w_non_normer,norme_w);
                  }
               };
              // variation des coefficients
              // cas de T: T = to(Vt) * S * ut, Vt étant la vitesse tangentielle et ut étant le vecteur directeur de Vt
              D_coef_frot_fluid.Change_taille(nbddl);
              if (frot_fluid != NULL) 
                { Courbe1D::ValDer valder = frot_fluid->Valeur_Et_derivee(nVt);
                  for (int i=1;i<=nbddl;i++) 
                    D_coef_frot_fluid(i) = valder.derivee * D_nV(i) * poids(ni) * (*ex.jacobien_tdt)
                                           + valder.valeur * poids(ni) * (*ex.d_jacobien_tdt)(i);
                 }; 
              // cas de Fn = poids_volu * fn(V) * S * (normale*u) * u, u unitaire suivant V
              D_coef_n.Change_taille(nbddl);
              if (coef_aero_n != NULL) 
                { Courbe1D::ValDer valder = coef_aero_n->Valeur_Et_derivee(nV);
                  for (int i=1;i<=nbddl;i++) 
                    D_coef_n(i) = poidvol *  poids(ni) * (
                                  valder.derivee * D_nV(i)  * (*ex.jacobien_tdt) * cos_directeur
                                  + valder.valeur * (*ex.d_jacobien_tdt)(i) * cos_directeur 
                                  + valder.valeur * (*ex.jacobien_tdt) * (D_u(i) * normale + u * D_normale(i))
                                                        );
                 };                                       
              // cas de Ft = poids_volu * ft(V) * S * (normale*u) * w, w unitaire suivant la normal à V
              D_coef_t.Change_taille(nbddl);
              if (coef_aero_t != NULL) 
               { Courbe1D::ValDer valder = coef_aero_t->Valeur_Et_derivee(nV);
                 for (int i=1;i<=nbddl;i++) 
                    D_coef_t(i) = poidvol *  poids(ni) * (
                                  valder.derivee * D_nV(i)  * (*ex.jacobien_tdt) * cos_directeur
                                  + valder.valeur * (*ex.d_jacobien_tdt)(i) * cos_directeur 
                                  + valder.valeur * (*ex.jacobien_tdt) * (D_u(i) * normale + u * D_normale(i))
                                                      );
                };
              };  //-- fin du if (avec_raid)
           // calcul du résidu
           int ix=1; int dimf = 3;   // ne fonctionne qu'en dim 3  
           // cas de T: T = to(Vt) * S * ut, Vt étant la vitesse tangentielle et ut étant le vecteur directeur de Vt
           double coef_frot_fluid=0.;
           if (frot_fluid != NULL) { coef_frot_fluid += frot_fluid->Valeur(nVt) * poids(ni) * (*ex.jacobien_tdt);}; 
           // cas de Fn = poids_volu * fn(V) * S * (normale*u) * u, u unitaire suivant V
           double coef_n=0.;
           if (coef_aero_n != NULL) 
              { coef_n += poidvol * coef_aero_n->Valeur(nV) *  poids(ni) * (*ex.jacobien_tdt) * cos_directeur;}; 
           // cas de Ft = poids_volu * ft(V) * S * (normale*u) * w, w unitaire suivant la normal à V
           double coef_t=0.;
           if (coef_aero_t != NULL) 
              { coef_t += poidvol * coef_aero_t->Valeur(nV) *  poids(ni) * (*ex.jacobien_tdt) * cos_directeur;};
           for (int ne =1; ne<= nbne; ne++) 
             for (int i=1;i<= dimf;i++,ix++)
               { SM(ix) += taphi(ni)(ne) * coef_mul * (coef_frot_fluid * Vt(i) + coef_n * u(i) + coef_t * w(i));
                 if (avec_raid)
                  {for (int j =1; j<= nbddl; j++)  
                     { //cout << "\n taphi(ni)(ne)= " << taphi(ni)(ne) << endl;
							  //cout << "\n D_coef_frot_fluid(j) " << D_coef_frot_fluid(j) << endl;
							  //cout << "\n Vt(i) " << Vt(i) << endl;
							  //cout << "\n D_Vt(j)(i) " << D_Vt(j)(i) << endl;
							  //cout << "\n D_coef_n(j) " << D_coef_n(j) << endl;
							  //cout << "\n u(i) " << u(i) << endl;
							  //cout << "\n D_u(j)(i) " << D_u(j)(i) << endl;
							  //cout << "\n D_coef_t(j) " << D_coef_t(j) << endl;
							  //cout << "\n w(i) " << w(i) << endl;
							  //cout << "\n D_w(j)(i) " << D_w(j)(i) << endl;
							  
							  KM(ix,j) -= taphi(ni)(ne) * coef_mul *(
                                     D_coef_frot_fluid(j) * Vt(i) + coef_frot_fluid * D_Vt(j)(i)
                                   + D_coef_n(j) * u(i) + coef_n * D_u(j)(i) 
                                   + D_coef_t(j) * w(i) + coef_t * D_w(j)(i)
                                                   ); 
                      };
                  };
                };  
          } // fin boucle sur les points d'intégrations
        break;   
      }; // fin du cas de la dimension 3
      case 2: // cas de la dimension 2
      { // dans le cas ou le poid volumique est nul, on retourne sans calcul
        if (poidvol < ConstMath::trespetit)
        {  el.res = (*res_extS)(numfront); el.raid = (*raid_extS)(numfront);       
           return el;
         };
        Vecteur& SM = (*((*res_extA)(numfront)));  //  pour simplifier
        Mat_pleine & KM = (*((*raid_extA)(numfront))); //  "      "       "
        el.res = (*res_extA)(numfront); el.raid = (*raid_extA)(numfront); // pour le retour
        int ni; // compteur globale de point d'integration
        Deformation & defA = *defArete(numfront); // pour simplifier l'écriture
        for (defA.PremierPtInteg(), ni = 1;defA.DernierPtInteg();defA.NevezPtInteg(),ni++)
         { // -- calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
           // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
           const Met_abstraite::Impli& ex = defA.Cal_implicit(premier_calcul); 
           // détermination de la normale à la surface qui  vaut le vecteur perpendiculaire au g1
           // dans le sens directe 
           Coordonnee normale;  normale(1) = -(*ex.giB_tdt)(1)(2); normale(2) = (*ex.giB_tdt)(1)(1);
           normale.Normer(); // on norme la normale
           // récupération de la vitesse
           const Coordonnee V = defA.VitesseM_tdt(); double nV = V.Norme(); 
           // dans le cas où la vitesse est trop faible, le résidu est nul
           if (nV < ConstMath::trespetit) return el;           
           Coordonnee Vt=V-(V*normale)*normale;// calcul de la vitesse tangente
           double nVt = Vt.Norme(); // norme de la vitesse tangente
           if (nVt < ConstMath::trespetit) { nVt=0.;Vt.Zero();} else {Vt /= nVt;}; // on norme la vitesse tangente
           Coordonnee u= V / nV; // on norme la vitesse, car ensuite on n'a besoin que de la direction unitaire de la vitesse
           double cos_directeur = u * normale; // permet de calculer surface normale à la vitesse
           Coordonnee w(dime); // comme on est en 2D, w ce calcul directement avec les composantes de u
           if (coef_aero_t != NULL) {w(1) = -u(2); w(2)=u(1);};
           // calcul éventuel des variations
           static Tableau <Coordonnee> D_Vt,D_u,D_w;
           static Tableau <double> D_coef_frot_fluid,D_coef_n,D_coef_t;
           if (avec_raid)
            { // calcul de la variation de la normale
              // 1) variation du produit vectoriel
              Tableau <Coordonnee > D_pasnormale(nbddl);
              for (int i=1;i<=nbddl;i++)
                {D_pasnormale(i)(1) = -(*ex.d_giB_tdt)(i)(1)(2); D_pasnormale(i)(2) = (*ex.d_giB_tdt)(i)(1)(1);};
              // 2) de la normale        
              Tableau <Coordonnee> D_normale = Util::VarUnVect_coor(normale,D_pasnormale,normale.Norme());
              // récupération de la variation de vitesse
              bool ddl_vitesse;
              const Tableau <Coordonnee> D_V = defA.Der_VitesseM_tdt(ddl_vitesse);
              D_Vt.Change_taille(nbddl); // calcul de la variation de la vitesse tangente
              for (int ii=1;ii<=nbddl;ii++)
               {D_Vt(ii)=D_V(ii) - (D_V(ii) * normale + V * D_normale(ii)) * normale
                                 - (V*normale) * D_normale(ii);};
              // variation de u
              D_u.Change_taille(nbddl); D_u = Util::VarUnVect_coor( V, D_V,nV);
              // variation de la norme de V
              Tableau <double> D_nV = Util::VarNorme(D_V,V,nV);
              // variation de w
              D_w.Change_taille(nbddl);
              if (coef_aero_t != NULL)
               { for (int i=1;i<=nbddl;i++) {D_w(i)(1)=-D_u(i)(2);D_w(i)(2)=D_u(i)(1);}; 
                } 
              // variation des coefficients
              // cas de T: T = to(Vt) * S * ut, Vt étant la vitesse tangentielle et ut étant le vecteur directeur de Vt
              D_coef_frot_fluid.Change_taille(nbddl);
              if (frot_fluid != NULL) 
                { Courbe1D::ValDer valder = frot_fluid->Valeur_Et_derivee(nVt);
                  for (int i=1;i<=nbddl;i++) 
                    D_coef_frot_fluid(i) = valder.derivee * D_nV(i) * poids(ni) * (*ex.jacobien_tdt)
                                           + valder.valeur * poids(ni) * (*ex.d_jacobien_tdt)(i);
                 }; 
              // cas de Fn = poids_volu * fn(V) * S * (normale*u) * u, u unitaire suivant V
              D_coef_n.Change_taille(nbddl);
              if (coef_aero_n != NULL) 
                { Courbe1D::ValDer valder = coef_aero_n->Valeur_Et_derivee(nV);
                  for (int i=1;i<=nbddl;i++) 
                    D_coef_n(i) = poidvol *  poids(ni) * (
                                  valder.derivee * D_nV(i)  * (*ex.jacobien_tdt) * cos_directeur
                                  + valder.valeur * (*ex.d_jacobien_tdt)(i) * cos_directeur 
                                  + valder.valeur * (*ex.jacobien_tdt) * (D_u(i) * normale + u * D_normale(i))
                                                        );
                 };                                       
              // cas de Ft = poids_volu * ft(V) * S * (normale*u) * w, w unitaire suivant la normal à V
              D_coef_t.Change_taille(nbddl);
              if (coef_aero_t != NULL) 
               { Courbe1D::ValDer valder = coef_aero_t->Valeur_Et_derivee(nV);
                 for (int i=1;i<=nbddl;i++) 
                    D_coef_t(i) = poidvol *  poids(ni) * (
                                  valder.derivee * D_nV(i)  * (*ex.jacobien_tdt) * cos_directeur
                                  + valder.valeur * (*ex.d_jacobien_tdt)(i) * cos_directeur 
                                  + valder.valeur * (*ex.jacobien_tdt) * (D_u(i) * normale + u * D_normale(i))
                                                      );
                };
              };  //-- fin du if (avec_raid)
           // calcul du résidu
           int ix=1;    
           // cas de T: T = to(Vt) * S * ut, Vt étant la vitesse tangentielle et ut étant le vecteur directeur de Vt
           double coef_frot_fluid=0.;
           if (frot_fluid != NULL) { coef_frot_fluid += frot_fluid->Valeur(nVt) * poids(ni) * (*ex.jacobien_tdt);}; 
           // cas de Fn = poids_volu * fn(V) * S * (normale*u) * u, u unitaire suivant V
           double coef_n=0.;
           if (coef_aero_n != NULL) 
              { coef_n += poidvol * coef_aero_n->Valeur(nV) *  poids(ni) * (*ex.jacobien_tdt) * cos_directeur;}; 
           // cas de Ft = poids_volu * ft(V) * S * (normale*u) * w, w unitaire suivant la normal à V
           double coef_t=0.;
           if (coef_aero_t != NULL) 
              { coef_t += poidvol * coef_aero_t->Valeur(nV) *  poids(ni) * (*ex.jacobien_tdt) * cos_directeur;};
           for (int ne =1; ne<= nbne; ne++) 
             for (int i=1;i<= dime;i++,ix++)
               { SM(ix) += taphi(ni)(ne) * coef_mul * (coef_frot_fluid * Vt(i) + coef_n * u(i) + coef_t * w(i));
                 if (avec_raid)
                  {for (int j =1; j<= nbddl; j++)  
                     {KM(ix,j) -= taphi(ni)(ne) * coef_mul * (
                                     D_coef_frot_fluid(j) * Vt(i) + coef_frot_fluid * D_Vt(j)(i)
                                   + D_coef_n(j) * u(i) + coef_n * D_u(j)(i) 
                                   + D_coef_t(j) * w(i) + coef_t * D_w(j)(i)
                                                   ); 
                      };
                  };
                };  
          } // fin boucle sur les points d'intégrations
        break; 
        };  
      case 1: // cas de la dimension 1
      { // dans le cas ou le poid volumique est nul, on retourne sans calcul
        if (poidvol < ConstMath::trespetit)
        {  Element::ResRaid el;    
           el.res = (*res_extN)(numfront); 
           el.raid = (*raid_extN)(numfront);       
           return el;
         };
        Vecteur& SM = (*((*res_extN)(numfront)));  //  pour simplifier
        Mat_pleine & KM = (*((*raid_extN)(numfront))); //  "      "       "
        el.res = (*res_extN)(numfront); el.raid = (*raid_extN)(numfront); // pour le retour
        // il n'y a pas plusieurs point d'intégration
        // en 1D il ne peut y avoir que deux extrémités, la première à une normale = -1 , la seconde 1
        Coordonnee normale(1); if (numfront == 1) {normale(1)=-1.;} else {normale(1)=1.;};
        // récupération de la vitesse
        Noeud* noe = (tabb(numfront)->TabNoeud())(1); // on récupére tout d'abord le noeud fontière
        const Coordonnee& V = def->VitesseM_tdt(noe); double nV = V.Norme();
        // dans le cas où la vitesse est trop faible, le résidu est nul
        if (nV < ConstMath::trespetit) return el;
        // dans le cas 1D il n'y a qu'un effort de trainée (car pas de vitesse tangente)           
        Coordonnee u= V / nV; // on norme la vitesse, car ensuite on n'a besoin que de la direction unitaire de la vitesse
        // u en fait = 1 ou -1 d'où variation = 0
        // calcul éventuel des variations
        static Tableau <double> D_coef_n;
        if (avec_raid)
            { // calcul de la variation de la normale
              // récupération de la variation de vitesse
              bool ddl_vitesse;
              const Tableau <Coordonnee> D_V = def->Der_VitesseM_tdt(ddl_vitesse,noe);
              // variation de la norme de V
              Tableau <double> D_nV = Util::VarNorme(D_V,V,nV);
              // variation des coefficients
              // cas de Fn = poids_volu * fn(V) * S * (normale*u) * u, u unitaire suivant V
              D_coef_n.Change_taille(nbddl);
              if (coef_aero_n != NULL) 
                { Courbe1D::ValDer valder = coef_aero_n->Valeur_Et_derivee(nV);
                  for (int i=1;i<=nbddl;i++) 
                    D_coef_n(i) = poidvol *  poids(1) * valder.derivee * D_nV(i);
                 };                                       
              };  //-- fin du if (avec_raid)
        // calcul du résidu
        int ix=1; int dimf = dime;     
        // cas de Fn = poids_volu * fn(V) * S * (normale*u) * u, u unitaire suivant V
        double coef_n=0.; 
        if (coef_aero_n != NULL) // le poids(1) contiend la section, et il n'y a pas de jacobien, cos directeur = 1
              { coef_n += poidvol * coef_mul * coef_aero_n->Valeur(nV) *  poids(1) ;}; 
        for (int ne =1; ne<= nbne; ne++)  // il n'y qu'un point d'intégration le noeud lui-même
          for (int i=1;i<= dimf;i++,ix++)
               { SM(ix) += taphi(1)(ne) * ( coef_n * u(i) );
                 if (avec_raid)
                  {for (int j =1; j<= nbddl; j++)  
                     {KM(ix,j) -= taphi(1)(ne) * coef_mul * (D_coef_n(j) * u(i));}; 
                  };
               }; 
        break;   
      }; // fin du cas de la dimension 1 */
      
    }; //-- fin du switch sur la dimension  
    // retour   
    return el;  
  };