// 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-2021 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 "Deformation.h"
#include "ConstMath.h"
#include "ParaGlob.h"
#include "MathUtil.h"
#include "Tenseur3.h"
#include "Tenseur2.h"
#include "Tenseur1.h"
# include "TypeConsTens.h"
#include "NevezTenseur.h"


//================== suite du fichier Deformation.cp  ======================
// calcul des valeurs propres, les projections propres, et les variations d'un tenseurs BH 
// utilisé dans  Cal_vite_rota_objectif
void Deformation::Val_et_projection_prop_tenseur(const TenseurBH & B_BH,const Tableau <TenseurBH * >& d_B_BH
                      ,Coordonnee& kii,Tableau <TenseurBH* >& Palpha_BH,bool variation
                      ,Tableau <Coordonnee >& d_kii, Tableau <Tableau <TenseurBH* > >& d_Palpha_BH
                      ,int& cas_valPropre)
{ double unTier=1./3.; double unSixieme=1./6.;// quelques constantes
  int dima=Abs(B_BH.Dimension()); // récup de la dimension de B_BH
  kii = B_BH.ValPropre(cas_valPropre); // calcul des valeurs propres de base
  // ----------------------------
  // | calcul des projecteurs   |
  // ----------------------------

  // -- traitement en fonction de la dimension
  switch (dima)
   { case 1:    // ---- cas de la dimension un ------
       // calcul des projections propres
      {switch (cas_valPropre)
        { case 1: // cas où l'on a 1 valeurs propres  correctement calculée
           { // calcul du projecteur qui ici est l'identité en BH donc delta_i^J 
             (*Palpha_BH(1))= IdBH1;
             break;
           }
          default:  
            cout << "\n erreur dans les calculs des valeurs propres, on ne peut pas calculer les projecteurs " 
                 << "\n Deformation::Val_et_projection_prop_tenseur(...";
            Sortie(1);
        };
       break;
      }   
     case 2:    // ---- cas de la dimension deux ------
       // calcul des projections propres
       // il nous faut distinguer les différents cas suivant que l'on a 1 ou 2 valeurs propres distinctes
      {switch (cas_valPropre)
        { case 1: // cas où l'on a 2 valeurs propres différentes
           { // on vérifie cependant la différence entre les valeurs propres
             if ( Dabs(kii(1)-kii(2)) > ConstMath::pasmalpetit)
              { // calcul des deux projections           
                (*Palpha_BH(1))=(1./(kii(1)-kii(2)))*(B_BH- kii(2)*IdBH2);
                (*Palpha_BH(2))=(1./(kii(2)-kii(1)))*(B_BH- kii(1)*IdBH2);
                break;
               }
             else   
             // sinon on continue sur le cas 0, 
              {cas_valPropre = 0; }; // pour les variations ce qui évitera un test   
            }
          case 0: // toutes les valeurs propres sont égales
           { //  -- tout d'abord redimentionne Palpha_BH
             // calcul du projecteur qui ici est l'identité en BH donc delta_i^J 
             (*Palpha_BH(1))= IdBH2;
             break;
           }
          default:  
            cout << "\n erreur dans les calculs des valeurs propres, on ne peut pas calculer les projecteurs " 
                 << "\n Deformation::Val_et_projection_prop_tenseur(...";
            Sortie(1);
        };
       break;
      }   
      case 3:  // ---- cas de la dimension trois ------
       // calcul des projections propres
       // il nous faut distinguer les différents cas suivant que l'on a 1 ou 2 ou 3 valeurs propres distinctes
      {switch (cas_valPropre)
        { case 1: // cas où l'on a 3 valeurs propres différentes
           { // calcul des trois projections           
             (*Palpha_BH(1))=(1./((kii(1)-kii(2))*(kii(1)-kii(3))))*(B_BH- kii(2)*IdBH3)*(B_BH- kii(3)*IdBH3);
             (*Palpha_BH(2))=(1./((kii(2)-kii(1))*(kii(2)-kii(3))))*(B_BH- kii(1)*IdBH3)*(B_BH- kii(3)*IdBH3);
             (*Palpha_BH(3))=(1./((kii(3)-kii(2))*(kii(3)-kii(1))))*(B_BH- kii(1)*IdBH3)*(B_BH- kii(2)*IdBH3);
             break;
            }
          case 0: // toutes les valeurs propres sont égales
           { //  -- tout d'abord redimentionne Palpha_BH
             // calcul du projecteur qui ici est l'identité en BH donc delta_i^J 
             (*Palpha_BH(1))= IdBH3;
             break;
           }
          case 2: case 3: // cas où les deux premières valeurs propres sont égales ou les deux dernières
           { // calcul des deux projecteurs
             (*Palpha_BH(1))=(1./(kii(1)-kii(2)))*(B_BH- kii(2)*IdBH3);
             (*Palpha_BH(2))=(1./(kii(2)-kii(1)))*(B_BH- kii(1)*IdBH3);
             break;
           }
          default:  
            cout << "\n erreur dans les calculs des valeurs propres, on ne peut pas calculer les projecteurs " 
                 << "\n Deformation::Val_et_projection_prop_tenseur(...";
            Sortie(1);
        };
       break;
      }   
      default:  
        cout << "\n dimension (premier) non encore pris en compte: " << dima
             << "\n Deformation::Val_et_projection_prop_tenseur(...";
        Sortie(1);
   };

  // ------------------------------------------------------------------
  // | calcul des variations des valeurs propres et des projecteurs   |
  // ------------------------------------------------------------------  
                 
  if(variation)
   {int nbvar=d_B_BH.Taille(); // récup de la dimension
    //  --- choix suivant la dimension
    switch (dima)
    { case 1:  // %%%%%  dimension 1 %%%%%%%
       {// 1) --- calcul de la variation des valeurs propres 
        // calcul de la variation de la valeur propre
        for (int i= 1;i<=nbvar;i++)
         { const Tenseur1BH& dB_BH = *((const Tenseur1BH*)(d_B_BH(i)) ) ;  // pour simplifier l'ecriture     
           // variation de la valeur propre
           d_kii(i)(1)= dB_BH.Trace();
           // 2) -- maintenant calcul des variations des projections propres 
           Tableau <Tenseur1BH* > & dPalpha_BH= (Tableau <Tenseur1BH* > &) d_Palpha_BH(i); // pour simplifier
           switch (cas_valPropre) // il n'y a que le cas = 1 qui est normal
            { case 1: // cas où l'on a 1 valeur propre normalement calculée
               { // le projecteur est l'identité donc la variation est nulle
                 dPalpha_BH(1)->Inita(0.);
                 break;
                }
             }; // -- fin du switch sur cas_valPropre; pas de cas défault car ça a déjà été vu au début 
          }; // -- fin de la boucle sur les ddl : nbvar
        break;
        }; // -- fin du case 1 pour les variations
      case 2:  // %%%%%  dimension 2 %%%%%%%
       {// 1) --- calcul de la variation des valeurs propres 
        // calcul des 2 invariants de base
        double Ib_B=B_BH.Trace();double IIb_B=0.5*B_BH.II();     
        // calcul des 2 invariants principaux de Cayley-Hamilton 
        // I_A = A_i^i, II_A=(I_A)**2/2-trace(A*A)/2, 
        double I_B_BH=Ib_B; double II_B_BH=0.5*Ib_B*Ib_B-IIb_B;
        // calcul des variations des valeurs propres
            // def de deux tenseurs de travail
        Tenseur2BH dB_BHetB_BH(B_BH);Tenseur2BH B_BHetdB_BH(B_BH);
            // def d'un tenseur de travail
        Tenseur2BH Z_ij_BH;
        for (int i= 1;i<=nbvar;i++)
         { const Tenseur2BH& dB_BH = *((const Tenseur2BH*)(d_B_BH(i)) ) ;  // pour simplifier l'ecriture     
           // variation du premier invariant de base
           double dIb_B=dB_BH.Trace();
           // 2) -- maintenant calcul des variations des valeurs propres et des projections propres   
           Tableau <Tenseur2BH* > & dPalpha_BH= (Tableau <Tenseur2BH* > &) d_Palpha_BH(i); // pour simplifier
           switch (cas_valPropre)
            { case 0: // toutes les valeurs propres sont égales
               { // on ne sait pas calculer la variation de la  valeur propre, on la met à 0
                 d_kii(i).Zero();
                 // le projecteur est l'identité on ne connait pas également sa variation, on la met à 0
                 dPalpha_BH(1)->Inita(0.);
                 break;
               }
              case 1: // cas où l'on a 2 valeurs propres différentes
               { // cas où les deux valeurs propres sont différentes
                 // calcul de la variation du second invariant de base
                 dB_BHetB_BH=dB_BH * B_BH;B_BHetdB_BH=B_BH * dB_BH;
                 double dIIb_B=0.5*( dB_BHetB_BH.Trace()+ B_BHetdB_BH.Trace());
                 // variation du second  invariant principal de Cayley-Hamilton (le premier est déjà calculé)          
                 double dII_B_BH = dIb_B*Ib_B  - dIIb_B;
                 // variation des valeurs propres
                 for (int j=1;j<=2;j++)
                       d_kii(i)(j)=(dIb_B * kii(j) - dII_B_BH )/(2.*kii(j)-Ib_B);
                 // calcul un seul tenseur intermédiaire qui sert deux fois
                 Z_ij_BH= (1./ (kii(1)-kii(2))) * dB_BH;
                 // variation des projections
                 (*dPalpha_BH(1)) = Z_ij_BH -  (d_kii(i)(2)*IdBH2) / (kii(1)-kii(2))
                                        - (*Palpha_BH(1)) * ((d_kii(i)(1) -d_kii(i)(2)) / (kii(1)-kii(2)));
                 (*dPalpha_BH(2)) = - Z_ij_BH -  (d_kii(i)(1)*IdBH2) / (kii(2)-kii(1))
                                        - (*Palpha_BH(2)) * ((d_kii(i)(2) -d_kii(i)(1)) / (kii(2)-kii(1)));
                 break;
               }
            }; // -- fin du switch sur cas_valPropre; pas de cas défault car ça a déjà été vu au début 
          }; // -- fin de la boucle sur les ddl : nbvar
        break;
        }; // -- fin du case 2 pour les variations
      case 3:    // %%%%%  dimension 3 %%%%%%%
       {// 1) --- calcul de la variation des valeurs propres 
        // calcul des 3 invariants de base
        double Ib_B=B_BH.Trace();double IIb_B=0.5*B_BH.II();double IIIb_B=unTier*B_BH.III();        
        // calcul des 3 invariants principaux de Cayley-Hamilton 
        // I_A = A_i^i, II_A=(I_A)**2/2-trace(A*A)/2, 
        // III_A=1/3*trace ((A*A)*A)-I_A*1/2*trace (A*A)+1/6*(I_A)**3
        double I_B_BH=Ib_B; double II_B_BH=0.5*Ib_B*Ib_B-IIb_B;
        double III_B_BH=IIIb_B-Ib_B*IIb_B+unSixieme*Ib_B*Ib_B*Ib_B;
        // calcul des variations des valeurs propres
            // def de deux tenseurs de travail
        Tenseur3BH dB_BHetB_BH(B_BH);Tenseur3BH B_BHetdB_BH(B_BH);
            // def de 6 tenseurs de travail
        Tableau <Tenseur3BH> Z_ij_BH(6);
        for (int i= 1;i<=nbvar;i++)
         { const Tenseur3BH& dB_BH = *((const Tenseur3BH*)(d_B_BH(i)) ) ;  // pour simplifier l'ecriture     
           // 2) -- maintenant calcul des variations des projections propres   
           Tableau <Tenseur3BH* > & dPalpha_BH= (Tableau <Tenseur3BH* > &) d_Palpha_BH(i); // pour simplifier
           switch (cas_valPropre)
            { case 0: // toutes les valeurs propres sont égales
               { // variation de la valeur propre // non !!!(on ne sait pas la calculer, on la met à 0.
                 //d_kii(i)(1)=0.;
                 // dans ce cas cela veut dire que le tenseur est sphérique, on suppose qu'il reste sphérique
                 // sinon ce n'est pas possible, dans ce cas on a toujours I_A = 3 * lambda 
                 // d'où d_lambda = 1/3 * d_I_A
                 d_kii(i)(1)=1./3. * dB_BH.Trace();
                 // le projecteur est l'identité donc la variation est nulle(c'est pas vrai mais on ne sans sert pas  après)
                 dPalpha_BH(1)->Inita(0.);
                 break;
               }
              case 1: // cas où l'on a 3 valeurs propres différentes
               { // variation des invariants de base
                 double dIb_B=dB_BH.Trace();
                 dB_BHetB_BH=dB_BH * B_BH;B_BHetdB_BH=B_BH * dB_BH;
                 double dIIb_B=0.5*( dB_BHetB_BH.Trace()+ B_BHetdB_BH.Trace());
                 double dIIIb_B=unTier*( (dB_BHetB_BH*B_BH).Trace()+(B_BHetdB_BH*B_BH).Trace()+ (B_BH*B_BHetdB_BH).Trace());
                 // variation des 2 invariants principaux de Cayley-Hamilton (le premier est déjà calculé)          
                 double dII_B_BH = dIb_B*Ib_B  - dIIb_B;
                 double dIII_B_BH = dIIIb_B - dIb_B*IIb_B - Ib_B*dIIb_B + 0.5*dIb_B*Ib_B*Ib_B;
                 // variation des 3 valeurs propres
                 for (int j=1;j<=3;j++)
                    d_kii(i)(j)=(dIb_B * kii(j) * kii(j) - dII_B_BH * kii(j) + dIII_B_BH)
                                /(3.*kii(j)*kii(j)-2.*I_B_BH*kii(j)+II_B_BH);
                 // calcul des tenseurs Z_ij_BH
                 // Z_12 -> 1 , Z_21 -> 2 , Z_13 -> 3 , Z_31 -> 4 , Z_23 -> 5 , Z_32 -> 6 ,
                 for (int m=1;m<=6;m++)
                  { int i1=ccdex.iii(m); int j1=ccdex.jjj(m); // ccdex est défini dans deformation
                    Z_ij_BH(m)=(1./((kii(i1)-kii(j1))*(kii(i1)-kii(j1))))
                               *((dB_BH-(d_kii(i))(j1)*IdBH3) * (kii(i1)-kii(j1)) 
                               - (B_BH-kii(j1)*IdBH3)*((d_kii(i))(i1)-(d_kii(i))(j1)) );
                  };
                 // variation des projections
                 (*dPalpha_BH(1)) = (Z_ij_BH(1) * (B_BH- kii(3)*IdBH3)) / (kii(1)-kii(3)) 
                                    + ((B_BH- kii(2)*IdBH3) * Z_ij_BH(3)) / (kii(1)-kii(2));
                 (*dPalpha_BH(2)) = (Z_ij_BH(2) * (B_BH- kii(3)*IdBH3)) / (kii(2)-kii(3)) 
                                    + ((B_BH- kii(1)*IdBH3) * Z_ij_BH(5)) / (kii(2)-kii(1));
                 (*dPalpha_BH(3)) = (Z_ij_BH(6) * (B_BH- kii(1)*IdBH3)) / (kii(3)-kii(1)) 
                                    + ((B_BH- kii(2)*IdBH3) * Z_ij_BH(4)) / (kii(3)-kii(2));
                 break;
               }
              case 2: case 3: // cas où deux valeurs propres sont égales
               { // variation des invariants de base
                 double dIb_B=dB_BH.Trace();
                 dB_BHetB_BH=dB_BH * B_BH;B_BHetdB_BH=B_BH * dB_BH;
                 double dIIb_B=0.5*( dB_BHetB_BH.Trace()+ B_BHetdB_BH.Trace());
                 double dIIIb_B=unTier*( (dB_BHetB_BH*B_BH).Trace()+(B_BHetdB_BH*B_BH).Trace()+ (B_BH*B_BHetdB_BH).Trace());
                 // variation des 2 invariants principaux de Cayley-Hamilton (le premier est déjà calculé)          
                 double dII_B_BH = dIb_B*Ib_B  - dIIb_B;
                 double dIII_B_BH = dIIIb_B - dIb_B*IIb_B - Ib_B*dIIb_B + 0.5*dIb_B*Ib_B*Ib_B;
                 // variation des 2 valeurs propres, comme le dénominateur de la formule générale est nulle (et oui comme
                 // dans le cas des 3 valeurs propres identiques) on utilise de nouvelles relations:
                 // on fait l'hypothèse que l'évolution va conservé le fait qu'il y a deux valeurs propres identiques  
                 // soit a la valeur propre double et b l'autre valeur propre, on a:
                 // I_A = 2*a+b, et II_A= 2*a*b + a^2, d'où en différentiant: d_I_A=2*d_a + d_b => d_b=d_I_A-2*d_a
                 // et d_II_a=2*a*d_b + 2*b*d_a + 2*a*d_a = 2*a*(d_I_A-2*d_a) + 2*b*d_a + 2*a*d_a d'où
                 // d_II_a - 2*a*d_I_A = 2*(-2*a + b + a)*d_a = 2*(b-a)*d_a => d_a = (d_II_a - 2*a*d_I_A)/(2*(b-a))
                 // on remarque que la formule n'est valable que pour a différent de b !!
                 if (cas_valPropre==2)
                  { // les deux premières valeurs propres sont identiques
                    double r = (dII_B_BH - 2.*kii(1)*dIb_B)/(2.*(kii(2)-kii(1))); // variable inter
                    d_kii(i)(1) = r;
                    d_kii(i)(2) = dIb_B - 2.* r;
                  }
                 else // cas ou les deux dernières valeurs propres sont identiques
                  { // les deux premières valeurs propres sont identiques
                    double r = (dII_B_BH - 2.*kii(2)*dIb_B)/(2.*(kii(1)-kii(2))); // variable inter
                    d_kii(i)(2) = r;
                    d_kii(i)(1) = dIb_B - 2.* r;
                  };
                 
        //        vieux calcul général qui ne marche pas quand deux valeurs propres sont identiques !! 
        //         for (int j=1;j<=2;j++)
        //            d_kii(i)(j)=(dIb_B * kii(j) * kii(j) - dII_B_BH * kii(j) + dIII_B_BH)
        //                         /(3.*kii(j)*kii(j)-2.*I_B_BH*kii(j)+II_B_BH);
        
                 // calcul un seul tenseur intermédiaire qui sert deux fois
                 Z_ij_BH(1)= (1./ (kii(1)-kii(2))) * dB_BH;
                 // variation des projections
                 (*dPalpha_BH(1)) = Z_ij_BH(1) -  (d_kii(i)(2)*IdBH3) / (kii(1)-kii(2))
                                    - (*Palpha_BH(1)) * ((d_kii(i)(1) -d_kii(i)(2)) / (kii(1)-kii(2)));
                 (*dPalpha_BH(2)) = - Z_ij_BH(1) -  (d_kii(i)(1)*IdBH3) / (kii(2)-kii(1))
                                    - (*Palpha_BH(2)) * ((d_kii(i)(2) -d_kii(i)(1)) / (kii(2)-kii(1)));
                 break;
               }
             }; // -- fin du switch sur cas_valPropre; pas de cas défault car ça a déjà été vu au début 
          }; // -- fin de la boucle sur les ddl : nbvar
        break;
        }; // -- fin du case 3 pour les variations
      default:  
        cout << "\n dimension non encore pris en compte: " << dima
             << "\n Deformation::Val_et_projection_prop_tenseur(...";
        Sortie(1);
      } }          
};           
   
 
// dimensionement ou vérif des dimensions des variables intermédiaires
// utilisées pour les mesures cumulées et/ou log 
void Deformation::DimensionementVarLog(int dima,bool avec_var,int nbvar)
{ if (dimensionnement_tableaux==0)
   { // les tableaux internes
     B_BH_tr=NevezTenseurBH(dima);
     ki.Change_dim(dima);
     Palpha_BH_tr = ExisteTabTenseur_BH_tr(dima,Palpha_BH_tr,dima); // allocation et vérif de Palpha_BH_tr
     if (avec_var)
      {d_B_BH_tr = ExisteTabTenseur_BH_tr(nbvar,d_B_BH_tr,dima); // allocation et vérif de d_B_BH_tr
       d_ki = ExisteTabCoor_tr(nbvar,d_ki,dima); // vérif et allocation de d_ki
       d_Palpha_BH_tr = ExisteTabTabTenseur_BH_tr(dima,d_Palpha_BH_tr,nbvar); // vérif et allocation de d_Palpha_BH_tr
       dimensionnement_tableaux=2; 
       }
     else
      {dimensionnement_tableaux=1;};   
    }
  else if ((dimensionnement_tableaux==1)&& (avec_var))
    // cas où le dimensionnement explicit a déjà été fait mais pas celui des variations
    {d_B_BH_tr = ExisteTabTenseur_BH_tr(nbvar,d_B_BH_tr,dima); // allocation et vérif de d_B_BH_tr
     d_ki = ExisteTabCoor_tr(nbvar,d_ki,dima); // vérif et allocation de d_ki
     d_Palpha_BH_tr = ExisteTabTabTenseur_BH_tr(dima,d_Palpha_BH_tr,nbvar); // vérif et allocation de d_Palpha_BH_tr
     dimensionnement_tableaux=2; 
    };
};
                    
//   ----------------------------------------------            
//      $$$$$$$$$$$$$  private: $$$$$$$$$$$$$$$$
//   ----------------------------------------------

    // verif si tableau de tenseur est bien dimensionné, ou alloue si nécessaire
Tableau <TenseurBH * >* Deformation::ExisteTabTenseur_BH_tr(int nbvar,Tableau <TenseurBH * >* tab_Tens,int dima) const
 { // on vérifie que tab_Tens pointe sur quelque chose de correcte
   if (tab_Tens  != NULL)
     { int tails= tab_Tens->Taille();
       if (tails== nbvar)
          { bool correcte = true;
            for (int i=1;i<=nbvar;i++)
               if ((*tab_Tens)(i)->Dimension()!= dima) correcte = false;
                   if (correcte) return tab_Tens; // c'est déjà ok
           }; 
        // sinon le pointage est mauvais
        tab_Tens=NULL;
        // on ne supprime pas l'objet pointé, car il peut servir à d'autre instance de déformation
        // et comme il est en static, globalement la place est faible, 
      };      
   // maintenant on regarde l'allocation éventuelle
   if (tab_Tens  == NULL)
      { // on parcours la liste pour voire s'il existe un bon tableau de tenseur
        list <Tableau <TenseurBH * > >::iterator il,ilfin= list_var_tens_BH.end();
        bool trouver=false;
        for(il=list_var_tens_BH.begin();il!=ilfin;il++)
           { if ((*il).Taille()== nbvar) 
                {bool correcte = true;
                 for (int i=1;i<=nbvar;i++)
                    if ((*il)(i)->Dimension() != dima) {correcte = false;break;};
                 if (correcte) // cas où l'on réutilise une place static qui existe déjà 
                       { tab_Tens = &(*il); return tab_Tens;};
                  }  
             };  
        if (!trouver) // cas où il faut créer une nouvelle place
           { tab_Tens = new Tableau <TenseurBH* >; tab_Tens->Change_taille(nbvar);
             for (int i=1;i<=nbvar;i++)
                  (*tab_Tens)(i)=NevezTenseurBH(dima);
             list_var_tens_BH.push_back(*tab_Tens);
            };
         }   
    // retour
    return tab_Tens;
 };
  
    // verif si le tableau de tableau de tenseur est bien dimensionné, ou alloue si nécessaire
    // dima: dimension des tenseur et du premier tableau, nbvar: dimension du tableau extérieur
Tableau <Tableau <TenseurBH* > >*  Deformation::ExisteTabTabTenseur_BH_tr
                           (int dima,Tableau <Tableau <TenseurBH* > >* tab_Tens,int nbvar) const
 { // on vérifie que tab_Tens pointe sur quelque chose de correcte
   if (tab_Tens  != NULL)
     { bool correcte = true;
       int tails= tab_Tens->Taille();
       if (tails!=nbvar) 
         {correcte = false;}
       else  
         {for (int j=1;j<=nbvar;j++)
             {int tailsj=(*tab_Tens)(j).Taille();
              if (tailsj != dima)
                { correcte = false;break;}
              else  
                { for (int i=1;i<=dima;i++)
                      if ((*tab_Tens)(j)(i)->Dimension()!= dima) 
                         {correcte = false;break;};
                   if(!correcte) break;                               
                  };
               };
           };           
       if (correcte) 
         {return tab_Tens;}; // c'est déjà ok
       // sinon le pointage est mauvais
       tab_Tens=NULL;
       // on ne supprime pas l'objet pointé, car il peut servir à d'autre instance de déformation
       // et comme il est en static, globalement la place est faible, 
      };      
   // maintenant on regarde l'allocation éventuelle
   if (tab_Tens  == NULL)
      { // on parcours la liste pour voire s'il existe un bon tableau de tableau de tenseur
        list <Tableau <Tableau <TenseurBH* > > >::iterator il,ilfin= list_var_var_tens_BH.end();
        for(il=list_var_var_tens_BH.begin();il!=ilfin;il++)
           { if ((*il).Taille()== nbvar) 
              { bool correcte = true;
                for (int j=1;j<=nbvar;j++)
                  { if ((*il)(j).Taille()== dima) 
                     { for (int i=1;i<=dima;i++)
                         if ((*il)(j)(i)->Dimension() != dima) {correcte = false;break;};
                       };
                    if(!correcte) break;
                    };                            
                if (correcte) // cas où l'on réutilise une place static qui existe déjà 
                   { tab_Tens = &(*il); return tab_Tens;};
               }  
             };
       };      
   // si on arrive ici c'est que l'on n'a rien trouver, il faut créer un nouveau tableau       
   tab_Tens = new Tableau <Tableau <TenseurBH* > >; tab_Tens->Change_taille(nbvar);
   for (int ic=1;ic<=nbvar;ic++) 
     { (*tab_Tens)(ic).Change_taille(dima);
       for (int it=1;it<=dima;it++)  (*tab_Tens)(ic)(it)=NevezTenseurBH(dima);
      }; 
   list_var_var_tens_BH.push_back(*tab_Tens);
    // retour
    return tab_Tens;
 };

    // verif si tableau de vecteur est bien dimensionné, ou alloue si nécessaire
Tableau <Coordonnee>* Deformation::ExisteTabCoor_tr(int nbvar,Tableau <Coordonnee>* tab_Coor,int dima) const 
 { // on vérifie que tab_Coor pointe sur quelque chose de correcte
   if (tab_Coor  != NULL)
     { int tails= tab_Coor->Taille();
       if (tails== nbvar)
          { bool correcte = true;
            for (int i=1;i<=nbvar;i++)
               if ((*tab_Coor)(i).Dimension()!= dima) correcte = false;
                   if (correcte) return tab_Coor; // c'est déjà ok
           }; 
        // sinon le pointage est mauvais
        tab_Coor=NULL;
        // on ne supprime pas l'objet pointé, car il peut servir à d'autre instance de déformation
        // et comme il est en static, globalement la place est faible, 
      };      
   // maintenant on regarde l'allocation éventuelle
   if (tab_Coor  == NULL)
      { // on parcours la liste pour voire s'il existe un bon tableau de vecteur
        list <Tableau <Coordonnee> >::iterator il,ilfin= list_tab_coor.end();
        bool trouver=false;
        for(il=list_tab_coor.begin();il!=ilfin;il++)
           { if ((*il).Taille()== nbvar) 
                {bool correcte = true;
                 for (int i=1;i<=nbvar;i++)
                    if ((*il)(i).Dimension() != dima) {correcte = false;break;};
                 if (correcte) // cas où l'on réutilise une place static qui existe déjà 
                       { tab_Coor = &(*il); return tab_Coor;};
                  }  
             };  
        if (!trouver) // cas où il faut créer une nouvelle place
           { tab_Coor = new Tableau <Coordonnee>; tab_Coor->Change_taille(nbvar);
             for (int i=1;i<=nbvar;i++)
                  (*tab_Coor)(i).Change_dim(dima);
             list_tab_coor.push_back(*tab_Coor);
            };
         }   
    // retour
    return tab_Coor;
 };