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(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_3d


OT_VECTEUR_3DD::OT_VECTEUR_3DD (double2 v[3])
{
    for (int i=0; i<3; i++)
        valeur[i]=v[i];
}

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

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


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

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

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

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

OT_VECTEUR_3DD::operator const double2* () const
{
    return valeur;
}

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

OT_VECTEUR_3DD& OT_VECTEUR_3DD::operator =(const OT_VECTEUR_3DD& v3d)
{
    if (this==&v3d) return *this;
    valeur [0]=v3d.get_x();
    valeur [1]=v3d.get_y();
    valeur [2]=v3d.get_z();
    return *this;
}
double2 OT_VECTEUR_3DD::norme()
{
    double2 nor= valeur [0]*valeur [0]+valeur [1]*valeur [1]+valeur [2]*valeur [2];
    nor=sqrt(nor);
    return nor;
}
OT_VECTEUR_3DD OT_VECTEUR_3DD::vecteur_norme()
{
    OT_VECTEUR_3DD nor;
    double2 zero=0.;
    double2 norme= this->norme() ;
    nor=(*this)*(1./norme);
    return nor;
}

double2 OT_VECTEUR_3DD::get_x(void) const
{
    return valeur[0];
}
double2 & OT_VECTEUR_3DD::x(void)
{
    return valeur[0];
}
double2 OT_VECTEUR_3DD::get_y(void) const
{
    return valeur[1];
}
double2 & OT_VECTEUR_3DD::y(void)
{
    return valeur[1];
}
double2 OT_VECTEUR_3DD::get_z(void) const
{
    return valeur[2];
}
double2 & OT_VECTEUR_3DD::z(void)
{
    return valeur[2];
}
void OT_VECTEUR_3DD::change_x(double2 x)
{
    valeur[0] = x;
}
void OT_VECTEUR_3DD::change_y(double2 y)
{
    valeur[1] = y;
}
void OT_VECTEUR_3DD::change_z(double2 z)
{
    valeur[2] = z;
}
double2* OT_VECTEUR_3DD::get_xyz(void)
{
    return valeur;
}

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

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

OT_VECTEUR_3DD operator* (const OT_VECTEUR_3DD & rkV, const double2 a)
{
    return OT_VECTEUR_3DD(
               a*rkV[0],a*rkV[1],a*rkV[2]);
}
OT_VECTEUR_3DD operator* (const double2 a, const OT_VECTEUR_3DD & rkV)
{
    return OT_VECTEUR_3DD(
               a*rkV[0],a*rkV[1],a*rkV[2]);
}

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


OT_VECTEUR_3DD operator/ (const double2 a, const OT_VECTEUR_3DD & rkV)
{
    double2 inv = 1/a;
    return OT_VECTEUR_3DD(
               inv*rkV[0],inv*rkV[1],inv*rkV[2]);
}
OT_VECTEUR_3DD OT_VECTEUR_3DD::operator&(const OT_VECTEUR_3DD& v2)
{
    OT_VECTEUR_3DD 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];
    return tmp;
}
double2 OT_VECTEUR_3DD::operator* (const OT_VECTEUR_3DD & a)
{
    return a[0]*valeur[0]+a[1]*valeur[1]+a[2]*valeur[2];
}


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

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


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

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

}

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


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

OT_VECTEUR_4DD OT_VECTEUR_4DD::vecteur_norme_3d()
{
    double2 norme= valeur [0]*valeur [0]+valeur [1]*valeur [1]+valeur [2]*valeur [2];
    norme=norme^0.5;
    OT_VECTEUR_4DD nor;
    double2 zero=0.;
    assert(norme!=zero) ;
    nor.change_x(valeur[0]/norme);
    nor.change_y(valeur[1]/norme);
    nor.change_z(valeur[2]/norme);
    nor.change_w(valeur[3]);
    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:	363
Committed:	Thu Oct 18 22:12:03 2012 UTC (12 years, 6 months ago) by francois
File size:	29685 byte(s)
Log Message:	Mise a jour pour publi sur la comparaison