
// 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 "Util.h"
#include "Coordonnee3.h"
#include "Coordonnee2.h"
#include "Coordonnee1.h"
#include "ConstMath.h"
#include "MathUtil.h"

#ifndef  Util_H_deja_inclus

//#ifndef MISE_AU_POINT
//  inline 
//#endif
//// PRODUIT VECTORIEL DE DEUX VECTEURS EN COORDONNEES ASBOLUS
//Vecteur Util::ProdVec( const Vecteur & v1, const Vecteur & v2)
//  { 
//   #ifdef MISE_AU_POINT
//      if (v1.Taille() != v2.Taille())
//			{ cout << "\nErreur : la taille des deux vecteur n'est pas identiques  !\n";
//			  cout << "Util::ProdVec(Vecteur & v1,Vecteur & v2) \n";
//			  Sortie(1);
//			};
//	  if (v1.Taille() != 3)
//			{ cout << "\nErreur : la taille des deux vecteur est differente de trois  !\n";
//			  cout << "Util::ProdVec(Vecteur & v1,Vecteur & v2) \n";
//			  Sortie(1);
//			};		
//  #endif 
//  Vecteur v(3);
//  v(1) = v1(2) * v2(3) - v1(3) * v2(2);
//  v(2) = v1(3) * v2(1) - v1(1) * v2(3);   
//  v(3) = v1(1) * v2(2) - v1(2) * v2(1);   
//  return v;   
//  };
//
  
#ifndef MISE_AU_POINT
  inline 
#endif
// idem en coordonnées absolues avec le type Coordonnee
Coordonnee Util::ProdVec_coor( const Coordonnee & v1, const Coordonnee & v2)
  { 
   #ifdef MISE_AU_POINT
      if (v1.Dimension() != v2.Dimension())
			{ cout << "\nErreur : la taille des deux Coordonnees n'est pas identiques  !\n";
			  cout << "Util::ProdVec(const Coordonnee & v1, const Coordonnee & v2) \n";
			  Sortie(1);
			};
	  if (v1.Dimension() != 3)
			{ cout << "\nErreur : la taille des deux vecteur est differente de trois  !\n";
			  cout << "Util::ProdVec(Vecteur & v1,Vecteur & v2) \n";
			  Sortie(1);
			};		
  #endif 
  Coordonnee v(3);
  v(1) = v1(2) * v2(3) - v1(3) * v2(2);
  v(2) = v1(3) * v2(1) - v1(1) * v2(3);   
  v(3) = v1(1) * v2(2) - v1(2) * v2(1);   
  return v;   
  };
  
// idem en coordonnées contravariantes avec le type CoordonneeH
#ifndef MISE_AU_POINT
  inline 
#endif
CoordonneeH Util::ProdVec_coorH( const CoordonneeH & v1, const CoordonneeH & v2)
  { 
   #ifdef MISE_AU_POINT
      if (v1.Dimension() != v2.Dimension())
			{ cout << "\nErreur : la taille des deux CoordonneesH n'est pas identiques  !\n";
			  cout << "Util::ProdVec_coorH( const CoordonneeH & v1, const CoordonneeH & v2) \n";
			  Sortie(1);
			};
	  if (v1.Dimension() != 3)
			{ cout << "\nErreur : la taille des deux CoordonneesH est differente de trois  !\n";
			  cout << "Util::ProdVec_coorH( const CoordonneeH & v1, const CoordonneeH & v2) \n";
			  Sortie(1);
			};		
  #endif 
  CoordonneeH v(3);
  v(1) = v1(2) * v2(3) - v1(3) * v2(2);
  v(2) = v1(3) * v2(1) - v1(1) * v2(3);   
  v(3) = v1(1) * v2(2) - v1(2) * v2(1);   
  return v;   
  };

// idem en coordonnées covariantes avec le type CoordonneeB
#ifndef MISE_AU_POINT
  inline 
#endif
CoordonneeB Util::ProdVec_coorB( const CoordonneeB & v1, const CoordonneeB & v2)
  { 
   #ifdef MISE_AU_POINT
      if (v1.Dimension() != v2.Dimension())
			{ cout << "\nErreur : la taille des deux CoordonneesB n'est pas identiques  !\n";
			  cout << "Util::ProdVec_coorB( const CoordonneeB & v1, const CoordonneeB & v2) \n";
			  Sortie(1);
			};
	  if (v1.Dimension() != 3)
			{ cout << "\nErreur : la taille des deux CoordonneesB est differente de trois  !\n";
			  cout << "Util::ProdVec_coorB( const CoordonneeB & v1, const CoordonneeB & v2) \n";
			  Sortie(1);
			};		
  #endif 
  CoordonneeB v(3);
  v(1) = v1(2) * v2(3) - v1(3) * v2(2);
  v(2) = v1(3) * v2(1) - v1(1) * v2(3);   
  v(3) = v1(1) * v2(2) - v1(2) * v2(1);   
  return v;   
  };
    
// idem en coordonnées covariantes avec le type CoordonneeB, et en retour un coordonnee normal
#ifndef MISE_AU_POINT
  inline 
#endif
Coordonnee Util::ProdVec_coorBN( const CoordonneeB & v1, const CoordonneeB & v2)
  { 
   #ifdef MISE_AU_POINT
      if (v1.Dimension() != v2.Dimension())
			{ cout << "\nErreur : la taille des deux CoordonneesB n'est pas identiques  !\n";
			  cout << "Util::ProdVec_coorBN( const CoordonneeB & v1, const CoordonneeB & v2) \n";
			  Sortie(1);
			};
	  if (v1.Dimension() != 3)
			{ cout << "\nErreur : la taille des deux CoordonneesB est differente de trois  !\n";
			  cout << "Util::ProdVec_coorBN( const CoordonneeB & v1, const CoordonneeB & v2) \n";
			  Sortie(1);
			};		
  #endif 
  Coordonnee v(3);
  v(1) = v1(2) * v2(3) - v1(3) * v2(2);
  v(2) = v1(3) * v2(1) - v1(1) * v2(3);   
  v(3) = v1(1) * v2(2) - v1(2) * v2(1);   
  return v;   
  };
    
//#ifndef MISE_AU_POINT
//  inline 
//#endif
//// CALCUL DU DETERMINANT DE TROIS VECTEURS
//double Util::Determinant( const Vecteur & v1, const Vecteur & v2, const Vecteur & v3)
//  {
//   #ifdef MISE_AU_POINT
//      if ((v1.Taille() != 3) || (v2.Taille() != 3) || (v3.Taille() != 3))
//			{ cout << "\nErreur : la taille des trois vecteur est differente de trois  !\n";
//			  cout << "Util::Determinant(Vecteur & v1,Vecteur & v2,Vecteur & v3)\n";
//			  Sortie(1);
//			}  	
//  #endif 
//  return (ProdVec(v1,v2) * v3);
//  };
  
// CALCUL DU DETERMINANT DE TROIS VECTEURS en coordonnées covariantes
#ifndef MISE_AU_POINT
  inline 
#endif
double Util::DeterminantB( const CoordonneeB & v1, const CoordonneeB & v2, const CoordonneeB & v3)
  {
   #ifdef MISE_AU_POINT
      if ((v1.Dimension() != 3) || (v2.Dimension() != 3) || (v3.Dimension() != 3))
			{ cout << "\nErreur : la taille des trois Coordonee est differente de trois  !\n";
			  cout << "Util::DeterminantB(CoordonneeB & v1,CoordonneeB & v2,CoordonneeB & v3)\n";
			  Sortie(1);
			}  	
  #endif 
  return ((ProdVec_coorB(v1,v2)).ScalBB(v3));
  };
  
// CALCUL DU DETERMINANT DE TROIS VECTEURS  en coordonnées absolues avec le type Coordonnee
#ifndef MISE_AU_POINT
  inline 
#endif
double Util::Determinant( const Coordonnee & v1, const Coordonnee & v2, const Coordonnee & v3)
  {
   #ifdef MISE_AU_POINT
      if ((v1.Dimension() != 3) || (v2.Dimension() != 3) || (v3.Dimension() != 3))
			{ cout << "\nErreur : la taille des trois Coordonee est differente de trois  !\n";
			  cout << "Util::Determinant(Coordonnee & v1,Coordonnee & v2,Coordonnee & v3)\n";
			  Sortie(1);
			}  	
  #endif 
  return (ProdVec_coor(v1,v2) * v3);
  };

//#ifndef MISE_AU_POINT
//  inline 
//#endif
//// CALCUL DU DETERMINANT DE DEUX VECTEURS
//double Util::Determinant( const Vecteur & v1, const Vecteur & v2)
//  {
//   #ifdef MISE_AU_POINT
//      if (((v1.Taille() != 3) && (v1.Taille() != 2)) 
//           || ((v2.Taille() != 3) && (v2.Taille() != 2)) )
//		{ cout << "\nErreur : la dimension des deux vecteur est differente de trois et deux  !\n";
//	      cout << "Util::Determinant(Vecteur & v1,Vecteur & v2)\n";
//		  Sortie(1);
//		}  	
//  #endif
//  double res; 
//  switch (v1.Taille())
//     { case 2 : res = v1(1) * v2(2) - v1(2) * v2(1); break;
//       case 3 : res = !(ProdVec(v1,v2)); break;
//       default : cout << "erreur Determinant(Vecteur & v1,Vecteur & v2)\n";
//		       Sortie(1);
//	  }	        
//  return res;
//  };
//  

// CALCUL DU DETERMINANT DE DEUX VECTEURS en coordonnées covariantes
#ifndef MISE_AU_POINT
  inline 
#endif
double Util::DeterminantB( const CoordonneeB & v1, const CoordonneeB & v2)
  {
   #ifdef MISE_AU_POINT
      if (((v1.Dimension() != 3) && (v1.Dimension() != 2)) 
           || ((v2.Dimension() != 3) && (v2.Dimension() != 2)) )
		{ cout << "\nErreur : la dimension des deux Coordonnees est differente de trois et deux  !\n";
	      cout << "Util::DeterminantB(CoordonneeB & v1,CoordonneeB & v2)\n";
		  Sortie(1);
		}  	
  #endif
  double res; 
  switch (v1.Dimension())
     { case 2 : res = v1(1) * v2(2) - v1(2) * v2(1); break;
       case 3 : res = (ProdVec_coorB(v1,v2)).Norme(); break;
       default : cout << "erreur DeterminantB(CoordonneeB & v1,CoordonneeB & v2)\n";
		       Sortie(1);
	  }	        
  return res;
  };
  
// CALCUL DU DETERMINANT DE DEUX VECTEURS en coordonnées absolues avec le type Coordonnee
#ifndef MISE_AU_POINT
  inline 
#endif
double Util::Determinant( const Coordonnee & v1, const Coordonnee & v2)
  {
   #ifdef MISE_AU_POINT
      if (((v1.Dimension() != 3) && (v1.Dimension() != 2)) 
           || ((v2.Dimension() != 3) && (v2.Dimension() != 2)) )
		{ cout << "\nErreur : la dimension des deux Coordonnees est differente de trois et deux  !\n";
	      cout << "Util::Determinant(Coordonnee & v1,Coordonnee & v2)\n";
		  Sortie(1);
		}  	
  #endif
  double res; 
  switch (v1.Dimension())
     { case 2 : res = v1(1) * v2(2) - v1(2) * v2(1); break;
       case 3 : res = (ProdVec_coor(v1,v2)).Norme(); break;
       default : cout << "erreur Determinant(Coordonnee & v1,Coordonnee & v2)\n";
		       Sortie(1);
	  }	        
  return res;
  };


//#ifndef MISE_AU_POINT
//  inline 
//#endif
//// CALCUL DU DETERMINANT DE UN VECTEUR
//double Util::Determinant( const Vecteur & v1)
//  {    
//  return v1.Norme(); //v1 * v1;
//  };
//

// CALCUL DU DETERMINANT DE UN VECTEUR en coordonnée covariante
#ifndef MISE_AU_POINT
  inline 
#endif
double Util::DeterminantB( const CoordonneeB & v1)
  {    
  return v1.Norme(); //v1 * v1;
  };

// CALCUL DU DETERMINANT DE UN VECTEUR  en coordonnées absolues avec le type Coordonnee
#ifndef MISE_AU_POINT
  inline 
#endif
double Util::Determinant( const Coordonnee & v1)
  {    
  return v1.Norme(); //v1 * v1;
  };

#ifndef MISE_AU_POINT
  inline 
#endif
// calcul de la variation d'un vecteur unitaire connaissant la variation du
// vecteur non norme
// v : le vecteur non norme, Dv : la variation de v, nor : la norme de v
// en retour : la variation de vecteur 
Vecteur Util::VarUnVect( const Vecteur & v,  const Vecteur & Dv, double  nor)
  { double nor3 = nor * nor * nor;
    Vecteur De = Dv / nor - v * ( (v * Dv) / nor3);
    return De;   
   };

#ifndef MISE_AU_POINT
  inline 
#endif
// idem avec coordonnee en retour
Coordonnee Util::VarUnVect_coor( const Coordonnee & v,  const Coordonnee & Dv, double  nor)
  { double nor3 = nor * nor * nor;
    Coordonnee De = Dv / nor - v * ( (v * Dv) / nor3);
    return De;   
   };
   

#ifndef MISE_AU_POINT
  inline 
#endif
// idem avec coordonneeB en in-out 
  //!!!!!!!!!!!!! très important: il doit s'agir de vecteur exprimé dans un repère orthonormé
  // ceci est vrai quelque soit la variance affichée: car ici on ne prend pas en compte la variation
  //  d'une métrique associée à un repère non orthonormé
  // ex: les variations des gi sont ok car en fait les gi représentent les coorddonées dans un repère absolu
CoordonneeB Util::VarUnVect_coorB( const CoordonneeB & v,  const CoordonneeB & Dv, double  nor)
  { double nor3 = nor * nor * nor;
    CoordonneeB De = Dv / nor - v * ( (v.ScalBB(Dv)) / nor3);
    return De;   
   };

#ifndef MISE_AU_POINT
  inline
#endif
// idem avec coordonneeH en in-out
  //!!!!!!!!!!!!! très important: il doit s'agir de vecteur exprimé dans un repère orthonormé
  // ceci est vrai quelque soit la variance affichée: car ici on ne prend pas en compte la variation
  //  d'une métrique associée à un repère non orthonormé
  // ex: les variations des gi sont ok car en fait les gi représentent les coorddonées dans un repère absolu
CoordonneeH Util::VarUnVect_coorH( const CoordonneeH & v,  const CoordonneeH & Dv, double  nor)
  { double nor3 = nor * nor * nor;
    CoordonneeH De = Dv / nor - v * ( (v.ScalHH(Dv)) / nor3);
    return De;
   };

#ifndef MISE_AU_POINT
  inline 
#endif
// calcul du tableau de variation d'un vecteur unitaire connaissant le tableau de variation du
// vecteur non norme
// v : le vecteur non norme, Dv : la variation de v, nor : la norme de v
// en retour : le tableau de variation 
Tableau <Vecteur> Util::VarUnVect( const Vecteur & v, const  Tableau <Vecteur>& Dv, double  nor)
  { Tableau <Vecteur> De(Dv.Taille());
    double nor3 = nor * nor * nor;
    int DvTaille=Dv.Taille();
    for (int i=1;i<= DvTaille;i++)
      { De(i) = Dv(i) / nor - v * ( (v * Dv(i)) / nor3);
       }
    return De;   
   };
// idem avec des coordonnees   
#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <Coordonnee> Util::VarUnVect_coor( const Coordonnee & v, const  Tableau <Coordonnee>& Dv, double  nor)
  { Tableau <Coordonnee> De(Dv.Taille());
    double nor3 = nor * nor * nor;
    int DvTaille=Dv.Taille();
    for (int i=1;i<= DvTaille;i++)
      { De(i) = Dv(i) / nor - v * ( (v * Dv(i)) / nor3);
       }
    return De;   
   };
// idem avec des coordonneesB   
#ifndef MISE_AU_POINT
  inline 
#endif
  //!!!!!!!!!!!!! très important: il doit s'agir de vecteur exprimé dans un repère orthonormé
  // ceci est vrai quelque soit la variance affichée: car ici on ne prend pas en compte la variation
  //  d'une métrique associée à un repère non orthonormé
  // ex: les variations des gi sont ok car en fait les gi représentent les coorddonées dans un repère absolu
Tableau <CoordonneeB> Util::VarUnVect_coorB( const CoordonneeB & v, const  Tableau <CoordonneeB>& Dv, double  nor)
  { Tableau <CoordonneeB> De(Dv.Taille());
    double nor3 = nor * nor * nor;
    int DvTaille=Dv.Taille();
    for (int i=1;i<= DvTaille;i++)
      { De(i) = Dv(i) / nor - v * ( ( v.ScalBB(Dv(i)) ) / nor3);
       }
    return De;   
   };
	
#ifndef MISE_AU_POINT
  inline 
#endif
// idem  et le tableau de retour passé en paramètre
// calcul du tableau de variation d'un vecteur unitaire connaissant le tableau de variation du
// vecteur non norme
// v : le vecteur non norme, Dv : la variation de v, nor : la norme de v
// en retour : le tableau de variation 
Tableau <Vecteur>& Util::VarUnVect( const Vecteur & v, const  Tableau <Vecteur>& Dv, double  nor, Tableau <Vecteur>& retour)
  { retour.Change_taille(Dv.Taille());
    double nor3 = nor * nor * nor;
    int DvTaille=Dv.Taille();
    for (int i=1;i<= DvTaille;i++)
      { retour(i) = Dv(i) / nor - v * ( (v * Dv(i)) / nor3);
       }
	 return retour;	 
   };
// idem avec des coordonnees   
#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <Coordonnee>& Util::VarUnVect_coor( const Coordonnee & v, const  Tableau <Coordonnee>& Dv, double  nor, Tableau <Coordonnee>& retour )
  { retour.Change_taille(Dv.Taille());
    double nor3 = nor * nor * nor;
    int DvTaille=Dv.Taille();
    for (int i=1;i<= DvTaille;i++)
      { retour(i) = Dv(i) / nor - v * ( (v * Dv(i)) / nor3);
       }
	 return retour;	 
   };
// idem avec des coordonneesB   
#ifndef MISE_AU_POINT
  inline 
#endif
  //!!!!!!!!!!!!! très important: il doit s'agir de vecteur exprimé dans un repère orthonormé
  // ceci est vrai quelque soit la variance affichée: car ici on ne prend pas en compte la variation
  //  d'une métrique associée à un repère non orthonormé
  // ex: les variations des gi sont ok car en fait les gi représentent les coorddonées dans un repère absolu
Tableau <CoordonneeB>& Util::VarUnVect_coorB( const CoordonneeB & v, const  Tableau <CoordonneeB>& Dv, double  nor,Tableau <CoordonneeB>& retour)
  { retour.Change_taille(Dv.Taille());
    double nor3 = nor * nor * nor;
    int DvTaille=Dv.Taille();
    for (int i=1;i<= DvTaille;i++)
      { retour(i) = Dv(i) / nor - v * ( ( v.ScalBB(Dv(i)) ) / nor3);
       }
	 return retour;	 
   };
// idem avec des coordonnees   
  // la variation du vecteur est supposé se trouver dans le premier vecteur de la base 
#ifndef MISE_AU_POINT
  inline 
#endif
  //!!!!!!!!!!!!! très important: il doit s'agir de vecteur exprimé dans un repère orthonormé
  // ceci est vrai quelque soit la variance affichée: car ici on ne prend pas en compte la variation
  //  d'une métrique associée à un repère non orthonormé
  // ex: les variations des gi sont ok car en fait les gi représentent les coorddonées dans un repère absolu
Tableau <Coordonnee>& Util::VarUnVect_coorBN( const CoordonneeB & v, const Tableau <BaseB>& Dv, double  nor,Tableau <Coordonnee>& retour)
  { retour.Change_taille(Dv.Taille());
    double nor3 = nor * nor * nor;
    int DvTaille=Dv.Taille();
    for (int i=1;i<= DvTaille;i++)
      { //retour(i) = (Dv(i) / nor - v * ( ( v.ScalBB(Dv(i)(1)) ) / nor3)).Coor_const();
		 const Coordonnee & dv = Dv(i).Coordo(1); // pour simplifier
		 retour(i) = (dv / nor - v.Coor_const() * ( v.Coor_const() * dv) / nor3);
		  
       }
	 return retour;	 
   };
	
	
	

//#ifndef MISE_AU_POINT
//  inline 
//#endif
//// calcul de la variation d'un produit vectoriel
//// vi et Dvi les vecteurs du produit vectoriel et leurs variations
//Tableau <Vecteur> Util::VarProdVect( const Vecteur & v1, const Vecteur & v2,
//                  const Tableau <Vecteur >& Dv1, const Tableau <Vecteur >& Dv2)
// { 
//   #ifdef MISE_AU_POINT
//      if ((v1.Taille() != v2.Taille()) || (Dv1.Taille() != Dv2.Taille())
//           || (Dv1.Taille() == 0) 
//           || (v1.Taille() != Dv1(1).Taille()) || (Dv1(1).Taille() != Dv1(2).Taille()) )
//		{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
//	      cout << "Util::VarProdVect(Vecteur & v1,Vecteur & v2," <<
//                 "Tableau <Vecteur >& Dv1,Tableau <Vecteur >& Dv2)\n";
//		  Sortie(1);
//		}  	
//   #endif
//   int Dv1Taille = Dv1.Taille();
//   Tableau <Vecteur> De(Dv1Taille);
//   for (int i=1;i<= Dv1Taille;i++)
//      { De(i) = ProdVec(Dv1(i),v2) + ProdVec(v1,Dv2(i));
//       }
//   return De;      
// };
 
#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <Coordonnee> Util::VarProdVect_coor( const Coordonnee & v1, const Coordonnee & v2,
                  const Tableau <Coordonnee >& Dv1, const Tableau <Coordonnee >& Dv2)
 { 
   #ifdef MISE_AU_POINT
      if ((v1.Dimension() != v2.Dimension()) || (Dv1.Taille() != Dv2.Taille())
           || (Dv1.Taille() == 0) 
           || (v1.Dimension() != Dv1(1).Dimension()) || (Dv1(1).Dimension() != Dv1(2).Dimension()) )
		{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
	      cout << "Util::VarProdVect_coor(Coordonnee & v1,Coordonnee & v2," <<
                 "Tableau <Coordonnee >& Dv1,Tableau <Coordonnee >& Dv2)\n";
		  Sortie(1);
		}  	
   #endif
   int Dv1taille=Dv1.Taille();
   Tableau <Coordonnee> De(Dv1taille);
   for (int i=1;i<= Dv1taille;i++)
      { De(i) = ProdVec_coor(Dv1(i),v2) + ProdVec_coor(v1,Dv2(i));
       }
   return De;      
 };

// idem que le précédent module mais avec un retour en coordonnée, et un tableau de BaseB
#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <Coordonnee> Util::VarProdVect_coor(const Coordonnee & v1, const Coordonnee & v2,
                  const Tableau <BaseB >& Dv)
 { 
   #ifdef MISE_AU_POINT
      if ((v1.Dimension() != v2.Dimension()) || (Dv(1)(1).Dimension() != v1.Dimension())
           || (Dv.Taille() == 0)) 
		{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
	      cout << "Util::VarProdVect( int toto, const Vecteur & v1, const Vecteur & v2,"
               << " const Tableau <BaseB >& Dv) \n";
		  Sortie(1);
		}  	
   #endif
   int DvTaille = Dv.Taille();
   Tableau <Coordonnee> De(DvTaille);
   for (int i=1;i<= DvTaille;i++)
      { De(i) = ProdVec_coor(Dv(i).Coordo(1),v2) + ProdVec_coor(v1,Dv(i).Coordo(2));
       }
   return De;      
 };

// idem que le précédent module mais avec un in-out  en coordonnéeB, et un tableau de BaseB
#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <CoordonneeB> Util::VarProdVect_coorB(const CoordonneeB & v1, const CoordonneeB & v2,
                  const Tableau <BaseB >& Dv)
 { 
   #ifdef MISE_AU_POINT
      if ((v1.Dimension() != v2.Dimension()) || (Dv(1)(1).Dimension() != v1.Dimension())
           || (Dv.Taille() == 0)) 
		{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
	      cout << "Util::VarProdVect( int toto, const Vecteur & v1, const Vecteur & v2,"
               << " const Tableau <BaseB >& Dv) \n";
		  Sortie(1);
		}  	
   #endif
   int DvTaille = Dv.Taille();
   Tableau <CoordonneeB> De(DvTaille);
   for (int i=1;i<= DvTaille;i++)
      { De(i) = ProdVec_coorB(Dv(i)(1),v2) + ProdVec_coorB(v1,Dv(i)(2));
       }
   return De;      
 };


// idem mais avec le retour passé en paramètre
#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <Coordonnee>& Util::VarProdVect_coor( const Coordonnee & v1, const Coordonnee & v2,
                  const Tableau <Coordonnee >& Dv1, const Tableau <Coordonnee >& Dv2,Tableau <Coordonnee>& retour)
 { 
   #ifdef MISE_AU_POINT
      if ((v1.Dimension() != v2.Dimension()) || (Dv1.Taille() != Dv2.Taille())
           || (Dv1.Taille() == 0) 
           || (v1.Dimension() != Dv1(1).Dimension()) || (Dv1(1).Dimension() != Dv1(2).Dimension()) )
		{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
	      cout << "Util::VarProdVect_coor(Coordonnee & v1,Coordonnee & v2," <<
                 "Tableau <Coordonnee >& Dv1,Tableau <Coordonnee >& Dv2)\n";
		  Sortie(1);
		}  	
   #endif
   int Dv1taille=Dv1.Taille();
	retour.Change_taille(Dv1taille);
   for (int i=1;i<= Dv1taille;i++)
      { retour(i) = ProdVec_coor(Dv1(i),v2) + ProdVec_coor(v1,Dv2(i));
       }
   return retour;      
 };

#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <Coordonnee>& Util::VarProdVect_coor(const Coordonnee & v1, const Coordonnee & v2,
                  const Tableau <BaseB >& Dv,Tableau <Coordonnee>& retour)
 { 
   #ifdef MISE_AU_POINT
      if ((v1.Dimension() != v2.Dimension()) || (Dv(1)(1).Dimension() != v1.Dimension())
           || (Dv.Taille() == 0)) 
		{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
	      cout << "Util::VarProdVect( int toto, const Vecteur & v1, const Vecteur & v2,"
               << " const Tableau <BaseB >& Dv) \n";
		  Sortie(1);
		}  	
   #endif
   int DvTaille = Dv.Taille();
	retour.Change_taille(DvTaille);
   for (int i=1;i<= DvTaille;i++)
      { retour(i) = ProdVec_coor(Dv(i).Coordo(1),v2) + ProdVec_coor(v1,Dv(i).Coordo(2));
       }
   return retour;      
 };

#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <CoordonneeB>& Util::VarProdVect_coorB(const CoordonneeB & v1, const CoordonneeB & v2,
                  const Tableau <BaseB >& Dv,Tableau <CoordonneeB>& retour)
 { 
   #ifdef MISE_AU_POINT
      if ((v1.Dimension() != v2.Dimension()) || (Dv(1)(1).Dimension() != v1.Dimension())
           || (Dv.Taille() == 0)) 
		{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
	      cout << "Util::VarProdVect( int toto, const Vecteur & v1, const Vecteur & v2,"
               << " const Tableau <BaseB >& Dv) \n";
		  Sortie(1);
		}  	
   #endif
   int DvTaille = Dv.Taille();
	retour.Change_taille(DvTaille);
   for (int i=1;i<= DvTaille;i++)
      { retour(i) = ProdVec_coorB(Dv(i)(1),v2) + ProdVec_coorB(v1,Dv(i)(2));
       }
   return retour;      
 };

#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <Coordonnee>& Util::VarProdVect_coorBN(const CoordonneeB & v1, const CoordonneeB & v2,
                  const Tableau <BaseB >& Dv,Tableau <Coordonnee>& retour)
 { 
   #ifdef MISE_AU_POINT
      if ((v1.Dimension() != v2.Dimension()) || (Dv(1)(1).Dimension() != v1.Dimension())
           || (Dv.Taille() == 0)) 
		{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
	      cout << "Util::VarProdVect( int toto, const Vecteur & v1, const Vecteur & v2,"
               << " const Tableau <BaseB >& Dv) \n";
		  Sortie(1);
		}  	
   #endif
   int DvTaille = Dv.Taille();
	retour.Change_taille(DvTaille);
   for (int i=1;i<= DvTaille;i++)
      { retour(i) = ProdVec_coorBN(Dv(i)(1),v2) + ProdVec_coorBN(v1,Dv(i)(2));
       }
   return retour;      
 };
 

#ifndef MISE_AU_POINT
  inline 
#endif
// retourne le numero du ddl recherche identifie par en, dans le tableau passé en paramètre
// s'il existe sinon 0
int Util::Existe(const Tableau<Ddl >& tab_ddl,Enum_ddl en)
  { static int posiddl; // initialisée par le système à zéro
    int tab_ddl_Taille = tab_ddl.Taille();
     // on commence par tester a partir de la derniere position sauvegardee
    for (int i=posiddl+1;i<= tab_ddl_Taille;i++)
      if (en == tab_ddl(i).Id_nom())
          { posiddl=i; return i; }
     // on continue la recherche   a partir du debut
    for (int i=1;i<= posiddl;i++)
      if (en == tab_ddl(i).Id_nom())
          { posiddl=i; return i; }
     // sinon retour 0     
    return 0;     
  };
   
#ifndef MISE_AU_POINT
  inline 
#endif
// calcul du produit mixte des  vecteurs d'une base naturelle en coordonnée covariante
double Util::ProduitMixte(const BaseB & tab_v )
 { int nbvecteur= tab_v.NbVecteur(); // le nombre de vecteur
   // on utilise le fait que les vecteurs de la base sont uniquement 1D 2D ou 3D pour optimiser
   switch (nbvecteur)
   { case 1: // on calcul la norme du vecteur qui est donc toujours positive
          return DeterminantB(tab_v(1));         
          break;
     case 2: // cas de deux vecteurs, on ramène la norme du produit vectoriel
          return DeterminantB(tab_v(1),tab_v(2));      
          break;
     case 3: // produit mixte normal qui peut-être négatif
          return DeterminantB(tab_v(1),tab_v(2),tab_v(3));   
          break;
     default: 
       cout << "\n nombre de vecteur de la base non correcte : " << nbvecteur
            << "\n double Util::ProduitMixte(const BaseB & tab_v )";
       Sortie(1);
    }
    // pour le compilo
    Sortie(1);
    return 0.; // ne doit jamais passer ici   
  };

  // calcul de la variation de la norme d'un vecteur, connaissant la variation du vecteur,  
  // le vecteur, et sa norme
#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <double> Util::VarNorme( const Tableau <Coordonnee >& Dv,const Coordonnee& V,const double& norme)
 { int nbddl = Dv.Taille();
   Tableau <double> tab_Nv(nbddl);
   if (norme < ConstMath::trespetit) {return tab_Nv;}; // retour dans le cas où la vitesse est nulle
   double unsurnorme = 1./norme;   
   for (int i=1;i<=nbddl;i++)
    { tab_Nv(i) = (V*Dv(i))*unsurnorme;};
   return tab_Nv;
  };

// idem mais avec le type vecteur                 
#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <double> Util::VarNorme( const Tableau <Vecteur >& Dv,const Vecteur& V,const double& norme)
 { int nbddl = Dv.Taille();
   Tableau <double> tab_Nv(nbddl);
   if (norme < ConstMath::trespetit) {return tab_Nv;}; // retour dans le cas où la vitesse est nulle
   double unsurnorme = 1./norme;   
   for (int i=1;i<=nbddl;i++)
    { tab_Nv(i) = (V*Dv(i))*unsurnorme;};
   return tab_Nv;
  };

  // calcul de l'inverse d'une matrice 3x3 donnée par ces coordonnées, en retour la matrice
  // d'entrée est remplacée par la matrice inverse
  //  1) cas d'une matrice non symétrique (quelconque), rangement des valeurs:
  // (1,1) ; (1,2) ; (1,3) ; (2,1) ; (2,2) ; (2,3) ; (3,1) ; (3,2) ; (3,3)
#ifndef MISE_AU_POINT
  inline 
#endif
void Util::Inverse_mat3x3(listdouble9Iter & i9Iter)

  // le codage est très basique, pour essayer d'avoir un calcul très efficace !!
  
{ // on commence par créer une matrice qui contient une copie de la matrice a inverser
  double t[9];
  double * tex = (i9Iter)->donnees; // idem pour la matrice passée en paramètre
  for ( int i=0; i< 9; i++)
        t[i] = tex[i];     // recopie des valeurs 
  // pour plus de clarté on va définir des identifiants explicites
  double*  M11 =&t[0];
  double*  M12 =&t[1];
  double*  M13 =&t[2];
  double*  M21 =&t[3];
  double*  M22 =&t[4];
  double*  M23 =&t[5];
  double*  M31 =&t[6];
  double*  M32 =&t[7];
  double*  M33 =&t[8];
  // on calcul un vecteur d'équilibrage V
  double V[3];
  V[0] = DabsMaX(*M11,DabsMaX(*M12, *M13));
  V[1] = DabsMaX(*M21,DabsMaX(*M22, *M23));
  V[2] = DabsMaX(*M31,DabsMaX(*M32, *M33));
  if ( (V[0] < ConstMath::trespetit ) || (V[1] < ConstMath::trespetit ) || (V[2] < ConstMath::trespetit ))
    { cout << "\n erreur *** la matrice 3x3 est singuliere et ne peut pas etre inversee"
	        << "\n m11= "<< *M11 <<" m12= "<< *M12 <<" m13= "<< *M13
	        << "\n m21= "<< *M21 <<" m22= "<< *M22 <<" m23= "<< *M23
	        << "\n m31= "<< *M31 <<" m32= "<< *M32 <<" m33= "<< *M33 ;
		if (ParaGlob::NiveauImpression() > 5)
		  cout << "\n Util::Inverse_mat3x3(listdouble9Iter & i9Iter)";
		Sortie(1);  	  
     };
   V[0] = 1./V[0]; V[1] = 1./V[1];V[2] = 1./V[2];
	// def d'un vecteur de permutation
	short int indi[3]; // par exemple indi[1]+1 = indique pour la 2ième ligne de  la nouvelle matrice alph_beta,
	                  //  à quelle ligne de l'ancienne matrice cela correspond 
	short int jdi[3]; // par exemple jdi[1]+1 = indique pour la 2ième ligne de l'ancienne matrice, 
						  //  à quelle ligne de la nouvelle matrice alpha_beta cela correspond
						  // permet de faire le retour inverse de la permutation 
	
//-----------------------------------------------------------------------
// --- A triangulation de la matrice avec équilibrage suivant les lignes
//-----------------------------------------------------------------------
// Beta11  beta12  beta13
// alpha21 beta22  beta23
// alpha31 alpha32 beta33

// 1) == cas de la première nouvelle colonne	
   // on commence par chercher la ligne telle que |Mi1|*Vi est maxi
	double mvi1 = Dabs(*M11) * V[0];
	double mvi2 = Dabs(*M21) * V[1]; double mvi3 = Dabs(*M31) * V[2];
	
	if (mvi1 >= mvi2) { if (mvi1 >= mvi3) { indi[0]=0;jdi[0]=0;} // mvi1 est le plus grand
	                    else              { indi[0]=2;jdi[2]=0;}; // mvi3 est le plus grand
							}  
	else    { if (mvi2 >= mvi3) { indi[0]=1;jdi[1]=0;}  // mvi2 est le plus grand
	          else  { indi[0]=2;jdi[2]=0;};   // mvi3 est le plus grand 
			  };
   // on permutte les lignes 1 et indi[0]+1 si indi[0] n'est pas égal à 0
	double inter = 0.;
   switch (indi[0]) 
	{ case 1: // permutation des lignes 1 et 2
			{inter = *M11; *M11 = *M21; *M21 = inter;
			 inter = *M12; *M12 = *M22; *M22 = inter;
			 inter = *M13; *M13 = *M23; *M23 = inter;
			 inter = V[0]; V[0] = V[1]; V[1] = inter;
			 }
			break;
	  case 2: // permutation des lignes 1 et 3
			{inter = *M11; *M11 = *M31; *M31 = inter;
			 inter = *M12; *M12 = *M32; *M32 = inter;
			 inter = *M13; *M13 = *M33; *M33 = inter;
			 inter = V[0]; V[0] = V[2]; V[2] = inter;
			 }
			break;
	};	
   // calcule de la première colonne
	// beta11 = M11 donc rien n'a faire pour le premier terme
	*M21 /= *M11; // alpha21 = M21 / beta11
	*M31 /= *M11; // alpha31 = M31 / beta11
	
// 2) == cas de la deuxième nouvelle colonne	
   //Beta12 = M12 : donc rien n'a faire pour Beta12
	
	// choix entre la 2ième et 3ième ligne
	double p22 = *M22 - *M21 * *M12;
	double p32 = *M32 - *M31 * *M12;
	if(V[1]*Dabs(p22) < V[2]*Dabs(p32))
	 { // cas ou l'on inverse les lignes 2 et 3
	   indi[1] = 2;jdi[2]=1;indi[2] = 1;jdi[1]=2;
		inter = *M21; *M21 = *M31; *M31 = inter;
		inter = *M22; *M22 = *M32; *M32 = inter;
		inter = *M23; *M23 = *M33; *M33 = inter;
		inter = p22; p22 = p32; p32 = inter;
	 }
	else // sinon c'est ok
	 { indi[1] = 1;jdi[1]=1;indi[2] = 2;jdi[2]=2;}; 
	
	// calcul de la deuxième colonne
	// beta12 = M12 donc rien n'a faire pour ce premier terme	
	*M22 = p22; // beta22 = M22 - alpha21 * beta12
	*M32 = p32 / *M22; // alpha32 = 1/beta22 * (M32 - alpha31*beta12)

// 3) == cas de la troisième colonne
   // beta13 = M13 donc rien n'a faire avec ce terme
	*M23 -= *M21 * *M13 ; // beta23 = M23 - alpha21*beta13
	*M33 -= *M31 * *M13 - *M32 * *M23 ; // beta33 = M33 - alpha31*beta13 - alpha32*beta23
//-----------------------------------------------------------------------
// --- B calcul de l'inverse de la matrice en résolvant 3 problèmes
//-----------------------------------------------------------------------
	// on a les formules suivantes:
	// [alpha](Z)=(Y)   sachant que alphaii=1
	// -> Z1=Y1; Z2=Y2-alpha21*Z1; Z3=Y3-alpha31*Z1-alpha32*Z2
	//
	// ensuite [beta](X)=(Z) d'où
	// -> X3=Z3/beta33; X2=1/beta22 * (Z2-beta23*X3); 
	//    X1=1./beta11 * (Z1-beta12*X2-beta13*X3)
   //
	// il faut tenir compte des permutations éventuelles Yi -> Y_i
	// du coup on ne sait pas exactement quelles sont les grandeurs nulles ou pas 

	//1) première colonne
	// Y1=1.,Y2=0.,Y3=0.;  et sans permutation on a X1=tex[0], X2 = tex[3] et X3 = tex[6]
	double Y[3];double * pt[3];
	{Y[0]=1.;Y[1]=0.;Y[2]=0.; 
	 double Z1 = Y[indi[0]]; // Z1 = Y_1
	 double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]]; 
	 double Z2= Y_2- *M21 * Z1; double Z3=Y_3- *M31 * Z1- *M32 * Z2;
	 pt[0] = &tex[0]; pt[1] = &tex[3]; pt[2] = &tex[6];
	 *(pt[jdi[2]]) = Z3/ *M33; *(pt[jdi[1]]) = (Z2 - *M23 * (*(pt[jdi[2]])) )/  *M22;
	 *(pt[jdi[0]]) = (Z1 - *M12 * (*(pt[jdi[1]])) - *M13 * (*(pt[jdi[2]])) )/  *M11;
	}  
	//2) deuxième colonne
	// Y1=0.,Y2=1.,Y3=0.;  et sans permutation on a  X1=tex[1], X2 = tex[4] et X3 = tex[7]
	{Y[0]=0.;Y[1]=1.;Y[2]=0.; 
	 double Z1 = Y[indi[0]]; // Z1 = Y_1
	 double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]]; 
	 double Z2= Y_2- *M21 * Z1; double Z3=Y_3- *M31 * Z1- *M32 * Z2;
	 pt[0] = &tex[1]; pt[1] = &tex[4]; pt[2] = &tex[7];
	 *(pt[jdi[2]]) = Z3/ *M33; *(pt[jdi[1]]) = (Z2 - *M23 * (*(pt[jdi[2]])) )/ *M22;
	 *(pt[jdi[0]]) = (Z1 - *M12 * (*(pt[jdi[1]])) - *M13 * (*(pt[jdi[2]])) )/  *M11;
	}  
	//2) troisième colonne
	// Y1=0.,Y2=0.,Y3=1.;  et sans permutation on a  X1=tex[2], X2 = tex[5] et X3 = tex[8]
	{Y[0]=0.;Y[1]=0.;Y[2]=1.; 
	 double Z1 = Y[indi[0]]; // Z1 = Y_1
	 double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]]; 
	 double Z2= Y_2- *M21 * Z1; double Z3=Y_3- *M31 * Z1- *M32 * Z2;
	 pt[0] = &tex[2]; pt[1] = &tex[5]; pt[2] = &tex[8];
	 *(pt[jdi[2]]) = Z3/ *M33; *(pt[jdi[1]]) = (Z2 - *M23 * (*(pt[jdi[2]])) )/ *M22;
	 *(pt[jdi[0]]) = (Z1 - *M12 * (*(pt[jdi[1]])) - *M13 * (*(pt[jdi[2]])) )/ *M11;
	}  


};
 
  // calcul de l'inverse d'une matrice 3x3 donnée par ces coordonnées, en retour la matrice
  // d'entrée est remplacée par la matrice inverse
  //  2) cas d'une matrice symétrique , rangement des valeurs:
  // (1,1) ; (2,2) ; (3,3) ; (2,1) = (1,2) ; (3,2) = (2,3) ; (3,1) = (1,3) ; 
#ifndef MISE_AU_POINT
  inline 
#endif
void Util::Inverse_mat3x3(listdouble6Iter & i6Iter)
  // le codage est très basique, pour essayer d'avoir un calcul très efficace !!
  
{ // on commence par créer une matrice qui contient une copie de la matrice a inverser
  double t[9];
  double * tex = (i6Iter)->donnees; // idem pour la matrice passée en paramètre
  // recopie des valeurs
  t[0] = tex[0]; t[4] = tex[1]; t[8] = tex[2]; 
  t[1] = t[3] = tex[3]; t[5] = t[7]= tex[4]; t[6] = t[2] = tex[5]; 
  // pour plus de clarté on va définir des identifiants explicites
  double*  A11 = &t[0];
  double*  A12 = &t[1];
  double*  A13 = &t[2];
  double*  A21 = &t[3];
  double*  A22 = &t[4];
  double*  A23 = &t[5];
  double*  A31 = &t[6];
  double*  A32 = &t[7];
  double*  A33 = &t[8];
//cout << "\n " << *A11 <<" "<< *A12 <<" "<< *A13 <<" "<< *A21 <<" "<< *A22 <<" "<< *A23 <<" "<< *A31 <<" "<< *A32 <<" "<< *A33 <<endl;
  // on calcul un vecteur d'équilibrage V
  double V[3];
  V[0] = DabsMaX(*A11,DabsMaX(*A12, *A13));
  V[1] = DabsMaX(*A21,DabsMaX(*A22, *A23));
  V[2] = DabsMaX(*A31,DabsMaX(*A32, *A33));
  if ( (V[0] < ConstMath::trespetit ) || (V[1] < ConstMath::trespetit ) || (V[2] < ConstMath::trespetit ))
    { cout << "\n erreur *** la matrice 3x3 est singuliere et ne peut pas etre inversee"
	        << "\n m11= "<< *A11 <<" m12= "<< *A12 <<" m13= "<< *A13
	        << "\n m21= "<< *A21 <<" m22= "<< *A22 <<" m23= "<< *A23
	        << "\n m31= "<< *A31 <<" m32= "<< *A32 <<" m33= "<< *A33 ;
		if (ParaGlob::NiveauImpression() > 5)
		  cout << "\n Util::Inverse_mat3x3(listdouble9Iter & i9Iter)";
		Sortie(1);  	  
     };
   V[0] = 1./V[0]; V[1] = 1./V[1];V[2] = 1./V[2];
	// def d'un vecteur de permutation
	short int indi[3]; // par exemple indi[1]+1 = indique pour la 2ième ligne de  la nouvelle matrice alph_beta,
	                  //  à quelle ligne de l'ancienne matrice cela correspond 
	short int jdi[3]; // par exemple jdi[1]+1 = indique pour la 2ième ligne de l'ancienne matrice, 
						  //  à quelle ligne de la nouvelle matrice alpha_beta cela correspond
						  // permet de faire le retour inverse de la permutation 
	
//-----------------------------------------------------------------------
// --- A triangulation de la matrice avec équilibrage suivant les lignes
//-----------------------------------------------------------------------
// Beta11  beta12  beta13
// alpha21 beta22  beta23
// alpha31 alpha32 beta33

// 1) == cas de la première nouvelle colonne	
   // on commence par chercher la ligne telle que |Mi1|*Vi est maxi
	double mvi1 = Dabs(*A11) * V[0];
	double mvi2 = Dabs(*A21) * V[1]; double mvi3 = Dabs(*A31) * V[2];
	
	if (mvi1 >= mvi2) { if (mvi1 >= mvi3) { indi[0]=0;jdi[0]=0;} // mvi1 est le plus grand
	                    else              { indi[0]=2;jdi[2]=0;}; // mvi3 est le plus grand
							}  
	else   {if (mvi2 >= mvi3) { indi[0]=1;jdi[1]=0;}  // mvi2 est le plus grand
	        else { indi[0]=2;jdi[2]=0;};   // mvi3 est le plus grand
			 };  
   // on permutte les lignes 1 et indi[0] si indi[0] n'est pas égal à 1
	double inter = 0.;
   switch (indi[0]) 
	{ case 1: // permutation des lignes 1 et 2
			{inter = *A11; *A11 = *A21; *A21 = inter;
			 inter = *A12; *A12 = *A22; *A22 = inter;
			 inter = *A13; *A13 = *A23; *A23 = inter;
			 inter = V[0]; V[0] = V[1]; V[1] = inter;
			 }
			break;
	  case 2: // permutation des lignes 1 et 3
			{inter = *A11; *A11 = *A31; *A31 = inter;
			 inter = *A12; *A12 = *A32; *A32 = inter;
			 inter = *A13; *A13 = *A33; *A33 = inter;
			 inter = V[0]; V[0] = V[2]; V[2] = inter;
			 }
			break;
	};	
   // calcule de la première colonne
	// beta11 = A11 donc rien n'a faire pour le premier terme
	*A21 /= *A11; // alpha21 = A21 / beta11
	*A31 /= *A11; // alpha31 = A31 / beta11
	
// 2) == cas de la deuxième nouvelle colonne	
   //Beta12 = A12 : donc rien n'a faire pour Beta12
	
	// choix entre la 2ième et 3ième ligne
	double p22 = *A22 - *A21 * *A12;
	double p32 = *A32 - *A31 * *A12;
	if(V[1]*Dabs(p22) < V[2]*Dabs(p32))
	 { // cas ou l'on inverse les lignes 2 et 3
	   indi[1] = 2;jdi[2]=1;indi[2] = 1;jdi[1]=2;
		inter = *A21; *A21 = *A31; *A31 = inter;
		inter = *A22; *A22 = *A32; *A32 = inter;
		inter = *A23; *A23 = *A33; *A33 = inter;
		inter = p22; p22 = p32; p32 = inter;
	 }
	else // sinon c'est ok
	 { indi[1] = 1;jdi[1]=1;indi[2] = 2;jdi[2]=2;}; 
	
	// calcul de la deuxième colonne
	// beta12 = A12 donc rien n'a faire pour ce premier terme	
	*A22 = p22; // beta22 = A22 - alpha21 * beta12
	*A32 = p32 / *A22; // alpha32 = 1/beta22 * (A32 - alpha31*beta12)

// 3) == cas de la troisième colonne
   // beta13 = A13 donc rien n'a faire avec ce terme
	*A23 -= *A21 * *A13 ; // beta23 = A23 - alpha21*beta13
	*A33 -= *A31 * *A13 - *A32 * *A23 ; // beta33 = A33 - alpha31*beta13 - alpha32*beta23
//-----------------------------------------------------------------------
// --- B calcul de l'inverse de la matrice en résolvant 3 problèmes
//-----------------------------------------------------------------------
	// on a les formules suivantes:
	// [alpha](Z)=(Y)   sachant que alphaii=1
	// -> Z1=Y1; Z2=Y2-alpha21*Z1; Z3=Y3-alpha31*Z1-alpha32*Z2
	//
	// ensuite [beta](X)=(Z) d'où
	// -> X3=Z3/beta33; X2=1/beta22 * (Z2-beta23*X3); 
	//    X1=1./beta11 * (Z1-beta12*X2-beta13*X3)
   //
	// il faut tenir compte des permutations éventuelles Yi -> Y_i
	// du coup on ne sait pas exactement quelles sont les grandeurs nulles ou pas 


	//1) première colonne
	// Y1=1.,Y2=0.,Y3=0.;  et sans permutation on a  X1=tex[0], X2 = tex[3] et X3 = tex[5]
	double Y[3];double * pt[3];
	{Y[0]=1.;Y[1]=0.;Y[2]=0.; 
	 double Z1 = Y[indi[0]]; // Z1 = Y_1
	 double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]]; 
	 double Z2= Y_2-  *A21 * Z1; double Z3=Y_3-  *A31 * Z1- *A32 * Z2;
	 pt[0] = &tex[0]; pt[1] = &tex[3]; pt[2] = &tex[5];
	 *(pt[jdi[2]]) = Z3/ *A33; *(pt[jdi[1]]) = (Z2 - *A23 * (*(pt[jdi[2]])) )/  *A22;
	 *(pt[jdi[0]]) = (Z1 - *A12 * (*(pt[jdi[1]])) - *A13 * (*(pt[jdi[2]])) )/  *A11;
	}  
	//2) deuxième colonne
	// Y1=0.,Y2=1.,Y3=0.;  et sans permutation on a X1=tex[3] déjà calculé, X2 = tex[1] et X3 = tex[4]
	// on pourrait sans doute faire des éconnomies due à la symétrie, cependant due aux permutations
	// ce n'est pas évident d'utiliser correctement les symétries, donc on recalcul tout
	{Y[0]=0.;Y[1]=1.;Y[2]=0.; 
	 double Z1 = Y[indi[0]]; // Z1 = Y_1
	 double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]]; 
	 double Z2= Y_2-  *A21 * Z1; double Z3=Y_3-  *A31 * Z1- *A32 * Z2;
	 pt[0] = &tex[3]; pt[1] = &tex[1]; pt[2] = &tex[4];
	 *(pt[jdi[2]]) = Z3/  *A33; *(pt[jdi[1]]) = (Z2 - *A23 * (*(pt[jdi[2]])) )/  *A22;
	 *(pt[jdi[0]]) = (Z1 - *A12 * (*(pt[jdi[1]])) - *A13 * (*(pt[jdi[2]])) )/  *A11;
	}  
	//2) troisième colonne
	// Y1=0.,Y2=0.,Y3=1.;  et sans permutation on a  X1=tex[5] déjà calculé, X2 = tex[4] déjà calculé et X3 = tex[2]
	{Y[0]=0.;Y[1]=0.;Y[2]=1.; 
	 double Z1 = Y[indi[0]]; // Z1 = Y_1
	 double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]]; 
	 double Z2= Y_2-  *A21 * Z1; double Z3=Y_3-  *A31 * Z1- *A32 * Z2;
	 pt[0] = &tex[5]; pt[1] = &tex[4]; pt[2] = &tex[2];
	 *(pt[jdi[2]]) = Z3/  *A33; *(pt[jdi[1]]) = (Z2 - *A23 * (*(pt[jdi[2]])) )/  *A22;
	 *(pt[jdi[0]]) = (Z1 - *A12 * (*(pt[jdi[1]])) - *A13 * (*(pt[jdi[2]])) )/  *A11;
	}  


};
  // calcul de l'inverse d'une matrice 2x2 donnée par ces coordonnées, en retour la matrice
  // d'entrée est remplacée par la matrice inverse
  //  1) cas d'une matrice non symétrique (quelconque), rangement des valeurs:
  // (1,1) ; (2,2) ; (2,1) ; (1,2) ; 
#ifndef MISE_AU_POINT
  inline 
#endif
void Util::Inverse_mat2x2(listdouble4Iter & i4Iter)

  // le codage est très basique, pour essayer d'avoir un calcul très efficace !!
  
{ // on commence par créer une matrice qui contient une copie de la matrice a inverser
  double t[4];
  double * tex = (i4Iter)->donnees; // idem pour la matrice passée en paramètre
  // recopie des valeurs
  t[0] = tex[0]; t[3] = tex[1]; t[2] = tex[2]; t[1] = tex[4];
  // pour plus de clarté on va définir des identifiants explicites
// *** finalement on n'utilise pas de define car cela donne une grandeur définie pour tous Herezh
// et dans certaine méthodes on utilise des variables qui ont ce nom là !!
// #define B11 t[0]
// #define B12 t[1]
// #define B21 t[2]
// #define B22 t[3]
 double*  B11 = &t[0];
 double*  B12 = &t[1];
 double*  B21 = &t[2];
 double*  B22 = &t[3];

//cout << "\n " << *B11 <<" "<< *B12 <<" "<< *B21 <<" "<< *B22 << endl;
  // on calcul un vecteur d'équilibrage V
  double V[2];
  V[0] = DabsMaX(*B11,*B12);
  V[1] = DabsMaX(*B21,*B22);
  if ( (V[0] < ConstMath::trespetit ) || (V[1] < ConstMath::trespetit ) )
    { cout << "\n erreur *** la matrice 2x2 est singuliere et ne peut pas etre inversee"
	        << "\n m11= "<< *B11 <<" m12= "<< *B12
	        << "\n m21= "<< *B21 <<" m22= "<< *B22 ;
		if (ParaGlob::NiveauImpression() > 5)
		  cout << "\n Util::Inverse_mat2x2(listdouble4Iter & i4Iter)";
		Sortie(1);  	  
     };
   V[0] = 1./V[0]; V[1] = 1./V[1];
	// def d'un indicateur de permutation
	bool permutation = false; // indique si oui ou non on permutte la première ligne avec la seconde
	
// 1) == cas de la première nouvelle colonne	
   // on commence par chercher la ligne telle que |Mi1|*Vi est maxi
	double mvi1 = Dabs(*B11) * V[0];
	double mvi2 = Dabs(*B21) * V[1];
	
	double inter = 0.;
	if (mvi1 < mvi2) 
	  { permutation = true; // mvi2 est le plus grand
       // permutation des lignes 1 et 2
		 inter = *B11; *B11 = *B21; *B21 = inter;
		 inter = *B12; *B12 = *B22; *B22 = inter;
		 inter = V[0]; V[0] = V[1]; V[1] = inter;
		};

   // calcule de la première colonne
	// beta11 = A11 donc rien n'a faire pour le premier terme
	*B21 /= *B11; // alpha21 = B21 / beta11
	
// 2) == cas de la deuxième nouvelle colonne	
   //Beta12 = B12 : donc rien n'a faire pour Beta12
	// calcul du deuxième terme de la deuxième colonne
	*B22 = *B22 - *B21 * *B12; // beta22 = B22 - alpha21 * beta12

//-----------------------------------------------------------------------
// --- B calcul de l'inverse de la matrice en résolvant 2 problèmes
//-----------------------------------------------------------------------
	// on a les formules suivantes:
	// [alpha](Z)=(Y)   sachant que alphaii=1
	// -> Z1=Y1; Z2=Y2-alpha21*Z1; 
	//
	// ensuite [beta](X)=(Z) d'où
	//    X2=Z2/beta22 ; 
	//    X1=1./beta11 * (Z1-beta12*X2)
   //
	// il faut tenir compte d'une permutation éventuelle 
	
	if (!permutation)
	 { //1) première colonne
	   // Y1=1.,Y2=0.,  et sans permutation on a  X1=tex[0], X2 = tex[2] 
	   double Y2=0.;
	   double Z1 = 1.;
    double Z2 = Y2 - *B21 * Z1;
    tex[2] = Z2 / *B22; tex[0] = (Z1 - *B12 * tex[2] )/ *B11;
      //2) deuxième colonne
      // Y1=0.,Y2=1.;  et sans permutation on a X1=tex[3], X2 = tex[1]
      // on pourrait faire des éconnomies due à la symétrie, cependant due aux permutations on recalcul tout
    Y2=1.;
    Z1 = 0.; // Z1 = Y_1
	   Z2 = Y2- *B21 * Z1;
		  tex[1] = Z2 / *B22; tex[3] = (Z1 - *B12 * tex[1] )/ *B11;
	  }
	else 
	 { //avec permutation: 1) première colonne
	   // Y1=0.,Y2=1.,  et avec permutation on a  X1=tex[2], X2 = tex[0] 
	   double Y2=1.;
	   double Z1 = 0.;
    double Z2 = Y2 - *B21 * Z1;
    tex[0] = Z2 / *B22; tex[2] = (Z1 - *B12 * tex[0] )/ *B11;
	   //2) deuxième colonne
	   // Y1=1.,Y2=0.;  et sans permutation on a X1=tex[1] , X2 = tex[3] 
	   // on pourrait faire des éconnomies due à la symétrie, cependant due aux permutations on recalcul tout
    Y2=0.;
    Z1 = 1.; // Z1 = Y_1
	   Z2 = Y2- *B21 * Z1;
  		tex[3] = Z2 / *B22; tex[1] = (Z1 - *B12 * tex[3] )/ *B11;
	  };

};
 
  // calcul de l'inverse d'une matrice 2x2 donnée par ces coordonnées, en retour la matrice
  // d'entrée est remplacée par la matrice inverse
  //  2) cas d'une matrice symétrique , rangement des valeurs:
  // (1,1) ; (2,2) ; (2,1) = (1,2) ;
#ifndef MISE_AU_POINT
  inline 
#endif
void Util::Inverse_mat2x2(listdouble3Iter & i3Iter)
  // le codage est très basique, pour essayer d'avoir un calcul très efficace !!
  
{ // on commence par créer une matrice qui contient une copie de la matrice a inverser
  double t[4];
  double * tex = (i3Iter)->donnees; // idem pour la matrice passée en paramètre
  // recopie des valeurs
  t[0] = tex[0]; t[3] = tex[1]; t[1] = t[2] = tex[2]; 
  // pour plus de clarté on va définir des identifiants explicites
  double*  C11 = &t[0];
  double*  C12 = &t[1];
  double*  C21 = &t[2];
  double*  C22 = &t[3];
  //cout << "\n " << *C11 <<" "<< *C12 <<" "<< *C21 <<" "<< *C22 << endl;
  // on calcul un vecteur d'équilibrage V
  double V[2];
  V[0] = DabsMaX(*C11,*C12);
  V[1] = DabsMaX(*C21,*C22);
  if ( (V[0] < ConstMath::trespetit ) || (V[1] < ConstMath::trespetit ) )
    { cout << "\n erreur *** la matrice 2x2 est singuliere et ne peut pas etre inversee"
	        << "\n m11= "<< *C11 <<" m12= "<< *C12
	        << "\n m21= "<< *C21 <<" m22= "<< *C22 ;
		if (ParaGlob::NiveauImpression() > 5)
		  cout << "\n Util::Inverse_mat2x2(listdouble3Iter & i3Iter)";
		Sortie(1);  	  
     };
   V[0] = 1./V[0]; V[1] = 1./V[1];
	// def d'un indicateur de permutation
	bool permutation = false; // indique si oui ou non on permutte la première ligne avec la seconde
	
//-----------------------------------------------------------------------
// --- A triangulation de la matrice avec équilibrage suivant les lignes
//-----------------------------------------------------------------------
// Beta11  beta12  
// alpha21 beta22  

// 1) == cas de la première nouvelle colonne	
   // on commence par chercher la ligne telle que |Mi1|*Vi est maxi
	double mvi1 = Dabs(*C11) * V[0];
	double mvi2 = Dabs(*C21) * V[1];
	
	double inter = 0.;
	if (mvi1 < mvi2) 
	  { permutation = true; // mvi2 est le plus grand
       // permutation des lignes 1 et 2
		 inter = *C11; *C11 = *C21; *C21 = inter;
		 inter = *C12; *C12 = *C22; *C22 = inter;
		 inter = V[0]; V[0] = V[1]; V[1] = inter;
		};
   // calcule de la première colonne
	// beta11 = A11 donc rien n'a faire pour le premier terme
	*C21 /= *C11; // alpha21 = C21 / beta11
	
// 2) == cas de la deuxième nouvelle colonne	
   //Beta12 = C12 : donc rien n'a faire pour Beta12
	// calcul du deuxième terme de la deuxième colonne
	*C22 = *C22 - *C21 * *C12; // beta22 = C22 - alpha21 * beta12

//-----------------------------------------------------------------------
// --- B calcul de l'inverse de la matrice en résolvant 2 problèmes
//-----------------------------------------------------------------------
	// on a les formules suivantes:
	// [alpha](Z)=(Y)   sachant que alphaii=1
	// -> Z1=Y1; Z2=Y2-alpha21*Z1; 
	//
	// ensuite [beta](X)=(Z) d'où
	//    X2=Z2/beta22 ; 
	//    X1=1./beta11 * (Z1-beta12*X2)
   //
	// il faut tenir compte d'une permutation éventuelle 
	
	if (!permutation)
	 { //1) première colonne
	   // Y1=1.,Y2=0.,  et sans permutation on a  X1=tex[0], X2 = tex[2] 
	   double Y2=0.;
	   double Z1 = 1.;
		double Z2 = Y2 - *C21 * Z1;
		tex[2] = Z2 / *C22; tex[0] = (Z1 - *C12 * tex[2] )/ *C11;
	   //2) deuxième colonne
	   // Y1=0.,Y2=1.;  et sans permutation on a X1=tex[2] déjà calculé, X2 = tex[1] 
	   // on pourrait faire des éconnomies due à la symétrie, cependant due aux permutations on recalcul tout
		Y2=1.; 
		Z1 = 0.; // Z1 = Y_1
	   Z2 = Y2-  *C21 * Z1;
		tex[1] = Z2 / *C22; tex[2] = (Z1 - *C12 * tex[1] )/  *C11;
	  }
	else 
	 { //avec permutation: 1) première colonne
	   // Y1=0.,Y2=1.,  et avec permutation on a  X1=tex[2], X2 = tex[0] 
	   double Y2=1.;
	   double Z1 = 0.;
		double Z2 = Y2 - *C21 * Z1;
		tex[0] = Z2 / *C22; tex[2] = (Z1 - *C12 * tex[0] )/  *C11;
	   //2) deuxième colonne
	   // Y1=1.,Y2=0.;  et sans permutation on a X1=tex[1] déjà calculé, X2 = tex[2] 
	   // on pourrait faire des éconnomies due à la symétrie, cependant due aux permutations on recalcul tout
		Y2=0.; 
		Z1 = 1.; // Z1 = Y_1
	   Z2 = Y2-  *C21 * Z1;
		tex[2] = Z2 / *C22; tex[1] = (Z1 - *C12 * tex[2] )/  *C11;
	  };

};



#endif
