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 <iomanip>
#include <math.h>
#include <string.h>
#include <assert.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];
}

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

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;
 }
}

OT_VECTEUR_3D OT_VECTEUR_3D::gram_shmidt(const OT_VECTEUR_3D& vint)
{
 OT_VECTEUR_3D v(0.,0.,0.);
 double alpha=(*this)*vint;
 v= *this-alpha*vint;
 return v;
}

OT_VECTEUR_3D OT_VECTEUR_3D::unite(int i)
{
 OT_VECTEUR_3D v(0.,0.,0.);
 v[i]=1.;
 return v; 
}

                                                              
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];
}

OT_MATRICE_3D::OT_MATRICE_3D(void)
 {
  for(int i=0;i<3;i++)
   {
     for(int j=0;j<3;j++)
      {
      vec[i][j]=0.;
      }

   }

 }
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];
}

OT_MATRICE_3D::OT_MATRICE_3D(double* t)
{
  for(int i=0;i<3;i++)
   {
     for(int j=0;j<3;j++)
      {
      vec[i][j]=t[i*3+j];
      }

   }
}

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::transpose() const
{
        OT_MATRICE_3D temp;
        for (unsigned i = 0; i<3; i++)
                for (unsigned j = 0; j < 3; j++)
                        temp.valeur(i,j) =valeur(j,i);

return temp;                        
}

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 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)= mdd1(i,j)+mdd2(i,j);
   }
   return tmp;
}


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


OT_VECTEUR_4DD::OT_VECTEUR_4DD (double2 v[4])
{
    for (int i=0; i<4; i++)
        valeur[i]=v[i];
}
                      
OT_VECTEUR_4DD::OT_VECTEUR_4DD (double2 x, double2 y, double2 z, double2 w)
{
    valeur[0] = x;
    valeur[1] = y;
    valeur[2] = z;
    valeur[3] = w;
}

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


double2 & OT_VECTEUR_4DD::operator()(int i)
{
    return valeur[i];
}

double2 OT_VECTEUR_4DD::operator()(int i) const
{
    return valeur[i];    
}

double2 & OT_VECTEUR_4DD::operator[](int i)
{
    return valeur[i];
}

double2 OT_VECTEUR_4DD::operator[](int i) const
{
    return valeur[i];    
}
               
OT_VECTEUR_4DD::operator const double2* () const
{
    return valeur;
}

OT_VECTEUR_4DD::operator double2* ()
{
    return valeur;
}

OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator =(const OT_VECTEUR_4DD& v4d)
{
    if (this==&v4d) return *this;
    valeur [0]=v4d.get_x();
    valeur [1]=v4d.get_y();
    valeur [2]=v4d.get_z();
    valeur [3]=v4d.get_w();
  return *this;
}
double2 OT_VECTEUR_4DD::norme()
{
double2 nor= valeur [0]*valeur [0]+valeur [1]*valeur [1]+valeur [2]*valeur [2]+valeur [3]*valeur [3];
nor=nor^0.5;
return nor;
}
OT_VECTEUR_4DD OT_VECTEUR_4DD::vecteur_norme()
{
OT_VECTEUR_4DD nor;
double2 zero=0.;
double2 norme= this->norme() ;
assert(norme!=zero) ;
nor=(*this)*(1./norme);
return nor;
}

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

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

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


OT_VECTEUR_4DD operator/ (const double2 a, const OT_VECTEUR_4DD & rkV)
{
    double2 inv = 1/a;
    return OT_VECTEUR_4DD(
        inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
    );
}
OT_VECTEUR_4DD OT_VECTEUR_4DD::operator^(const OT_VECTEUR_4DD& v2)
{
OT_VECTEUR_4DD tmp;
  tmp[0]= valeur[1]*v2[2]-valeur[2]*v2[1];
  tmp[1]= valeur[2]*v2[0]-valeur[0]*v2[2];
  tmp[2]= valeur[0]*v2[1]-valeur[1]*v2[0];
  tmp[3]= 0;
  //tmp[3]= valeur[1]*v2[2]-valeur[2]*v2[1];
  return tmp;
}
double2 OT_VECTEUR_4DD::operator* (const OT_VECTEUR_4DD & a)
{
    return a[0]*valeur[0]+a[1]*valeur[1]+a[2]*valeur[2]+a[3]*valeur[3];
}


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

OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator-= (const OT_VECTEUR_4DD& rkV)
{
    valeur[0] = valeur[0] -rkV.valeur[0];
    valeur[1] = valeur[1] -rkV.valeur[1];
    valeur[2] = valeur[2] -rkV.valeur[2];
    valeur[3] = valeur[3] -rkV.valeur[3];
    return *this;
}
OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator*= (double2 fScalar)
{
    valeur[0] =valeur[0] * fScalar;
    valeur[1] =valeur[1] * fScalar;
    valeur[2] =valeur[2] * fScalar;
    valeur[3] =valeur[3] * fScalar;
    return *this;
}
OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator/= (double2 fScalar)
{               
    double2 a = 1/fScalar;
    valeur[0] =valeur[0] * a;
    valeur[1] =valeur[1] * a;
    valeur[2] =valeur[2] * a;
    valeur[3] =valeur[3] * a;
    return *this;
}


std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_4DD & __vec)
{
        __os << "[ "<<setw(20)<< __vec[0]<<setw(4)<<", "<<setw(20)<<__vec[1]<<", "<<setw(20)<<__vec[2]<<setw(4)<<", "<<setw(20)<<__vec[3]<<setw(4)<<"] ";
        return __os;
}

int  OT_VECTEUR_4DD::compare_valeur (const OT_VECTEUR_4DD& mdd) const
{
 if (valeur[0]==mdd.valeur[0]&&valeur[1]==mdd.valeur[1] &&valeur[2]==mdd.valeur[2]&&valeur[3]==mdd.valeur[3])
 return 1;
 else return 0;

}

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


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:	169
Committed:	Fri Feb 13 23:04:14 2009 UTC (16 years, 3 months ago) by francois
Original Path:	*magic/lib/outil/outil/src/ot_mathematique.cpp*
File size:	25660 byte(s)
Log Message:	Resolution de bug avec la version de gcc 4.3.2 --> mise la norme de c++