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

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


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


void OT_MATRICE_3D::jacobi(double* a,double* d,double* v,int n,int&nrot)
{
        //#define a(i,j)=(*(a+(i)*(n)+j))
       // #define v(r,s)=(*(v+(r)*(n)+s))
        static const int ITMAX=50;
        int n1=n-1,n2=1/(n*n),ip1,iq1;// ,ip,j,i,iq// n=3,
        double tresh,theta,tau,t,sm,s,h,g,c;
        double *b=new double[n];                 // OT_VECTEUR_3D b,z;
        double *z=new double[n];

        for(int ip=0;ip<n;ip++)
                {
                        for(int iq=0;iq<n;iq++)
                                v[ip*n+iq]=0.0;
                        v[ip*n+ip]=1.0;
                        b[(unsigned int)ip]=d[(unsigned int)ip]=a[ip*n+ip];
                        z[(unsigned int)ip]=0.0;
                }
        nrot=0;
        for(int i=1;i<=ITMAX;i++)
                {
                        sm=0.0;
                        for(int ip=0;ip<n1;ip++)
                                for(int iq=ip+1;iq<n;iq++)
                                        sm+=fabs(a[ip*n+iq]);
                        if(sm==0.0)
                                return;
                        if(i<4)
                                tresh=0.2*sm*n2;
                        else
                                tresh=0.0;
                        for(int ip=0;ip<n1;ip++)
                                for(int iq=ip+1;iq<n;iq++)
                                        {
                                                g=100*fabs(a[ip*n+iq]);
                                                if(i>4&&(fabs(d[(unsigned int)ip])+g)==fabs(d[(unsigned int)ip])&&(fabs(d[(unsigned int)iq])+g)==fabs(d[(unsigned int)iq]))
                                                        a[ip*n+iq]=0.0;
                                                else if(fabs(a[ip*n+iq])>tresh)
                                                        {
                                                                h=d[(unsigned int)iq]-d[(unsigned int)ip];
                                                                if((fabs(h)+g)==fabs(h))
                                                                        t=a[ip*n+iq]/h;
                                                                else
                                                                        {
                                                                                theta=0.5*h/a[ip*n+iq];
                                                                                t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
                                                                                if(theta<0.0)t=-t;
                                                                        }
                                                                c=1.0/sqrt(1+t*t);
                                                                s=t*c;
                                                                tau=s/(1.0+c);
                                                                h=t*a[ip*n+iq];
                                                                z[(unsigned int)ip]-=h;
                                                                z[(unsigned int)iq]+=h;
                                                                d[(unsigned int)ip]-=h;
                                                                d[(unsigned int)iq]+=h;
                                                                a[ip*n+iq]=0.0;
                                                                ip1=ip-1;
                                                                iq1=iq-1;
                                                                for(int j=0;j<=ip1;j++)
                                                                        {
                                                                                g=a[j*n+ip];
                                                                                h=a[j*n+iq];
                                                                                a[j*n+ip]=g-s*(h+g*tau);
                                                                                a[j*n+iq]=h+s*(g-h*tau);
                                                                        }
                                                                for(int j=ip+1;j<=iq1;j++)
                                                                        {
                                                                                g=a[ip,j];
                                                                                h=a[j,iq];
                                                                                a[ip*n+j]=g-s*(h+g*tau);
                                                                                a[j*n+iq]=h+s*(g-h*tau);
                                                                        }       
                                                                for(int j=iq+1;j<n;j++)
                                                                        {
                                                                                g=a[ip*n+j];
                                                                                h=a[iq*n+j];
                                                                                a[ip*n+j]=g-s*(h+g*tau);
                                                                                a[iq*n+j]=h+s*(g-h*tau);
                                                                        }       
                                                                for(int j=0;j<n;j++)
                                                                        {
                                                                                g=v[j*n+ip];
                                                                                h=v[j*n+iq];
                                                                                v[j*n+ip]=g-s*(h+g*tau);
                                                                                v[j*n+iq]=h+s*(g-h*tau);
                                                                        }
                                                                ++nrot;
                                                        }
                                        }
                        for(int ip=0;ip<n;ip++)
                                {
                                        b[(unsigned int)ip]+=z[(unsigned int)ip];
                                        d[(unsigned int)ip]=b[(unsigned int)ip];
                                        z[(unsigned int)ip]=0.0;
                                }
                }


delete [] b;
delete [] z;
//#undef a
//#undef v
}


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

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:	58
Committed:	Fri Sep 28 21:00:16 2007 UTC (17 years, 7 months ago) by souaissa
Original Path:	*magic/lib/outil/outil/src/ot_mathematique.cpp*
File size:	22173 byte(s)
Error occurred while calculating annotation data.
Log Message:	vectorisation de geometrie pour comparaison