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 <fstream>
#include "ot_mathematique.h"


OT_HISTOGRAMME::OT_HISTOGRAMME(void)
{
  m_i_min=0;
  m_i_max=0;
  m_init=false;
  m_largeur_colonne=0;
}

OT_HISTOGRAMME::OT_HISTOGRAMME(double largeur_colonne):m_largeur_colonne(largeur_colonne)
{
  m_i_min=0;
  m_i_max=0;
  m_init=false;
}

OT_HISTOGRAMME::OT_HISTOGRAMME(OT_HISTOGRAMME& mdd)
{
  m_largeur_colonne=mdd.m_largeur_colonne;
  m_i_min=mdd.m_i_min;
  m_i_max=mdd.m_i_max;
  m_init=mdd.m_init;
  m_map_colonne=mdd.m_map_colonne;
}

OT_HISTOGRAMME::~OT_HISTOGRAMME(void)
{
  m_map_colonne.clear();
}

void OT_HISTOGRAMME::fixe_largeur_colonne(double largeur_colonne)
{
  m_largeur_colonne=largeur_colonne;
}

double OT_HISTOGRAMME::get_largeur_colonne(void)
{
  return m_largeur_colonne;
}

void OT_HISTOGRAMME::ajouter_valeur(double x,double val)
{
  long i = (long)(x/m_largeur_colonne);
  ajouter_valeur(i,val);
}

void OT_HISTOGRAMME::ajouter_valeur(long int i, double val)
{
  if(!m_init)
  {
    m_i_min=i;
    m_i_max=i;
    m_map_colonne.insert(std::pair<long,double>(i,val));
    m_init=true;
  }  
  else if(i>m_i_max)
  {
    for(long j=m_i_max+1;j<i;j++)    
      m_map_colonne.insert(std::pair<long,double>(j,0.0));
    m_i_max=i;
    m_map_colonne.insert(std::pair<long,double>(i,val));
  }
  else if(i<m_i_min)
  {
    for(long j=i+1;j<m_i_min;j++)    
      m_map_colonne.insert(std::pair<long,double>(j,0.0));
    m_i_min=i;
    m_map_colonne.insert(std::pair<long,double>(i,val));
  }
  else
  {
    m_map_colonne[i]+=val;
  }
}

double OT_HISTOGRAMME::get_valeur(double x)
{
  long i = (long)(x/m_largeur_colonne);
  std::map<long,double>::iterator it;
  it=m_map_colonne.find(i);
  if(it==m_map_colonne.end())
  {
    std::cerr << "*** ERREUR: OT_HISTOGRAMME::get hors de portee ***" << std::endl;
    return 1e300;
  }
  return it->second;  
}

int OT_HISTOGRAMME::get_premiere_valeur(std::map< long, double >::iterator& it,double &valeur)
{
  it=m_map_colonne.begin();
  if(it==m_map_colonne.end()) return 0;  
  valeur = it->second;
  return 1;
}

int OT_HISTOGRAMME::get_suivante_valeur(std::map< long, double >::iterator& it,double &valeur)
{
  it++;
  if(it==m_map_colonne.end()) return 0;  
  valeur = it->second;
  return 1;
}

int OT_HISTOGRAMME::get_premiere_valeur(std::map< long, double >::iterator& it,double &valeur,double &borne_min,double &borne_max)
{
  it=m_map_colonne.begin();
  if(it==m_map_colonne.end()) return 0;  
  valeur = it->second;
  borne_min=(it->first)*m_largeur_colonne;
  borne_max=(it->first)*m_largeur_colonne+m_largeur_colonne;
  return 1;
}

int OT_HISTOGRAMME::get_suivante_valeur(std::map< long, double >::iterator& it,double &valeur,double &borne_min,double &borne_max)
{
  it++;
  if(it==m_map_colonne.end()) return 0;  
  valeur = it->second;
  borne_min=(it->first)*m_largeur_colonne;
  borne_max=(it->first)*m_largeur_colonne+m_largeur_colonne;
  return 1;
}

int OT_HISTOGRAMME::get_premiere_valeur(std::map< long, double >::iterator& it, double &valeur, long int& i)
{
  it=m_map_colonne.begin();
  if(it==m_map_colonne.end()) return 0;  
  valeur = it->second;
  i = it->first;
  return 1;
}

int OT_HISTOGRAMME::get_suivante_valeur(std::map< long, double >::iterator& it, double &valeur, long int& i)
{
  it++;
  if(it==m_map_colonne.end()) return 0;  
  valeur = it->second;
  i = it->first;
  return 1;
}

long int OT_HISTOGRAMME::get_nb_colonne(void)
{
  return m_map_colonne.size();
}

double OT_HISTOGRAMME::get_x_min(void)
{
  return m_i_min*m_largeur_colonne;
}

double OT_HISTOGRAMME::get_x_max(void)
{
  return m_i_max*m_largeur_colonne+m_largeur_colonne;
}

OT_HISTOGRAMME& OT_HISTOGRAMME::operator=(OT_HISTOGRAMME& mdd)
{
  m_largeur_colonne=mdd.m_largeur_colonne;
  m_i_min=mdd.m_i_min;
  m_i_max=mdd.m_i_max;
  m_init=mdd.m_init;
  m_map_colonne=mdd.m_map_colonne;
}

void OT_HISTOGRAMME::enregistrer_bin(std::ofstream& ofstr)
{
  ofstr.write((char*)&m_init,sizeof(bool));
  ofstr.write((char*)&m_i_min,sizeof(long));
  ofstr.write((char*)&m_i_max,sizeof(long));
  ofstr.write((char*)&m_largeur_colonne,sizeof(double));
  long nb_colonne = get_nb_colonne();
  ofstr.write((char*)&nb_colonne,sizeof(long));
  std::map<long,double>::iterator it;
  for(it=m_map_colonne.begin();it!=m_map_colonne.end();it++)
  {
    long i = it->first;
    ofstr.write((char*)&i,sizeof(long));
    double val = it->second;
    ofstr.write((char*)&val,sizeof(double));
  }
}

void OT_HISTOGRAMME::ouvrir_bin(std::ifstream& ifstr)
{
  ifstr.read((char*)&m_init,sizeof(bool));
  ifstr.read((char*)&m_i_min,sizeof(long));
  ifstr.read((char*)&m_i_max,sizeof(long));
  ifstr.read((char*)&m_largeur_colonne,sizeof(double));
  long nb_colonne;
  ifstr.read((char*)&nb_colonne,sizeof(long));
  for(long j=0;j<nb_colonne;j++)
  {
    long k;
    ifstr.read((char*)&k,sizeof(long));
    double val;
    ifstr.read((char*)&val,sizeof(double));
    ajouter_valeur(k,val);
  }
}

void OT_HISTOGRAMME::exporter(std::ofstream& ofstrm)
{
  ofstrm << "x_min x_max x_milieu valeur" << std::endl;
  std::map<long,double>::iterator it_his;
  double val;
  double borne_min,borne_max;
  int k=get_premiere_valeur(it_his,val,borne_min,borne_max);
  while(k!=0)
  {
    double borne_mid = borne_min+m_largeur_colonne/2.;
    ofstrm << borne_min << " " << borne_max << " " << borne_mid << " " << val << std::endl;
    k=get_suivante_valeur(it_his,val,borne_min,borne_max);
  }
}


OT_QUATERNION::OT_QUATERNION(void)
{
  valeur[0]=0;
  valeur[1]=0;
  valeur[2]=0;
  valeur[3]=0;
}

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

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

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

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

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

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

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

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

double OT_QUATERNION::get_w(void)
{
  return valeur[3];
}

void OT_QUATERNION::change_w(double w)
{
  valeur[3]=w;
}

OT_QUATERNION& OT_QUATERNION::operator=(OT_QUATERNION& mdd)
{
  valeur[0]=mdd.valeur[0];
  valeur[1]=mdd.valeur[1];
  valeur[2]=mdd.valeur[2];
  valeur[3]=mdd.valeur[3];
  return *this;
}


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(double2 *xyz1,double2 *xyz2)
{
    valeur[0]=xyz2[0]-xyz1[0];
    valeur[1]=xyz2[1]-xyz1[1];
    valeur[2]=xyz2[2]-xyz1[2];
}

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_longueur(void) 
{
    return sqrt(valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
}

double2  OT_VECTEUR_3DD::get_longueur2(void)
{
    return (valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
}

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 << "[ "<<std::setw(20)<< __vec[0]<<std::setw(4)<<", "<<std::setw(20)<<__vec[1]<<", "<<std::setw(20)<<__vec[2]<<std::setw(4)<<", "<<std::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 << "[ "<<std::setw(20)<< __vec[0]<<std::setw(4)<<", "<<std::setw(20)<<__vec[1]<<", "<<std::setw(20)<<__vec[2]<<std::setw(4)<<", "<<std::setw(20)<<__vec[3]<<std::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);
}


bool OT_VECTEUR_4DD::est_nul(void)
{
double2 ZERO(0.);
if (valeur[0]==ZERO&&valeur[1]==ZERO &&valeur[2]==ZERO&&valeur[3]==ZERO)
        return true;
return false;
}

bool OT_VECTEUR_4DD::est_nul_3d(void)
{
double2 ZERO(0.);
if (valeur[0]==ZERO&&valeur[1]==ZERO &&valeur[2]==ZERO)
        return true;
return false;
}


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

int OPERATEUR::egal(double* xyz1, double* xyz2, double eps)
{
  if(fabs(xyz1[0]-xyz2[0])<eps)
    if(fabs(xyz1[1]-xyz2[1])<eps)
      if(fabs(xyz1[2]-xyz2[2])<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=std::max(d12,d13);
    dmax=std::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=std::max(dab,dac);
    dmax=std::max(dmax,dcb);
    dmax=std::max(dmax,ddb);
    dmax=std::max(dmax,dda);
    dmax=std::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;
}

int OPERATEUR::estdansletetra(double *xyz1,double *xyz2,double *xyz3,double *xyz4,double x,double y, double z)
{
OT_VECTEUR_3D v1(xyz2[0]-xyz1[0],xyz2[1]-xyz1[1],xyz2[2]-xyz1[2]);
OT_VECTEUR_3D v2(xyz3[0]-xyz1[0],xyz3[1]-xyz1[1],xyz3[2]-xyz1[2]);
OT_VECTEUR_3D v3(xyz4[0]-xyz1[0],xyz4[1]-xyz1[1],xyz4[2]-xyz1[2]);
OT_VECTEUR_3D v4(x-xyz1[0],y-xyz1[1],z-xyz1[2]);
OT_MATRICE_3D mat(v1,v2,v3);
OT_MATRICE_3D mat1(v4,v2,v3);
OT_MATRICE_3D mat2(v1,v4,v3);
OT_MATRICE_3D mat3(v1,v2,v4);
double det=mat.get_determinant();
double xsi=mat1.get_determinant()/det;
double eta=mat2.get_determinant()/det;
double dseta=mat3.get_determinant()/det;
int reponse1=1;
int reponse2=1;
double eps=0.0000001;
if (xsi<-eps) reponse1=0;
if (eta<-eps) reponse1=0;
if (dseta<-eps) reponse1=0;
if (xsi+eta+dseta>1.+eps) reponse1=0;
if (xsi<eps) reponse2=0;
if (eta<eps) reponse2=0;
if (dseta<eps) reponse2=0;
if (xsi+eta+dseta>1.-eps) reponse2=0;
int face1=0,face2=0,face3=0,face4=0;
if ((reponse1==1) && (reponse2==0))
    {
    if (dseta<eps) face1=1;  
    if (eta<eps) face2=1;  
    if (xsi<eps) face3=1;  
    if (xsi+eta+dseta>1.-eps) face4=1;
    }
return reponse1+2*reponse2+4*face1+8*face2+16*face3+32*face4;
}


int OPERATEUR::projeteestdansletriangle(double *xyz1,double *xyz2,double *xyz3,double x,double y, double z)
{
OT_VECTEUR_3D vec1(xyz1,xyz2);  
OT_VECTEUR_3D vec2(xyz1,xyz3);  
OT_VECTEUR_3D n=vec1&vec2;
n.norme();
double d=-(n.get_x()*xyz1[0]+n.get_y()*xyz1[1]+n.get_z()*xyz1[2]);
double t=-(n.get_x()*x+n.get_y()*y+n.get_z()*z+d);
double xx=x+t*n.get_x();
double yy=y+t*n.get_y();
double zz=z+t*n.get_z();
return estdansletriangle(xyz1,xyz2,xyz3,xx,yy,zz); 
}


int OPERATEUR::estdansletriangle(double *xyz1,double *xyz2,double *xyz3,double x,double y, double z)
{
OT_VECTEUR_3D vec1(xyz1,xyz2);  
OT_VECTEUR_3D vec2(xyz1,xyz3);  
OT_VECTEUR_3D vecb(x-xyz1[0],y-xyz1[1],z-xyz1[2]);      
OT_MATRICE_3D m(vec1,vec2,vecb);
double det=m.get_determinant();
double eps=1e-10*((vec1&vec2)*(vec1&vec2));
if (fabs(det)>eps) return 0;
double ab[3]={vec1.get_x(),vec1.get_y(),vec1.get_z()};
double ac[3]={vec2.get_x(),vec2.get_y(),vec2.get_z()};
double b[3]={vecb.get_x(),vecb.get_y(),vecb.get_z()};
det=ab[0]*ac[1]-ab[1]*ac[0];
int n1,n2;
if (fabs(det)>eps)
  {
  n1=0;
  n2=1;  
  }
else
  {
  det=ab[0]*ac[2]-ab[2]*ac[0];
  if (fabs(det)>eps)
    {
    n1=0;
    n2=2;  
    }
  else 
    {
    det=ab[1]*ac[2]-ab[2]*ac[1];
    n1=1;
    n2=2;
    }
  }
double xsi=(b[n1]*ac[n2]-b[n2]*ac[n1])/det;
double eta=(ab[n1]*b[n2]-ab[n2]*b[n1])/det;
int reponse1=1;
int reponse2=1;
eps=0.0000001;
if (xsi<-eps) reponse1=0;
if (eta<-eps) reponse1=0;
if (xsi+eta>1.+eps) reponse1=0;
if (xsi<eps) reponse2=0;
if (eta<eps) reponse2=0;
if (xsi+eta>1.-eps) reponse2=0;
int arete1=0,arete2=0,arete3=0;
if ((reponse1==1) && (reponse2==0))
    {
    if (eta<eps) arete1=1;  
    if (xsi<eps) arete3=1;  
    if (xsi+eta>1.-eps) arete2=1;
    }
return reponse1+2*reponse2+4*arete1+8*arete2+16*arete3;
}


int OPERATEUR::compare_etat_tetra(int etat,int valeur)
{
int faceinte=etat&1;
int inte=(etat&2)>>1;
int face1=(etat&4)>>2;
int face2=(etat&8)>>3;
int face3=(etat&16)>>4;
int face4=(etat&32)>>5;
if ((faceinte==1) && (valeur==INTERIEUR)) return 1;
if ((inte==1) && (valeur==INTERIEUR)) return 1;
if (/*(face==0) && */(inte==1) && (valeur==STRICTINTERIEUR)) return 1;
if ((faceinte==1) && (inte==0) && (valeur==SUR_FACE)) return 1;
if ((face1==1) && (valeur==FACE1)) return 1;
if ((face2==1) && (valeur==FACE2)) return 1;
if ((face3==1) && (valeur==FACE3)) return 1;
if ((face4==1) && (valeur==FACE4)) return 1;
if ((face1==1) && (face2==1) && (valeur==ARETE1)) return 1;
if ((face1==1) && (face3==1) && (valeur==ARETE2)) return 1;
if ((face2==1) && (face3==1) && (valeur==ARETE3)) return 1;
if ((face1==1) && (face4==1) && (valeur==ARETE4)) return 1;
if ((face2==1) && (face4==1) && (valeur==ARETE5)) return 1;
if ((face3==1) && (face4==1) && (valeur==ARETE6)) return 1;
if ((face1==1) && (face2==1) && (face3==1) && (valeur==SOMMET1)) return 1;
if ((face1==1) && (face2==1) && (face4==1) && (valeur==SOMMET2)) return 1;
if ((face1==1) && (face3==1) && (face4==1) && (valeur==SOMMET3)) return 1;
if ((face2==1) && (face3==1) && (face4==1) && (valeur==SOMMET4)) return 1;
return 0;
}

int OPERATEUR::compare_etat_triangle(int etat,int valeur)
{
int areteinte=etat&1;
int inte=(etat&2)>>1;
int arete1=(etat&4)>>2;
int arete2=(etat&8)>>3;
int arete3=(etat&16)>>4;
if ((areteinte==1) && (valeur==INTERIEUR)) return 1;
if ((inte==1) && (valeur==INTERIEUR)) return 1;
if ((inte==1) && (valeur==STRICTINTERIEUR)) return 1;
if ((areteinte==1) && (inte==0) && (valeur==SUR_ARETE)) return 1;
if ((arete1==1) && (valeur==ARETE1)) return 1;
if ((arete2==1) && (valeur==ARETE2)) return 1;
if ((arete3==1) && (valeur==ARETE3)) return 1;
if ((arete1==1) && (arete2==1) && (valeur==SOMMET2)) return 1;
if ((arete2==1) && (arete3==1) && (valeur==SOMMET3)) return 1;
if ((arete3==1) && (arete1==1) && (valeur==SOMMET1)) return 1;
return 0;
}

double OPERATEUR::taille_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=std::max(dab,dac);
  dmax=std::max(dmax,dcb);
  dmax=std::max(dmax,ddb);
  dmax=std::max(dmax,dda);
  dmax=std::max(dmax,ddc);
  return dmax;
}

double OPERATEUR::taille_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=std::max(d12,d13);
  dmax=std::max(dmax,d23);
  return dmax;
}


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:	951
Committed:	Fri Aug 10 15:17:17 2018 UTC (6 years, 9 months ago) by couturad
File size:	43125 byte(s)
Log Message:	-> Ajout de Project Chrono (voir CMakeLists.txt). -> Ajout d'un générateur de microstructure basé sur la dynamique des corps rigides (MSTRUCT_GENERATEUR_DCR). -> Ajout d'un opérateur de décallage de la topologie (MG_CG_OP_TRANSF_DECALLAGE). -> Retrait de «using namespace std» (conflit avec namespace chrono) et modification des fichiers affectés. -> Modification de mailleur2d.cpp afin d'enregistrer un fichier MAGiC (void.magic) lorsque le nombre d'itération dépasse la valeur maximale.