outil/src/ot_mathematique.cpp

//------------------------------------------------------------
//------------------------------------------------------------
//  MAGiC
//  Jean Christophe Cuillière et Vincent FRANCOIS
//  Département de Génie Mécanique - UQTR
//------------------------------------------------------------
//  Le projet MAGIC est un projet de recherche du département
//  de génie mécanique de l'Université du Québec à
//  Trois Rivières
//  Les librairies ne peuvent être utilisées sans l'accord
//  des auteurs (contact : francois@uqtr.ca)
//------------------------------------------------------------
//------------------------------------------------------------
//
//       ot_mathematique.cpp
//
//------------------------------------------------------------
//------------------------------------------------------------
//  COPYRIGHT 2000
//  Version du 02/03/2006 à 11H23
//------------------------------------------------------------
//------------------------------------------------------------


#include "gestionversion.h"

#include <math.h>
#include <algorithm>
#include "ot_mathematique.h"

using namespace std;


OT_VECTEUR_3D::OT_VECTEUR_3D(double x,double y,double z)
{
valeur[0]=x;
valeur[1]=y;
valeur[2]=z;
}


OT_VECTEUR_3D::OT_VECTEUR_3D(double *xyz)
{
valeur[0]=xyz[0];
valeur[1]=xyz[1];
valeur[2]=xyz[2];
}

OT_VECTEUR_3D::OT_VECTEUR_3D(double *xyz1,double *xyz2)
{
valeur[0]=xyz2[0]-xyz1[0];
valeur[1]=xyz2[1]-xyz1[1];
valeur[2]=xyz2[2]-xyz1[2];
}

OT_VECTEUR_3D::OT_VECTEUR_3D(const OT_VECTEUR_3D& mdd)
{
valeur[0]=mdd.valeur[0];
valeur[1]=mdd.valeur[1];
valeur[2]=mdd.valeur[2];
}

OT_VECTEUR_3D::OT_VECTEUR_3D(OT_VECTEUR_3D& mdd)
{
valeur[0]=mdd.valeur[0];
valeur[1]=mdd.valeur[1];
valeur[2]=mdd.valeur[2];
}


double OT_VECTEUR_3D::get_x(void) const 
{
return valeur[0];
}


double OT_VECTEUR_3D::get_y(void) const 
{
return valeur[1];
}


double OT_VECTEUR_3D::get_z(void) const
{
return valeur[2];
}

void OT_VECTEUR_3D::change_x(double x)
{
valeur[0]=x;
}

void OT_VECTEUR_3D::change_y(double y)
{
valeur[1]=y;
}

void OT_VECTEUR_3D::change_z(double z)
{
valeur[2]=z;
}

double* OT_VECTEUR_3D::get_xyz(void)
{
return valeur;
}
                       
OT_VECTEUR_3D::operator const double* () const
{
    return valeur;
}

OT_VECTEUR_3D::operator double* ()
{
    return valeur;
}

double OT_VECTEUR_3D::operator[] (int i) const
{
    return valeur[i];
}

double & OT_VECTEUR_3D::operator[] (int i)
{
    return valeur[i];
}

double OT_VECTEUR_3D::operator() (int i) const
{
    return valeur[i];
}
double & OT_VECTEUR_3D::operator() (int i) 
{
    return valeur[i];
}

double OT_VECTEUR_3D::get_longueur(void) const 
{
    return sqrt(valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
}

double  OT_VECTEUR_3D::get_longueur2(void) const 
{
return (valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
}


void OT_VECTEUR_3D::norme(void)
{
double longueur=get_longueur();
if (longueur != 0)
{
    valeur[0]=valeur[0]/longueur;
    valeur[2]=valeur[2]/longueur;
    valeur[1]=valeur[1]/longueur;
}
}
                                                              
int  OT_VECTEUR_3D::compare_valeur (const OT_VECTEUR_3D& mdd) const
{
        return memcmp(valeur, mdd.valeur, 3*sizeof(double));
}

bool  OT_VECTEUR_3D::operator== (const OT_VECTEUR_3D& mdd) const
{
        return compare_valeur(mdd) == 0;
}

bool  OT_VECTEUR_3D::operator!= (const OT_VECTEUR_3D& mdd) const  
{                               
        return compare_valeur(mdd) != 0;
}

bool  OT_VECTEUR_3D::operator<  (const OT_VECTEUR_3D& mdd) const
{                                   
        return compare_valeur(mdd) < 0;
}

bool  OT_VECTEUR_3D::operator<= (const OT_VECTEUR_3D& mdd) const    
{                                
        return compare_valeur(mdd) <= 0;        
}

bool  OT_VECTEUR_3D::operator>  (const OT_VECTEUR_3D& mdd) const   
{                                     
        return compare_valeur(mdd) > 0;        
}

bool  OT_VECTEUR_3D::operator>= (const OT_VECTEUR_3D& mdd) const  
{
        return compare_valeur(mdd) >= 0;
}
                             
std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_3D & __vec)
{
        __os << "{ "<< __vec[0]<<", "<<__vec[1]<<", "<<__vec[2]<<"}, ";
        return __os;
}

OT_VECTEUR_3D operator/ (const OT_VECTEUR_3D& mdd1, double diviseur)
{
    OT_VECTEUR_3D kQuot;

    if ( diviseur != (double)0.0 )
    {
        double facteur = ((double)1.0)/diviseur;
        kQuot.valeur[0] = facteur*mdd1.valeur[0];
        kQuot.valeur[1] = facteur*mdd1.valeur[1];
        kQuot.valeur[2] = facteur*mdd1.valeur[2];
    }
    else
    {
        kQuot.valeur[0] = 1E300;
        kQuot.valeur[1] = 1E300;
        kQuot.valeur[2] = 1E300;
    }

    return kQuot;
}

OT_VECTEUR_3D operator- (const OT_VECTEUR_3D& mdd1) 
{
    return OT_VECTEUR_3D(
        -mdd1.valeur[0],
        -mdd1.valeur[1],
        -mdd1.valeur[2]);
}

double operator*(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
{
double tmp;
tmp=mdd1.valeur[0]*mdd2.valeur[0]+mdd1.valeur[1]*mdd2.valeur[1]+mdd1.valeur[2]*mdd2.valeur[2];
return tmp;
}

OT_VECTEUR_3D operator*(double mdd1,const OT_VECTEUR_3D& mdd2)
{
OT_VECTEUR_3D tmp;
tmp.valeur[0]=mdd1*mdd2.valeur[0];
tmp.valeur[1]=mdd1*mdd2.valeur[1];
tmp.valeur[2]=mdd1*mdd2.valeur[2];
return tmp;
}

OT_VECTEUR_3D& OT_VECTEUR_3D::operator+= (const OT_VECTEUR_3D& mdd)
{
    valeur[0] += mdd.valeur[0];
    valeur[1] += mdd.valeur[1];
    valeur[2] += mdd.valeur[2];
    return *this;
}
OT_VECTEUR_3D& OT_VECTEUR_3D::operator-= (const OT_VECTEUR_3D& mdd)
{
    valeur[0] -= mdd.valeur[0];
    valeur[1] -= mdd.valeur[1];
    valeur[2] -= mdd.valeur[2];
    return *this;
}                    
OT_VECTEUR_3D& OT_VECTEUR_3D::operator*= (double facteur)
{
    valeur[0] *= facteur;
    valeur[1] *= facteur;
    valeur[2] *= facteur;
    return *this;
}
OT_VECTEUR_3D& OT_VECTEUR_3D::operator/= (double diviseur)
{
    if ( diviseur != (double)0.0 )
    {
        double facteur = ((double)1.0)/diviseur;
        valeur[0] *= facteur;
        valeur[1] *= facteur;
        valeur[2] *= facteur;
    }
    else
    {
        valeur[0] = 1E300;
        valeur[1] = 1E300;
        valeur[2] = 1E300;
    }

    return *this;
}

OT_VECTEUR_3D operator&(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
{
OT_VECTEUR_3D tmp;
tmp.valeur[0]=mdd1.valeur[1]*mdd2.valeur[2]-mdd1.valeur[2]*mdd2.valeur[1];
tmp.valeur[1]=mdd1.valeur[2]*mdd2.valeur[0]-mdd1.valeur[0]*mdd2.valeur[2];
tmp.valeur[2]=mdd1.valeur[0]*mdd2.valeur[1]-mdd1.valeur[1]*mdd2.valeur[0];
return tmp;
}

double OT_VECTEUR_3D::diff(void)
{
return 0.33333333333333*(fabs(valeur[0])+fabs(valeur[1])+fabs(valeur[2]));
}


OT_MATRICE_3D::OT_MATRICE_3D(OT_VECTEUR_3D& colonne1,OT_VECTEUR_3D& colonne2,OT_VECTEUR_3D& colonne3)
{
vec[0]=colonne1;  
vec[1]=colonne2;
vec[2]=colonne3;
}

OT_MATRICE_3D::OT_MATRICE_3D(const OT_MATRICE_3D& mdd)
{
vec[0]=mdd.vec[0];
vec[1]=mdd.vec[1];
vec[2]=mdd.vec[2];
}


std::ostream & operator << ( std::ostream & __os, const OT_MATRICE_3D & __mat)
{

        __os << "{ ";
        for (unsigned i=0; i<3; i++)
        {
                __os << "[ "<< __mat(i,0)<<" ]\t [ "<<__mat(i,1) <<" ]\t [ "<<__mat(i,2);
                if (i+1 < 3)
                        __os << std::endl << "  ";
        }
        __os << " }" << std::endl;
        return __os;
}

OT_VECTEUR_3D & OT_MATRICE_3D::operator [](int i)
{
return vec[i];
}
                           
double OT_MATRICE_3D::valeur(int iLigne, int iCol) const
{
    return vec[iCol][iLigne];
}
double & OT_MATRICE_3D::valeur(int iLigne, int iCol)
{
    return vec[iCol][iLigne];
}
double OT_MATRICE_3D::operator() (int iLigne, int iCol) const
{
    return vec[iCol][iLigne];
}
double & OT_MATRICE_3D::operator() (int iLigne, int iCol)
{
    return vec[iCol][iLigne];
}
double OT_MATRICE_3D::get_determinant()
{
    double coFacteur0 = valeur(1,1)*valeur(2,2) - valeur(1,2)*valeur(2,1);
    double coFacteur1 = valeur(1,2)*valeur(2,0) - valeur(1,0)*valeur(2,2);
    double coFacteur2 = valeur(1,0)*valeur(2,1) - valeur(1,1)*valeur(2,0);
    double determinant = valeur(0,0)*coFacteur0 + valeur(0,1)*coFacteur1 + valeur(0,2)*coFacteur2;
    return determinant;
}
void OT_MATRICE_3D::identite(void)
{
        for (unsigned i = 0; i<3; i++)
                for (unsigned j = 0; j < 3; j++)
                        vec[i][j] = (i==j) ? 1.0 : 0.0 ;
}

void OT_MATRICE_3D::transpose(OT_MATRICE_3D& res) const
{
        for (unsigned i = 0; i<3; i++)
                for (unsigned j = 0; j < 3; j++)
                        res.vec[i][j] = vec[j][i];
}

OT_MATRICE_3D OT_MATRICE_3D::inverse () const
{
    // Invert a 3x3 using cofactors.  This is faster than using a generic
    // Gaussian elimination because of the loop overhead of such a method.

    OT_MATRICE_3D matInverse;

    matInverse.valeur(0,0) = valeur(1,1)*valeur(2,2) - valeur(1,2)*valeur(2,1);
    matInverse.valeur(0,1) = valeur(0,2)*valeur(2,1) - valeur(0,1)*valeur(2,2);
    matInverse.valeur(0,2) = valeur(0,1)*valeur(1,2) - valeur(0,2)*valeur(1,1);
    matInverse.valeur(1,0) = valeur(1,2)*valeur(2,0) - valeur(1,0)*valeur(2,2);
    matInverse.valeur(1,1) = valeur(0,0)*valeur(2,2) - valeur(0,2)*valeur(2,0);
    matInverse.valeur(1,2) = valeur(0,2)*valeur(1,0) - valeur(0,0)*valeur(1,2);
    matInverse.valeur(2,0) = valeur(1,0)*valeur(2,1) - valeur(1,1)*valeur(2,0);
    matInverse.valeur(2,1) = valeur(0,1)*valeur(2,0) - valeur(0,0)*valeur(2,1);
    matInverse.valeur(2,2) = valeur(0,0)*valeur(1,1) - valeur(0,1)*valeur(1,0);

    double det = valeur(0,0)*matInverse(0,0) + valeur(0,1)*matInverse(1,0)+
        valeur(0,2)*matInverse(2,0);

    if ( fabs(det) <= 1E-100 )
    {
        matInverse.identite();
        return matInverse;
    }

    double invDet = 1/det;

    for (int i=0; i<3; i++)
    for (int j=0; j<3; j++)
        matInverse.vec[i][j] *= invDet;

    return matInverse;
}

void OT_MATRICE_3D::change_vecteur1(OT_VECTEUR_3D v)
{
vec[0]=v;
}
void OT_MATRICE_3D::change_vecteur2(OT_VECTEUR_3D v)
{
vec[1]=v;
}
void OT_MATRICE_3D::change_vecteur3(OT_VECTEUR_3D v)
{
vec[2]=v;
}
OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur1(void)
{
return vec[0];
}
OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur2(void)
{
return vec[1];
}
OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur3(void)
{
return vec[2];
}


OT_MATRICE_3D OT_MATRICE_3D::VecteursPropres(OT_VECTEUR_3D& Val_Propre,double Eps,int ItMax) 
{
OT_MATRICE_3D m=*this;
OT_VECTEUR_3D ZERO(0.,0.,0.) ;
OT_MATRICE_3D vct_prpres(ZERO,ZERO,ZERO);
OT_VECTEUR_3D V ;

int Dim=3;
int I,J,K;
OT_VECTEUR_3D U, Up;
double N,T;
//-----
for(int vi=0;vi<3;vi++)
 {
//=================== 
Up=vec[0]; // il faut l'initialiser d'une autre maniere
  K=0;
  do{
    U=Up;
    for(int I=0 ;I<Dim;I++){
      T=0;
      for(int J=0; J<Dim;J++) T=T+valeur(I,J)*U[J];
      Up[I]=T  ;
     }
    N=Up.get_longueur();
    if( Up[0]<0)  N=-N;
    T=0;
    for(int I=0;I<Dim;I++){
      Up[I]=Up[I]/N;
      T=T+fabs(Up[I]-U[I]);
    }
    K=K+1;
    }
  while((K<ItMax)||(T>Eps)) ;
  V=Up  ;
  vct_prpres[vi]=V;
  //================
 double N,D;
  N=0.0; D=0.0;
  for(int I=0;I<Dim;I++){
  
    D=D+V[I]*V[I];
    for(int J=0 ;J<Dim;J++) N=N+V[I]*valeur(I,J)*V[J];
    }
  double Val=N/D;
  Val_Propre[vi]=Val; 
  //===============
    for(int I=0; I<Dim ;I++){

    for(int J=0;J<Dim;J++) m(I,J)=m(I,J)-Val*V[I]*V[J];
  }
  
  //===============
   
  }
  //-----------------

  
  return vct_prpres;
}

  /*
double OT_MATRICE_3D::ValeurPropre(double Eps,int ItMax)
{
double Val_Propre;
OT_VECTEUR_3D V=VecteurPropre( Eps, ItMax);
int I,J,Dim;
double N,D;
  N=0.0; D=0.0; Dim=3;
  for(int I=0;I<Dim;I++){
  
    D=D+V[I]*V[I];
    for(int J=0 ;J<Dim;J++) N=N+V[I]*valeur(I,J)*V[J];
    }
  Val_Propre=N/D;
return Val_Propre;
}
   
void OT_MATRICE_3D:: Deflation(OT_VECTEUR_3D& V,double val_prpre)
{
int Dim=3;
  for(int I=0; I<Dim ;I++){

    for(int J=0;J<Dim;J++) valeur(I,J)=valeur(I,J)-val_prpre*V[I]*V[J];
  }
  
}   */

 
OT_VECTEUR_3D operator*(const OT_MATRICE_3D& mdd1,const OT_VECTEUR_3D& mdd2)
{
OT_VECTEUR_3D tmp;
tmp(0)=mdd1(0,0)*mdd2(0)+mdd1(0,1)*mdd2(1)+mdd1(0,2)*mdd2(2);
tmp(1)=mdd1(1,0)*mdd2(0)+mdd1(1,1)*mdd2(1)+mdd1(1,2)*mdd2(2);
tmp(2)=mdd1(2,0)*mdd2(0)+mdd1(2,1)*mdd2(1)+mdd1(2,2)*mdd2(2);
return tmp;}

OT_VECTEUR_3D operator+(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
{
OT_VECTEUR_3D tmp;
tmp(0)=mdd1(0)+mdd2(0);
tmp(1)=mdd1(1)+mdd2(1);
tmp(2)=mdd1(2)+mdd2(2);
return tmp;
}

OT_VECTEUR_3D operator-(const class OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
{
OT_VECTEUR_3D tmp;
tmp.valeur[0]=mdd1.valeur[0]-mdd2.valeur[0];
tmp.valeur[1]=mdd1.valeur[1]-mdd2.valeur[1];
tmp.valeur[2]=mdd1.valeur[2]-mdd2.valeur[2];
return tmp;
}
                 
OT_MATRICE_3D operator*(const OT_MATRICE_3D& mdd1,const OT_MATRICE_3D& mdd2)
{
    OT_MATRICE_3D tmp;
    for (int i=0; i<3; i++)
    for (int j=0; j<3; j++)
    {
    tmp(i,j) = 0;
    for (int k=0; k<3; k++)
        tmp(i,j) += mdd1(i,k)*mdd2(k,j);
    }
    return tmp;
}

#ifdef ot_vecteur_4d
OT_VECTEUR_4D::OT_VECTEUR_4D (double v[4])
{
    for (int i=0; i<4; i++)
        valeur[i]=v[i];
}
                      
OT_VECTEUR_4D::OT_VECTEUR_4D (double x, double y, double z, double w)
{
    valeur[0] = x;
    valeur[1] = y;
    valeur[2] = z;
    valeur[3] = w;
}

OT_VECTEUR_4D::OT_VECTEUR_4D (const OT_VECTEUR_4D & v)
{
    for (int i=0; i<4; i++)
        valeur[i]=v[i];
}


double & OT_VECTEUR_4D::operator()(int i)
{
    return valeur[i];
}

double OT_VECTEUR_4D::operator()(int i) const
{
    return valeur[i];    
}

double & OT_VECTEUR_4D::operator[](int i)
{
    return valeur[i];
}

double OT_VECTEUR_4D::operator[](int i) const
{
    return valeur[i];    
}
               
OT_VECTEUR_4D::operator const double* () const
{
    return valeur;
}

OT_VECTEUR_4D::operator double* ()
{
    return valeur;
}

double OT_VECTEUR_4D::get_x(void) const
{ return valeur[0]; }              
double & OT_VECTEUR_4D::x(void)
{ return valeur[0]; }
double OT_VECTEUR_4D::get_y(void) const
{ return valeur[1]; }
double & OT_VECTEUR_4D::y(void)
{ return valeur[1]; }
double OT_VECTEUR_4D::get_z(void) const
{ return valeur[2]; }
double & OT_VECTEUR_4D::z(void)
{ return valeur[2]; }
double OT_VECTEUR_4D::get_w(void) const
{ return valeur[3]; }
double & OT_VECTEUR_4D::w(void) 
{ return valeur[3]; }
void OT_VECTEUR_4D::change_x(double x)
{ valeur[0] = x; }
void OT_VECTEUR_4D::change_y(double y) 
{ valeur[1] = y; }
void OT_VECTEUR_4D::change_z(double z)
{ valeur[2] = z; }
void OT_VECTEUR_4D::change_w(double w)
{ valeur[3] = w; }
double* OT_VECTEUR_4D::get_xyzw(void)
{ return valeur; }

OT_VECTEUR_4D OT_VECTEUR_4D::operator+ (const OT_VECTEUR_4D & a)
{
    return OT_VECTEUR_4D(
        a(0)+valeur[0],a(1)+valeur[1],a(2)+valeur[2],a(3)+valeur[3]
    );
}
     
OT_VECTEUR_4D OT_VECTEUR_4D::operator- (const OT_VECTEUR_4D & a)
{
    return OT_VECTEUR_4D(
        valeur[0]-a(0),valeur[1]-a(1),valeur[2]-a(2),valeur[3]-a(3)
    );
}

OT_VECTEUR_4D operator* (const OT_VECTEUR_4D & rkV, const double a)
{
    return OT_VECTEUR_4D(
        a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
    );
}
OT_VECTEUR_4D operator* (const double a, const OT_VECTEUR_4D & rkV)
{                  
    return OT_VECTEUR_4D(
        a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
    );
}
                           
OT_VECTEUR_4D operator/ (const OT_VECTEUR_4D & rkV, const double a)
{
    double inv = 1/a;
    return OT_VECTEUR_4D(
        inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
    );
}


OT_VECTEUR_4D operator/ (const double a, const OT_VECTEUR_4D & rkV)
{
    double inv = 1/a;
    return OT_VECTEUR_4D(
        inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
    );
}

double OT_VECTEUR_4D::operator* (const OT_VECTEUR_4D & a)
{
    return a[0]*valeur[0]+a[1]*valeur[1]+a[2]*valeur[2]+a[3]*valeur[3];
}


OT_VECTEUR_4D& OT_VECTEUR_4D::operator+= (const OT_VECTEUR_4D& rkV)
{
    valeur[0] += rkV.valeur[0];
    valeur[1] += rkV.valeur[1];
    valeur[2] += rkV.valeur[2];
    valeur[3] += rkV.valeur[3];
    return *this;
}

OT_VECTEUR_4D& OT_VECTEUR_4D::operator-= (const OT_VECTEUR_4D& rkV)
{
    valeur[0] -= rkV.valeur[0];
    valeur[1] -= rkV.valeur[1];
    valeur[2] -= rkV.valeur[2];
    valeur[3] -= rkV.valeur[3];
    return *this;
}
OT_VECTEUR_4D& OT_VECTEUR_4D::operator*= (double fScalar)
{
    valeur[0] *= fScalar;
    valeur[1] *= fScalar;
    valeur[2] *= fScalar;
    valeur[3] *= fScalar;
    return *this;
}
OT_VECTEUR_4D& OT_VECTEUR_4D::operator/= (double fScalar)
{               
    double a = 1/fScalar;
    valeur[0] *= a;
    valeur[1] *= a;
    valeur[2] *= a;
    valeur[3] *= a;
    return *this;
}
OT_VECTEUR_4D &OT_VECTEUR_4D::operator= (double a)
{
    valeur[0]=valeur[1]=valeur[2]=valeur[3]=a;
    return *this;
}
                          
std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_4D & __vec)
{
        __os << "{ "<< __vec[0]<<", "<<__vec[1]<<", "<<__vec[2]<<", "<<__vec[3]<<"}, ";
        return __os;
}
#endif // ot_vecteur_4d

int OPERATEUR::egal(double a,double b,double eps)
{
double res=fabs(a-b);
if (res<eps) return(1);
return(0);
}


double  OPERATEUR::qualite_triangle(double* noeud1,double* noeud2,double* noeud3)
{
OT_VECTEUR_3D vec12(noeud1,noeud2);
OT_VECTEUR_3D vec13(noeud1,noeud3);
OT_VECTEUR_3D vec23(noeud2,noeud3);
double d12=vec12.get_longueur();
double d13=vec13.get_longueur();
double d23=vec23.get_longueur();
double dmax=max(d12,d13);
dmax=max(dmax,d23);
double p=0.5*(d12+d13+d23);
double crit=(p-d12)*(p-d13)*(p-d23)/p;
if (crit<0) crit=0.;
crit=sqrt(crit);
crit=crit/dmax;
crit=crit*3.4641101614;
return crit;
}

void OPERATEUR::doubleto2int(double val,int& val1,int& val2)
{
int *p=(int*)(&(val));
val1=(*p);
val2=(*(p+1));
}

double  OPERATEUR::qualite_tetra(double* noeud1,double* noeud2,double* noeud3,double *noeud4)
{
OT_VECTEUR_3D ab(noeud1,noeud2);
OT_VECTEUR_3D ac(noeud1,noeud3);
OT_VECTEUR_3D cb(noeud3,noeud2);
OT_VECTEUR_3D db(noeud4,noeud2);
OT_VECTEUR_3D da(noeud4,noeud1);
OT_VECTEUR_3D dc(noeud4,noeud3);
double dab=ab.get_longueur();
double dac=ac.get_longueur();
double dcb=cb.get_longueur();
double ddb=db.get_longueur();
double dda=da.get_longueur();
double ddc=dc.get_longueur();

double dmax=max(dab,dac);
dmax=max(dmax,dcb);
dmax=max(dmax,ddb);
dmax=max(dmax,dda);
dmax=max(dmax,ddc);

OT_VECTEUR_3D pvec=ab&ac;
double vol=pvec*da;
if (vol>0.) return 0.;
vol=fabs(vol);
double som=pvec.get_longueur()/2.;
pvec=ab&da;
som=som+pvec.get_longueur()/2.;
pvec=ac&dc;
som=som+pvec.get_longueur()/2.;
pvec=cb&db;
som=som+pvec.get_longueur()/2.;
double crit=vol/som/dmax/0.408249;
return crit;
}


template <int N> DOUBLEN<N>::DOUBLEN()
{
tab[0]=-1;
}

template <int N> DOUBLEN<N>::DOUBLEN(DOUBLEN& mdd)
{
for (int i=0;i<N;i++)
        tab[i]=mdd.tab[i];
}

template <int N> DOUBLEN<N>::DOUBLEN(double *v)
{
for (int i=0;i<N;i++)
        tab[i]=v[i];
}

template <int N> DOUBLEN<N>::~DOUBLEN()
{
}

template <int N> double DOUBLEN<N>::get_valeur(int num)
{
return tab[num];
}

template <int N> void DOUBLEN<N>::change_valeur(int num,double val)
{
tab[num]=val;
}

template <int N> DOUBLEN<N> & DOUBLEN<N>::operator=(DOUBLEN<N> & mdd)
{
for (int i=0;i<N;i++)
        tab[i]=mdd.tab[i];
return *this;
}

template <int N> DOUBLEN<N> & DOUBLEN<N>::operator=(const DOUBLEN<N> & mdd)
{
for (int i=0;i<N;i++)
        tab[i]=mdd.tab[i];
return *this;
}


#ifdef BORLANDCPP
       // instanciation pour Windows (Borland c++)
       template class __export DOUBLEN<1>;
       template class __export DOUBLEN<4>;
       template class __export DOUBLEN<10>;
#else
       #ifdef GCC
       // instanciation pour Unix (g++)
       template class DOUBLEN<1>;
       template class DOUBLEN<4>;
       template class DOUBLEN<10>;
       #else
       // instanciation pour Windows (Visual Studio)
       #endif
#endif


Revision:	27
Committed:	Thu Jul 5 15:26:40 2007 UTC (17 years, 10 months ago) by foucault
Original Path:	*magic/lib/outil/outil/src/ot_mathematique.cpp*
File size:	19804 byte(s)
Log Message:
#	User	Rev	Content
1	foucault	27	//------------------------------------------------------------
2			//------------------------------------------------------------
3			// MAGiC
4			// Jean Christophe Cuillière et Vincent FRANCOIS
5			// Département de Génie Mécanique - UQTR
6			//------------------------------------------------------------
7			// Le projet MAGIC est un projet de recherche du département
8			// de génie mécanique de l'Université du Québec à
9			// Trois Rivières
10			// Les librairies ne peuvent être utilisées sans l'accord
11			// des auteurs (contact : francois@uqtr.ca)
12			//------------------------------------------------------------
13			//------------------------------------------------------------
14			//
15			// ot_mathematique.cpp
16			//
17			//------------------------------------------------------------
18			//------------------------------------------------------------
19			// COPYRIGHT 2000
20			// Version du 02/03/2006 à 11H23
21			//------------------------------------------------------------
22			//------------------------------------------------------------
23
24
25			#include "gestionversion.h"
26
27			#include <math.h>
28			#include <algorithm>
29			#include "ot_mathematique.h"
30
31			using namespace std;
32
33
34			OT_VECTEUR_3D::OT_VECTEUR_3D(double x,double y,double z)
35			{
36			valeur[0]=x;
37			valeur[1]=y;
38			valeur[2]=z;
39			}
40
41
42			OT_VECTEUR_3D::OT_VECTEUR_3D(double *xyz)
43			{
44			valeur[0]=xyz[0];
45			valeur[1]=xyz[1];
46			valeur[2]=xyz[2];
47			}
48
49			OT_VECTEUR_3D::OT_VECTEUR_3D(double xyz1,double xyz2)
50			{
51			valeur[0]=xyz2[0]-xyz1[0];
52			valeur[1]=xyz2[1]-xyz1[1];
53			valeur[2]=xyz2[2]-xyz1[2];
54			}
55
56			OT_VECTEUR_3D::OT_VECTEUR_3D(const OT_VECTEUR_3D& mdd)
57			{
58			valeur[0]=mdd.valeur[0];
59			valeur[1]=mdd.valeur[1];
60			valeur[2]=mdd.valeur[2];
61			}
62
63			OT_VECTEUR_3D::OT_VECTEUR_3D(OT_VECTEUR_3D& mdd)
64			{
65			valeur[0]=mdd.valeur[0];
66			valeur[1]=mdd.valeur[1];
67			valeur[2]=mdd.valeur[2];
68			}
69
70
71			double OT_VECTEUR_3D::get_x(void) const
72			{
73			return valeur[0];
74			}
75
76
77			double OT_VECTEUR_3D::get_y(void) const
78			{
79			return valeur[1];
80			}
81
82
83			double OT_VECTEUR_3D::get_z(void) const
84			{
85			return valeur[2];
86			}
87
88			void OT_VECTEUR_3D::change_x(double x)
89			{
90			valeur[0]=x;
91			}
92
93			void OT_VECTEUR_3D::change_y(double y)
94			{
95			valeur[1]=y;
96			}
97
98			void OT_VECTEUR_3D::change_z(double z)
99			{
100			valeur[2]=z;
101			}
102
103			double* OT_VECTEUR_3D::get_xyz(void)
104			{
105			return valeur;
106			}
107
108			OT_VECTEUR_3D::operator const double* () const
109			{
110			return valeur;
111			}
112
113			OT_VECTEUR_3D::operator double* ()
114			{
115			return valeur;
116			}
117
118			double OT_VECTEUR_3D::operator[] (int i) const
119			{
120			return valeur[i];
121			}
122
123			double & OT_VECTEUR_3D::operator[] (int i)
124			{
125			return valeur[i];
126			}
127
128			double OT_VECTEUR_3D::operator() (int i) const
129			{
130			return valeur[i];
131			}
132			double & OT_VECTEUR_3D::operator() (int i)
133			{
134			return valeur[i];
135			}
136
137			double OT_VECTEUR_3D::get_longueur(void) const
138			{
139			return sqrt(valeur[0]valeur[0]+valeur[1]valeur[1]+valeur[2]*valeur[2]);
140			}
141
142			double OT_VECTEUR_3D::get_longueur2(void) const
143			{
144			return (valeur[0]valeur[0]+valeur[1]valeur[1]+valeur[2]*valeur[2]);
145			}
146
147
148			void OT_VECTEUR_3D::norme(void)
149			{
150			double longueur=get_longueur();
151			if (longueur != 0)
152			{
153			valeur[0]=valeur[0]/longueur;
154			valeur[2]=valeur[2]/longueur;
155			valeur[1]=valeur[1]/longueur;
156			}
157			}
158
159			int OT_VECTEUR_3D::compare_valeur (const OT_VECTEUR_3D& mdd) const
160			{
161			return memcmp(valeur, mdd.valeur, 3*sizeof(double));
162			}
163
164			bool OT_VECTEUR_3D::operator== (const OT_VECTEUR_3D& mdd) const
165			{
166			return compare_valeur(mdd) == 0;
167			}
168
169			bool OT_VECTEUR_3D::operator!= (const OT_VECTEUR_3D& mdd) const
170			{
171			return compare_valeur(mdd) != 0;
172			}
173
174			bool OT_VECTEUR_3D::operator< (const OT_VECTEUR_3D& mdd) const
175			{
176			return compare_valeur(mdd) < 0;
177			}
178
179			bool OT_VECTEUR_3D::operator<= (const OT_VECTEUR_3D& mdd) const
180			{
181			return compare_valeur(mdd) <= 0;
182			}
183
184			bool OT_VECTEUR_3D::operator> (const OT_VECTEUR_3D& mdd) const
185			{
186			return compare_valeur(mdd) > 0;
187			}
188
189			bool OT_VECTEUR_3D::operator>= (const OT_VECTEUR_3D& mdd) const
190			{
191			return compare_valeur(mdd) >= 0;
192			}
193
194			std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_3D & __vec)
195			{
196			__os << "{ "<< __vec[0]<<", "<<__vec[1]<<", "<<__vec[2]<<"}, ";
197			return __os;
198			}
199
200			OT_VECTEUR_3D operator/ (const OT_VECTEUR_3D& mdd1, double diviseur)
201			{
202			OT_VECTEUR_3D kQuot;
203
204			if ( diviseur != (double)0.0 )
205			{
206			double facteur = ((double)1.0)/diviseur;
207			kQuot.valeur[0] = facteur*mdd1.valeur[0];
208			kQuot.valeur[1] = facteur*mdd1.valeur[1];
209			kQuot.valeur[2] = facteur*mdd1.valeur[2];
210			}
211			else
212			{
213			kQuot.valeur[0] = 1E300;
214			kQuot.valeur[1] = 1E300;
215			kQuot.valeur[2] = 1E300;
216			}
217
218			return kQuot;
219			}
220
221			OT_VECTEUR_3D operator- (const OT_VECTEUR_3D& mdd1)
222			{
223			return OT_VECTEUR_3D(
224			-mdd1.valeur[0],
225			-mdd1.valeur[1],
226			-mdd1.valeur[2]);
227			}
228
229			double operator*(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
230			{
231			double tmp;
232			tmp=mdd1.valeur[0]mdd2.valeur[0]+mdd1.valeur[1]mdd2.valeur[1]+mdd1.valeur[2]*mdd2.valeur[2];
233			return tmp;
234			}
235
236			OT_VECTEUR_3D operator*(double mdd1,const OT_VECTEUR_3D& mdd2)
237			{
238			OT_VECTEUR_3D tmp;
239			tmp.valeur[0]=mdd1*mdd2.valeur[0];
240			tmp.valeur[1]=mdd1*mdd2.valeur[1];
241			tmp.valeur[2]=mdd1*mdd2.valeur[2];
242			return tmp;
243			}
244
245			OT_VECTEUR_3D& OT_VECTEUR_3D::operator+= (const OT_VECTEUR_3D& mdd)
246			{
247			valeur[0] += mdd.valeur[0];
248			valeur[1] += mdd.valeur[1];
249			valeur[2] += mdd.valeur[2];
250			return *this;
251			}
252			OT_VECTEUR_3D& OT_VECTEUR_3D::operator-= (const OT_VECTEUR_3D& mdd)
253			{
254			valeur[0] -= mdd.valeur[0];
255			valeur[1] -= mdd.valeur[1];
256			valeur[2] -= mdd.valeur[2];
257			return *this;
258			}
259			OT_VECTEUR_3D& OT_VECTEUR_3D::operator*= (double facteur)
260			{
261			valeur[0] *= facteur;
262			valeur[1] *= facteur;
263			valeur[2] *= facteur;
264			return *this;
265			}
266			OT_VECTEUR_3D& OT_VECTEUR_3D::operator/= (double diviseur)
267			{
268			if ( diviseur != (double)0.0 )
269			{
270			double facteur = ((double)1.0)/diviseur;
271			valeur[0] *= facteur;
272			valeur[1] *= facteur;
273			valeur[2] *= facteur;
274			}
275			else
276			{
277			valeur[0] = 1E300;
278			valeur[1] = 1E300;
279			valeur[2] = 1E300;
280			}
281
282			return *this;
283			}
284
285			OT_VECTEUR_3D operator&(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
286			{
287			OT_VECTEUR_3D tmp;
288			tmp.valeur[0]=mdd1.valeur[1]mdd2.valeur[2]-mdd1.valeur[2]mdd2.valeur[1];
289			tmp.valeur[1]=mdd1.valeur[2]mdd2.valeur[0]-mdd1.valeur[0]mdd2.valeur[2];
290			tmp.valeur[2]=mdd1.valeur[0]mdd2.valeur[1]-mdd1.valeur[1]mdd2.valeur[0];
291			return tmp;
292			}
293
294			double OT_VECTEUR_3D::diff(void)
295			{
296			return 0.33333333333333*(fabs(valeur[0])+fabs(valeur[1])+fabs(valeur[2]));
297			}
298
299
300			OT_MATRICE_3D::OT_MATRICE_3D(OT_VECTEUR_3D& colonne1,OT_VECTEUR_3D& colonne2,OT_VECTEUR_3D& colonne3)
301			{
302			vec[0]=colonne1;
303			vec[1]=colonne2;
304			vec[2]=colonne3;
305			}
306
307			OT_MATRICE_3D::OT_MATRICE_3D(const OT_MATRICE_3D& mdd)
308			{
309			vec[0]=mdd.vec[0];
310			vec[1]=mdd.vec[1];
311			vec[2]=mdd.vec[2];
312			}
313
314
315			std::ostream & operator << ( std::ostream & __os, const OT_MATRICE_3D & __mat)
316			{
317
318			__os << "{ ";
319			for (unsigned i=0; i<3; i++)
320			{
321			__os << "[ "<< __mat(i,0)<<" ]\t [ "<<__mat(i,1) <<" ]\t [ "<<__mat(i,2);
322			if (i+1 < 3)
323			__os << std::endl << " ";
324			}
325			__os << " }" << std::endl;
326			return __os;
327			}
328
329			OT_VECTEUR_3D & OT_MATRICE_3D::operator [](int i)
330			{
331			return vec[i];
332			}
333
334			double OT_MATRICE_3D::valeur(int iLigne, int iCol) const
335			{
336			return vec[iCol][iLigne];
337			}
338			double & OT_MATRICE_3D::valeur(int iLigne, int iCol)
339			{
340			return vec[iCol][iLigne];
341			}
342			double OT_MATRICE_3D::operator() (int iLigne, int iCol) const
343			{
344			return vec[iCol][iLigne];
345			}
346			double & OT_MATRICE_3D::operator() (int iLigne, int iCol)
347			{
348			return vec[iCol][iLigne];
349			}
350			double OT_MATRICE_3D::get_determinant()
351			{
352			double coFacteur0 = valeur(1,1)valeur(2,2) - valeur(1,2)valeur(2,1);
353			double coFacteur1 = valeur(1,2)valeur(2,0) - valeur(1,0)valeur(2,2);
354			double coFacteur2 = valeur(1,0)valeur(2,1) - valeur(1,1)valeur(2,0);
355			double determinant = valeur(0,0)coFacteur0 + valeur(0,1)coFacteur1 + valeur(0,2)*coFacteur2;
356			return determinant;
357			}
358			void OT_MATRICE_3D::identite(void)
359			{
360			for (unsigned i = 0; i<3; i++)
361			for (unsigned j = 0; j < 3; j++)
362			vec[i][j] = (i==j) ? 1.0 : 0.0 ;
363			}
364
365			void OT_MATRICE_3D::transpose(OT_MATRICE_3D& res) const
366			{
367			for (unsigned i = 0; i<3; i++)
368			for (unsigned j = 0; j < 3; j++)
369			res.vec[i][j] = vec[j][i];
370			}
371
372			OT_MATRICE_3D OT_MATRICE_3D::inverse () const
373			{
374			// Invert a 3x3 using cofactors. This is faster than using a generic
375			// Gaussian elimination because of the loop overhead of such a method.
376
377			OT_MATRICE_3D matInverse;
378
379			matInverse.valeur(0,0) = valeur(1,1)valeur(2,2) - valeur(1,2)valeur(2,1);
380			matInverse.valeur(0,1) = valeur(0,2)valeur(2,1) - valeur(0,1)valeur(2,2);
381			matInverse.valeur(0,2) = valeur(0,1)valeur(1,2) - valeur(0,2)valeur(1,1);
382			matInverse.valeur(1,0) = valeur(1,2)valeur(2,0) - valeur(1,0)valeur(2,2);
383			matInverse.valeur(1,1) = valeur(0,0)valeur(2,2) - valeur(0,2)valeur(2,0);
384			matInverse.valeur(1,2) = valeur(0,2)valeur(1,0) - valeur(0,0)valeur(1,2);
385			matInverse.valeur(2,0) = valeur(1,0)valeur(2,1) - valeur(1,1)valeur(2,0);
386			matInverse.valeur(2,1) = valeur(0,1)valeur(2,0) - valeur(0,0)valeur(2,1);
387			matInverse.valeur(2,2) = valeur(0,0)valeur(1,1) - valeur(0,1)valeur(1,0);
388
389			double det = valeur(0,0)matInverse(0,0) + valeur(0,1)matInverse(1,0)+
390			valeur(0,2)*matInverse(2,0);
391
392			if ( fabs(det) <= 1E-100 )
393			{
394			matInverse.identite();
395			return matInverse;
396			}
397
398			double invDet = 1/det;
399
400			for (int i=0; i<3; i++)
401			for (int j=0; j<3; j++)
402			matInverse.vec[i][j] *= invDet;
403
404			return matInverse;
405			}
406
407			void OT_MATRICE_3D::change_vecteur1(OT_VECTEUR_3D v)
408			{
409			vec[0]=v;
410			}
411			void OT_MATRICE_3D::change_vecteur2(OT_VECTEUR_3D v)
412			{
413			vec[1]=v;
414			}
415			void OT_MATRICE_3D::change_vecteur3(OT_VECTEUR_3D v)
416			{
417			vec[2]=v;
418			}
419			OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur1(void)
420			{
421			return vec[0];
422			}
423			OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur2(void)
424			{
425			return vec[1];
426			}
427			OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur3(void)
428			{
429			return vec[2];
430			}
431
432
433
434			OT_MATRICE_3D OT_MATRICE_3D::VecteursPropres(OT_VECTEUR_3D& Val_Propre,double Eps,int ItMax)
435			{
436			OT_MATRICE_3D m=*this;
437			OT_VECTEUR_3D ZERO(0.,0.,0.) ;
438			OT_MATRICE_3D vct_prpres(ZERO,ZERO,ZERO);
439			OT_VECTEUR_3D V ;
440
441			int Dim=3;
442			int I,J,K;
443			OT_VECTEUR_3D U, Up;
444			double N,T;
445			//-----
446			for(int vi=0;vi<3;vi++)
447			{
448			//===================
449			Up=vec[0]; // il faut l'initialiser d'une autre maniere
450			K=0;
451			do{
452			U=Up;
453			for(int I=0 ;I<Dim;I++){
454			T=0;
455			for(int J=0; J<Dim;J++) T=T+valeur(I,J)*U[J];
456			Up[I]=T ;
457			}
458			N=Up.get_longueur();
459			if( Up[0]<0) N=-N;
460			T=0;
461			for(int I=0;I<Dim;I++){
462			Up[I]=Up[I]/N;
463			T=T+fabs(Up[I]-U[I]);
464			}
465			K=K+1;
466			}
467			while((K<ItMax)\|\|(T>Eps)) ;
468			V=Up ;
469			vct_prpres[vi]=V;
470			//================
471			double N,D;
472			N=0.0; D=0.0;
473			for(int I=0;I<Dim;I++){
474
475			D=D+V[I]*V[I];
476			for(int J=0 ;J<Dim;J++) N=N+V[I]valeur(I,J)V[J];
477			}
478			double Val=N/D;
479			Val_Propre[vi]=Val;
480			//===============
481			for(int I=0; I<Dim ;I++){
482
483			for(int J=0;J<Dim;J++) m(I,J)=m(I,J)-ValV[I]V[J];
484			}
485
486			//===============
487
488			}
489			//-----------------
490
491
492			return vct_prpres;
493			}
494
495			/*
496			double OT_MATRICE_3D::ValeurPropre(double Eps,int ItMax)
497			{
498			double Val_Propre;
499			OT_VECTEUR_3D V=VecteurPropre( Eps, ItMax);
500			int I,J,Dim;
501			double N,D;
502			N=0.0; D=0.0; Dim=3;
503			for(int I=0;I<Dim;I++){
504
505			D=D+V[I]*V[I];
506			for(int J=0 ;J<Dim;J++) N=N+V[I]valeur(I,J)V[J];
507			}
508			Val_Propre=N/D;
509			return Val_Propre;
510			}
511
512			void OT_MATRICE_3D:: Deflation(OT_VECTEUR_3D& V,double val_prpre)
513			{
514			int Dim=3;
515			for(int I=0; I<Dim ;I++){
516
517			for(int J=0;J<Dim;J++) valeur(I,J)=valeur(I,J)-val_prpreV[I]V[J];
518			}
519
520			} */
521
522
523
524			OT_VECTEUR_3D operator*(const OT_MATRICE_3D& mdd1,const OT_VECTEUR_3D& mdd2)
525			{
526			OT_VECTEUR_3D tmp;
527			tmp(0)=mdd1(0,0)mdd2(0)+mdd1(0,1)mdd2(1)+mdd1(0,2)*mdd2(2);
528			tmp(1)=mdd1(1,0)mdd2(0)+mdd1(1,1)mdd2(1)+mdd1(1,2)*mdd2(2);
529			tmp(2)=mdd1(2,0)mdd2(0)+mdd1(2,1)mdd2(1)+mdd1(2,2)*mdd2(2);
530			return tmp;}
531
532			OT_VECTEUR_3D operator+(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
533			{
534			OT_VECTEUR_3D tmp;
535			tmp(0)=mdd1(0)+mdd2(0);
536			tmp(1)=mdd1(1)+mdd2(1);
537			tmp(2)=mdd1(2)+mdd2(2);
538			return tmp;
539			}
540
541			OT_VECTEUR_3D operator-(const class OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
542			{
543			OT_VECTEUR_3D tmp;
544			tmp.valeur[0]=mdd1.valeur[0]-mdd2.valeur[0];
545			tmp.valeur[1]=mdd1.valeur[1]-mdd2.valeur[1];
546			tmp.valeur[2]=mdd1.valeur[2]-mdd2.valeur[2];
547			return tmp;
548			}
549
550			OT_MATRICE_3D operator*(const OT_MATRICE_3D& mdd1,const OT_MATRICE_3D& mdd2)
551			{
552			OT_MATRICE_3D tmp;
553			for (int i=0; i<3; i++)
554			for (int j=0; j<3; j++)
555			{
556			tmp(i,j) = 0;
557			for (int k=0; k<3; k++)
558			tmp(i,j) += mdd1(i,k)*mdd2(k,j);
559			}
560			return tmp;
561			}
562
563			#ifdef ot_vecteur_4d
564			OT_VECTEUR_4D::OT_VECTEUR_4D (double v[4])
565			{
566			for (int i=0; i<4; i++)
567			valeur[i]=v[i];
568			}
569
570			OT_VECTEUR_4D::OT_VECTEUR_4D (double x, double y, double z, double w)
571			{
572			valeur[0] = x;
573			valeur[1] = y;
574			valeur[2] = z;
575			valeur[3] = w;
576			}
577
578			OT_VECTEUR_4D::OT_VECTEUR_4D (const OT_VECTEUR_4D & v)
579			{
580			for (int i=0; i<4; i++)
581			valeur[i]=v[i];
582			}
583
584
585			double & OT_VECTEUR_4D::operator()(int i)
586			{
587			return valeur[i];
588			}
589
590			double OT_VECTEUR_4D::operator()(int i) const
591			{
592			return valeur[i];
593			}
594
595			double & OT_VECTEUR_4D::operator[](int i)
596			{
597			return valeur[i];
598			}
599
600			double OT_VECTEUR_4D::operator[](int i) const
601			{
602			return valeur[i];
603			}
604
605			OT_VECTEUR_4D::operator const double* () const
606			{
607			return valeur;
608			}
609
610			OT_VECTEUR_4D::operator double* ()
611			{
612			return valeur;
613			}
614
615			double OT_VECTEUR_4D::get_x(void) const
616			{ return valeur[0]; }
617			double & OT_VECTEUR_4D::x(void)
618			{ return valeur[0]; }
619			double OT_VECTEUR_4D::get_y(void) const
620			{ return valeur[1]; }
621			double & OT_VECTEUR_4D::y(void)
622			{ return valeur[1]; }
623			double OT_VECTEUR_4D::get_z(void) const
624			{ return valeur[2]; }
625			double & OT_VECTEUR_4D::z(void)
626			{ return valeur[2]; }
627			double OT_VECTEUR_4D::get_w(void) const
628			{ return valeur[3]; }
629			double & OT_VECTEUR_4D::w(void)
630			{ return valeur[3]; }
631			void OT_VECTEUR_4D::change_x(double x)
632			{ valeur[0] = x; }
633			void OT_VECTEUR_4D::change_y(double y)
634			{ valeur[1] = y; }
635			void OT_VECTEUR_4D::change_z(double z)
636			{ valeur[2] = z; }
637			void OT_VECTEUR_4D::change_w(double w)
638			{ valeur[3] = w; }
639			double* OT_VECTEUR_4D::get_xyzw(void)
640			{ return valeur; }
641
642			OT_VECTEUR_4D OT_VECTEUR_4D::operator+ (const OT_VECTEUR_4D & a)
643			{
644			return OT_VECTEUR_4D(
645			a(0)+valeur[0],a(1)+valeur[1],a(2)+valeur[2],a(3)+valeur[3]
646			);
647			}
648
649			OT_VECTEUR_4D OT_VECTEUR_4D::operator- (const OT_VECTEUR_4D & a)
650			{
651			return OT_VECTEUR_4D(
652			valeur[0]-a(0),valeur[1]-a(1),valeur[2]-a(2),valeur[3]-a(3)
653			);
654			}
655
656			OT_VECTEUR_4D operator* (const OT_VECTEUR_4D & rkV, const double a)
657			{
658			return OT_VECTEUR_4D(
659			arkV[0],arkV[1],arkV[2],arkV[3]
660			);
661			}
662			OT_VECTEUR_4D operator* (const double a, const OT_VECTEUR_4D & rkV)
663			{
664			return OT_VECTEUR_4D(
665			arkV[0],arkV[1],arkV[2],arkV[3]
666			);
667			}
668
669			OT_VECTEUR_4D operator/ (const OT_VECTEUR_4D & rkV, const double a)
670			{
671			double inv = 1/a;
672			return OT_VECTEUR_4D(
673			invrkV[0],invrkV[1],invrkV[2],invrkV[3]
674			);
675			}
676
677
678			OT_VECTEUR_4D operator/ (const double a, const OT_VECTEUR_4D & rkV)
679			{
680			double inv = 1/a;
681			return OT_VECTEUR_4D(
682			invrkV[0],invrkV[1],invrkV[2],invrkV[3]
683			);
684			}
685
686			double OT_VECTEUR_4D::operator* (const OT_VECTEUR_4D & a)
687			{
688			return a[0]valeur[0]+a[1]valeur[1]+a[2]valeur[2]+a[3]valeur[3];
689			}
690
691
692			OT_VECTEUR_4D& OT_VECTEUR_4D::operator+= (const OT_VECTEUR_4D& rkV)
693			{
694			valeur[0] += rkV.valeur[0];
695			valeur[1] += rkV.valeur[1];
696			valeur[2] += rkV.valeur[2];
697			valeur[3] += rkV.valeur[3];
698			return *this;
699			}
700
701			OT_VECTEUR_4D& OT_VECTEUR_4D::operator-= (const OT_VECTEUR_4D& rkV)
702			{
703			valeur[0] -= rkV.valeur[0];
704			valeur[1] -= rkV.valeur[1];
705			valeur[2] -= rkV.valeur[2];
706			valeur[3] -= rkV.valeur[3];
707			return *this;
708			}
709			OT_VECTEUR_4D& OT_VECTEUR_4D::operator*= (double fScalar)
710			{
711			valeur[0] *= fScalar;
712			valeur[1] *= fScalar;
713			valeur[2] *= fScalar;
714			valeur[3] *= fScalar;
715			return *this;
716			}
717			OT_VECTEUR_4D& OT_VECTEUR_4D::operator/= (double fScalar)
718			{
719			double a = 1/fScalar;
720			valeur[0] *= a;
721			valeur[1] *= a;
722			valeur[2] *= a;
723			valeur[3] *= a;
724			return *this;
725			}
726			OT_VECTEUR_4D &OT_VECTEUR_4D::operator= (double a)
727			{
728			valeur[0]=valeur[1]=valeur[2]=valeur[3]=a;
729			return *this;
730			}
731
732			std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_4D & __vec)
733			{
734			__os << "{ "<< __vec[0]<<", "<<__vec[1]<<", "<<__vec[2]<<", "<<__vec[3]<<"}, ";
735			return __os;
736			}
737			#endif // ot_vecteur_4d
738
739			int OPERATEUR::egal(double a,double b,double eps)
740			{
741			double res=fabs(a-b);
742			if (res<eps) return(1);
743			return(0);
744			}
745
746
747			double OPERATEUR::qualite_triangle(double* noeud1,double* noeud2,double* noeud3)
748			{
749			OT_VECTEUR_3D vec12(noeud1,noeud2);
750			OT_VECTEUR_3D vec13(noeud1,noeud3);
751			OT_VECTEUR_3D vec23(noeud2,noeud3);
752			double d12=vec12.get_longueur();
753			double d13=vec13.get_longueur();
754			double d23=vec23.get_longueur();
755			double dmax=max(d12,d13);
756			dmax=max(dmax,d23);
757			double p=0.5*(d12+d13+d23);
758			double crit=(p-d12)(p-d13)(p-d23)/p;
759			if (crit<0) crit=0.;
760			crit=sqrt(crit);
761			crit=crit/dmax;
762			crit=crit*3.4641101614;
763			return crit;
764			}
765
766			void OPERATEUR::doubleto2int(double val,int& val1,int& val2)
767			{
768			int p=(int)(&(val));
769			val1=(*p);
770			val2=(*(p+1));
771			}
772
773			double OPERATEUR::qualite_tetra(double* noeud1,double* noeud2,double* noeud3,double *noeud4)
774			{
775			OT_VECTEUR_3D ab(noeud1,noeud2);
776			OT_VECTEUR_3D ac(noeud1,noeud3);
777			OT_VECTEUR_3D cb(noeud3,noeud2);
778			OT_VECTEUR_3D db(noeud4,noeud2);
779			OT_VECTEUR_3D da(noeud4,noeud1);
780			OT_VECTEUR_3D dc(noeud4,noeud3);
781			double dab=ab.get_longueur();
782			double dac=ac.get_longueur();
783			double dcb=cb.get_longueur();
784			double ddb=db.get_longueur();
785			double dda=da.get_longueur();
786			double ddc=dc.get_longueur();
787
788			double dmax=max(dab,dac);
789			dmax=max(dmax,dcb);
790			dmax=max(dmax,ddb);
791			dmax=max(dmax,dda);
792			dmax=max(dmax,ddc);
793
794			OT_VECTEUR_3D pvec=ab&ac;
795			double vol=pvec*da;
796			if (vol>0.) return 0.;
797			vol=fabs(vol);
798			double som=pvec.get_longueur()/2.;
799			pvec=ab&da;
800			som=som+pvec.get_longueur()/2.;
801			pvec=ac&dc;
802			som=som+pvec.get_longueur()/2.;
803			pvec=cb&db;
804			som=som+pvec.get_longueur()/2.;
805			double crit=vol/som/dmax/0.408249;
806			return crit;
807			}
808
809
810
811			template <int N> DOUBLEN<N>::DOUBLEN()
812			{
813			tab[0]=-1;
814			}
815
816			template <int N> DOUBLEN<N>::DOUBLEN(DOUBLEN& mdd)
817			{
818			for (int i=0;i<N;i++)
819			tab[i]=mdd.tab[i];
820			}
821
822			template <int N> DOUBLEN<N>::DOUBLEN(double *v)
823			{
824			for (int i=0;i<N;i++)
825			tab[i]=v[i];
826			}
827
828			template <int N> DOUBLEN<N>::~DOUBLEN()
829			{
830			}
831
832			template <int N> double DOUBLEN<N>::get_valeur(int num)
833			{
834			return tab[num];
835			}
836
837			template <int N> void DOUBLEN<N>::change_valeur(int num,double val)
838			{
839			tab[num]=val;
840			}
841
842			template <int N> DOUBLEN<N> & DOUBLEN<N>::operator=(DOUBLEN<N> & mdd)
843			{
844			for (int i=0;i<N;i++)
845			tab[i]=mdd.tab[i];
846			return *this;
847			}
848
849			template <int N> DOUBLEN<N> & DOUBLEN<N>::operator=(const DOUBLEN<N> & mdd)
850			{
851			for (int i=0;i<N;i++)
852			tab[i]=mdd.tab[i];
853			return *this;
854			}
855
856
857			#ifdef BORLANDCPP
858			// instanciation pour Windows (Borland c++)
859			template class __export DOUBLEN<1>;
860			template class __export DOUBLEN<4>;
861			template class __export DOUBLEN<10>;
862			#else
863			#ifdef GCC
864			// instanciation pour Unix (g++)
865			template class DOUBLEN<1>;
866			template class DOUBLEN<4>;
867			template class DOUBLEN<10>;
868			#else
869			// instanciation pour Windows (Visual Studio)
870			#endif
871			#endif
872
873