rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/herezh_pp/contact/Droite.cc

346 lines
13 KiB
C++
Executable file
Raw Permalink Blame History

//#include "Debug.h"
// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2021 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
#include "Droite.h"
#include "ConstMath.h"
#include "MathUtil.h"
#include "Util.h"
//----------------------------------------------------------------
// def des donnees commune a toutes les droites
//----------------------------------------------------------------
int Droite::alea = 0;
// CONSTRUCTEURS :
// par defaut
// par defaut définit une droit // à x et passant par 0, et de dimension celle de l'espace de travail
Droite::Droite () :
A(ParaGlob::Dimension()), U(ParaGlob::Dimension())
{ // on met une valeur 1 par défaut sur U
U(1)=1.;
};
// avec les datas
Droite::Droite ( const Coordonnee& B, const Coordonnee& vec) :
A(B), U(vec)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != vec.Dimension())
{ cout << "\nErreur : les dimensions du point et du vecteur ne sont pas identique !";
cout <<"\ndim point = " <<B.Dimension() <<", dim vecteur =" << vec.Dimension();
cout << "\nDroite::Droite (Coordonnee& B,Coordonnee& vec)" << endl;
Sortie(1);
};
#endif
double d = U.Norme();
if (d <= ConstMath::trespetit)
{ cout << "\nErreur : la norme du vecteur est trop petite !";
cout <<"\nnorme = " << d;
cout << "\nDroite::Droite (Coordonnee& B,Coordonnee& vec)" << endl;
Sortie(1);
};
U /= d; // U est ainsi norme
};
// avec la dimension
// def par défaut d'une droite // à x
Droite::Droite (int dim) :
A(dim), U(dim)
{ // on met une valeur 1 par défaut sur U
U(1)=1.;
};
// de copie
Droite::Droite ( const Droite& a) :
A(a.A), U(a.U)
{};
// DESTRUCTEUR :
Droite::~Droite () {};
// surcharge des operator
Droite& Droite::operator = ( const Droite & P)
{ this->A = P.A; this->U = P.U;
return *this;
};
// METHODES PUBLIQUES :
// change le point de ref de la droite
void Droite::Change_ptref( const Coordonnee& B)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != U.Dimension())
{ cout << "\nErreur : les dimensions du point et du vecteur ne sont pas identique !";
cout <<"\ndim point = " <<B.Dimension() <<", dim vecteur =" << U.Dimension();
cout << "\nDroite::Change_ptref(Coordonnee& B)" << endl;
Sortie(1);
};
#endif
A = B;
};
// change le vecteur de la droite
void Droite::Change_vect( const Coordonnee& vec)
{
#ifdef MISE_AU_POINT
if (A.Dimension() != vec.Dimension())
{ cout << "\nErreur : les dimensions du point et du vecteur ne sont pas identique !";
cout <<"\ndim point = " <<A.Dimension() <<", dim vecteur =" << vec.Dimension();
cout << "\nDroite::Change_vect(Coordonnee& vec)" << endl;
Sortie(1);
};
#endif
double d = vec.Norme();
if (d <= ConstMath::trespetit)
{ cout << "\nErreur : la norme du vecteur est trop petite !";
cout <<"\nnorme = " << d;
cout << "\nDroite::Droite (Coordonnee& B,Coordonnee& vec)" << endl;
Sortie(1);
};
U = vec / d; // U est ainsi norme
};
// change toutes les donnees
void Droite::change_donnees( const Coordonnee& B, const Coordonnee& vec)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != vec.Dimension())
{ cout << "\nErreur : les dimensions du point et du vecteur ne sont pas identique !";
cout <<"\ndim point = " <<B.Dimension() <<", dim vecteur =" << vec.Dimension();
cout << "\nDroite::Change_vect(Coordonnee& vec)" << endl;
Sortie(1);
};
#endif
A = B;
double d = vec.Norme();
if (d <= ConstMath::trespetit)
{ cout << "\nErreur : la norme du vecteur est trop petite !";
cout <<"\nnorme = " << d;
cout << "\nDroite::Droite (Coordonnee& B,Coordonnee& vec)" << endl;
Sortie(1);
};
U = vec / d; // U est ainsi norme
};
// calcul l'intersection M de la droite avec une autre droite, ramene 0 s'il n'y
// a pas d'intersection, ramene -1 si l'intersection ne peut pas etre calculee
// et 1 s'il y a un point d'intersection
int Droite::Intersection( const Droite & D,Coordonnee& M)
{ const Coordonnee & B = D.PointDroite(); // par commodite d'ecriture
const Coordonnee & V = D.VecDroite(); // ""
Coordonnee AB = B - A;double Nab = AB.Norme();
// test en fonction de la dimension (a priori 2 ou 3 !!)
switch (A.Dimension())
{case 3:
{ Coordonnee W = Util::ProdVec_coor(U,V);
double d = W.Norme();
// if (Dabs(d) <= ConstMath::pasmalpetit)
if (Dabs(d) <= ConstMath::petit)
{ // cas de droite // ou confondues
// if (Dabs(Nab) <= ConstMath::pasmalpetit)
if (Dabs(Nab) <= ConstMath::petit)
// les droites sont confondues et ont meme point de reference
return -1; // il n'y a pas d'Intersection identifiee
AB /= Nab; // norme AB
Coordonnee V1 = Util::ProdVec_coor(AB,V);double Nv1= V1.Norme();
if (Dabs(Nv1) <= ConstMath::petit)
// produit vectoriel nul: les droites sont confondues, tous les points conviennent
// mais on ne peut pas calculer un seul point d'intersection, on choisit par exemple le point de ref:
{ M= B; return -1;}
else
// les droites sont distinctes: // mais avec des points de reference differents
{return -1;} // il n'y a pas d'interection identifiee
}
else
// on calcul la distance entre les droites
{ W /= d;
double h = AB * W; // h = la distance
if (Dabs(h) >= ConstMath::petit) //pasmalpetit)
// trop de distance entre les deux droites -> pas d'Intersection
return 0;
}
}
break;
case 2:
{ // on regarde si les droites ne sont pas //
if ( Dabs(1.-Dabs(U*V)) < ConstMath::petit) //pasmalpetit)
// soit les droites sont distinctes soit elles sont confondues
// mais dans tous les cas, on ne peut pas calculer un seul point d'intersection
return 0;
}
break;
case 1:
{ // en 1D on considère que l'on ne peut pas calculer un seul point d'intersection
return 0;
}
break;
default:
cout << "\n erreur, on ne devrait pas passer par ici !!! "
<< "\n *** Droite::Intersection( const Droite & D,Coordonnee& M) " << endl;
Sortie(1);
};
// cas ou l'Intersection existe
// détail du calcul du point M = intersection
// on a: vec(AB) = vec(AM) + vec(MB) = alpha U + beta V
// d'où : (1): AB*U = alpha - beta U*V et (2): AB*V = alpha U*V - beta
// on fait (1) * (-U*V) + (2) -> - (AB*U) (U*V) + (AB*V) = beta ((U*V)^2 -1)
// et BM = beta V
double uv = U * V;
double beta = -((AB * U) * uv - (AB * V)) / (uv * uv -1.);
M = B + beta * V;
////--- debug
//cout << "\n debug Droite::Intersection( ";
//A.Affiche(cout, 10); M.Affiche(cout, 10); cout << "\n" << endl;
//cout << "\n uv= " << uv << ", V: " << V(1) << " " << V(2) << " V= " << sqrt(V(1)*V(1)+V(2)*V(2))
// << " B: " << B(1) << " " << B(2) << ", AB: " << AB(1) << " " << AB(2) << ", beta= " << beta;
//double diff = (B-M)*(A-M); cout << " (B-M)*(A-M) "<< diff;
//
//cout << "\n A: " << A(1) << " " << A(2) << ", U: " << U(1) << " " << U(2)
// << ", M: " << M(1) << " " << M(2) << endl;
////--- fin debug
return 1;
};
// calcul si un point donné appartient ou non à la droite (à la précision donnée)
bool Droite::Appartient_a_la_droite(const Coordonnee& B)
{ // le point B appartient à la droite lorsque AB est // à U donc que le produit
// vectoriel de ces deux vecteurs est nul en 3D, en 2D c'est le déterminant et en 1D c'est toujours ok
Coordonnee AB = B-A;
switch (AB.Dimension())
{ case 1 : // la collinéarité est obligatoire ici
return true; break;
case 2 : // examen du déterminant
// if (Abs(Util::Determinant(U,AB)) <= ConstMath::pasmalpetit) return true; else return false;
if (Abs(Util::Determinant(U,AB)) <= ConstMath::petit) return true; else return false;
break;
case 3 :
// if((Util::ProdVec_coor(U,AB)).Norme() <= ConstMath::pasmalpetit) return true; else return false;
if((Util::ProdVec_coor(U,AB)).Norme() <= ConstMath::petit) return true; else return false;
break;
default :
cout << "\nErreur : la dimension de AB est nulle !, suite des calculs impossible\n";
cout << "Droite::Appartient_a_la_droite... \n";
Sortie(1);
};
// ne doit jamais passer ici mais pour taire le compilo ....
Sortie(1);
return false;
};
// ramene une normale, en 2D il y a une seule solution, en 3D a chaque appel on ramene
// une nouvelle normale calculée aleatoirement
Coordonnee Droite::UneNormale()
{ // on considere la dimension
switch (A.Dimension())
{ case 2:
// cas ou l'on est en dimension 2
{ Coordonnee ve(2);
ve(1) = -U(2);
ve(2) = U(1);
return ve;
}
break;
case 3: // cas où l'on est en dimension 3
{ Coordonnee ve(3);
Coordonnee v1(3),v2(3);
// on cherche un vecteur du plan normal a la droite : v1
v2(1) = 1.;
v1 = Util::ProdVec_coor(U,v2);
double r = v1.Norme();
if (r <= ConstMath::petit)
// cas U // v2
{ v2(2) = 1.; v2(1) = 0;
v1 = Util::ProdVec_coor(U,v2);
};
// maintenant le second
v2 = Util::ProdVec_coor(U,v1);
// on genere un angle aleatoirement -> un cos
double cos = ((double)rand())/RAND_MAX;
// et le vecteur correspondant
ve = cos * v1 + (1.-cos) * v2;
// update le nombre aleatoire
alea++; srand((unsigned) alea);
return ve;
}
break;
}
};
// calcul la distance d'un point à la droite
double Droite::Distance_a_la_droite(const Coordonnee& M) const
{ int dima = ParaGlob::Dimension();
Coordonnee AM(M-A);
Coordonnee HM(AM-(AM * U)*U);
return HM.Norme();
};
// fonction valable uniquement en 2D !!!, sinon erreur
// ramène true si les deux points sont du même coté de la droite, false sinon
bool Droite::DuMemeCote(const Coordonnee& M1, const Coordonnee& M2) const
{ if (ParaGlob::Dimension() == 2)
{
// erreur 12 avril 2011: correction, il faut utiliser le déterminant et non le produit vectoriel
// return ((Util::ProdVec_coor(U,(M1-A)) * Util::ProdVec_coor(U,(M2-A))) >= 0.);
return ((Util::Determinant(U,(M1-A)) * Util::Determinant(U,(M2-A))) >= 0.);
}
else
{ cout << "\n erreur, la fonction Droite::DuMemeCote(.., n'est pas utilisable pour une dimension "
<< " differente de 2D ";
Sortie(1);
return false; // pour éviter un warning
}
};
// surcharge de l'operateur de lecture
istream & operator >> (istream & entree, Droite & dr)
{ // vérification du type
string nom;
entree >> nom;
#ifdef MISE_AU_POINT
if (nom != "_droite_")
{ cout << "\nErreur, en lecture d'une instance Droite "
<< " on attendait _droite_ et on a lue: " << nom ;
cout << "istream & operator >> (istream & entree, Droite & dr)\n";
Sortie(1);
};
#endif
// puis lecture des différents éléments
entree >> nom >> dr.A >> nom >> dr.U;
return entree;
};
// surcharge de l'operateur d'ecriture
ostream & operator << ( ostream & sort,const Droite & dr)
{ // tout d'abord un indicateur donnant le type
sort << " _droite_ " ;
// puis les différents éléments
sort << "\n A= " << dr.A << " U= " << dr.U << " ";
return sort;
};
Reference in a new issue View git blame Copy permalink