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=i-1;j>m_i_max;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
  {
    std::map<long,double>::iterator it;
    it=m_map_colonne.find(i);
    if(it==m_map_colonne.end())
    {
      m_map_colonne.insert(std::pair<long,double>(i,val));
    }
    else it->second+=val;
  }
}

int OT_HISTOGRAMME::get_valeur(double x, double& val)
{
  long i = (long)(x/m_largeur_colonne);
  return get_valeur(i,val);
}

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


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

OT_HISTOGRAMME & OT_HISTOGRAMME::operator=(const 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:	968
Committed:	Sun Sep 16 15:27:49 2018 UTC (6 years, 8 months ago) by couturad
File size:	43545 byte(s)
Log Message:	Ajout d'une condition de sortie et d'un renvoi d'erreur pour le mailleur FEM. Subdivision des fichiers mstruct_analyse.h/.cpp en sous fichiers pour une meilleure lisibilite. Ajout d'une analyse des modules d'elasticite. Ajout d'une analyse de l'energie. Reconfiguration du main de microstructure.exe (suppression d'actions obsolètes). Reconfiguration des fichiers generer_nb_ves, post_process.