Herezh_dev/herezh_pp/contact/Cercle.cc

417 lines
14 KiB
C++
Executable file

// 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 "Cercle.h"
#include "MathUtil.h"
#include "ParaGlob.h"
// CONSTRUCTEURS :
// par defaut
Cercle::Cercle () :
centre(),rayon(0.),N(NULL)
{ if (ParaGlob::Dimension()== 1)
{ cout << "\nErreur : la dimension est 1 et on veut utiliser un cercle ??? ";
cout << "\nCercle::Cercle ()" << endl;
Sortie(1);
};
};
// avec les datas
Cercle::Cercle ( const Coordonnee& B, const double& r , const Coordonnee* No):
centre(B),rayon(r),N(NULL)
{ if (ParaGlob::Dimension()== 1)
{ cout << "\nErreur : la dimension est 1 et on veut utiliser un cercle ??? ";
cout << "\nCercle::Cercle ()" << endl;
Sortie(1);
};
if (rayon < 0.)
{ cout << "\n erreur dans la construction d'un Cercle, le rayon " << rayon
<< " est negatif !! ";
Sortie(1);
};
if (ParaGlob::Dimension() == 3)
if (No != NULL)
{ N = new Coordonnee(*No);
double norme = N->Norme();
if (norme < ConstMath::petit)
{ cout << "\n erreur dans la construction d'un Cercle, le vecteur normal "
<< " a une norme (=" << norme <<") trop petite !! " ;
Sortie(1);
};
(*N) /= norme; // on normalise le vecteur
};
// sinon la normale n'existe pas, il faudra la définir
};
// avec la dimension
Cercle::Cercle (int dim):
centre(dim),rayon(0.),N(NULL)
{ if (ParaGlob::Dimension()== 1)
{ cout << "\nErreur : la dimension est 1 et on veut utiliser un cercle ??? ";
cout << "\nCercle::Cercle ()" << endl;
Sortie(1);
};
};
// de copie
Cercle::Cercle ( const Cercle& a):
centre(a.centre),rayon(a.rayon),N(NULL)
{ if (ParaGlob::Dimension() == 3)
if (a.N != NULL) // a priori il est normé
{ N = new Coordonnee(*(a.N));};
// sinon la normale n'existe pas, il faudra la définir
};
// DESTRUCTEUR :
Cercle::~Cercle ()
{ if (N != NULL)
{delete N;N=NULL;};
};
// surcharge des operator
Cercle& Cercle::operator = ( const Cercle & cer)
{ centre = cer.centre; rayon = cer.rayon;
if (ParaGlob::Dimension() == 3)
if (cer.N != NULL) // a priori la normale est normée
{ N = new Coordonnee(*(cer.N));};
// sinon la normale n'existe pas, il faudra la définir
return *this;
};
// METHODES PUBLIQUES :
// change le centre du Cercle
void Cercle::Change_centre( const Coordonnee& B)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != centre.Dimension())
{ cout << "\nErreur : les dimensions du centre et du vecteur ne sont pas identique !";
cout <<"\ndim B = " <<B.Dimension() <<", dim centre =" << centre.Dimension();
cout << "\nCercle::Change_centre( const Coordonnee& B)" << endl;
Sortie(1);
};
#endif
centre = B;
};
// change le rayon du Cercle
void Cercle::Change_rayon( const double & r)
{if (r < 0.)
{ cout << "\n erreur dans la construction d'un Cercle, le rayon " << r
<< " est negatif !! "
<< "\n Cercle::Change_rayon( const double & r) ";
};
rayon = r;
};
// change la normale au Cercle en 3D, sinon c'est inutile, warning
void Cercle::Change_Normale( const Coordonnee& No)
{ if (ParaGlob::Dimension() != 3)
{ cout << "\nwarning : on n'est pas en 3D, la normale n'est pas utile pour le cercle !";
cout << "\nCercle::Change_Normale( const Coordonnee& No)" << endl;
}
else
{
#ifdef MISE_AU_POINT
if (No.Dimension() != 3)
{ cout << "\nErreur : la dimension de la normale devrait etre 3 alors qu'elle est : "
<< No.Dimension() <<" ";
cout << "\nCercle::Change_Normale( const Coordonnee& No)" << endl;
Sortie(1);
};
#endif
if (N == NULL)
{ N = new Coordonnee (No);}
else
{ *N = No;};
double norme = N->Norme();
if (norme < ConstMath::petit)
{ cout << "\n erreur dans la construction d'un Cercle, le vecteur normal "
<< " a une norme (=" << norme <<") trop petite !! " ;
cout << "\nCercle::Change_Normale( const Coordonnee& No)" << endl;
Sortie(1);
};
(*N) /= norme; // on normalise le vecteur
};
};
// change toutes les donnees
void Cercle::change_donnees( const Coordonnee& B, const double& r,const Coordonnee* No)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != centre.Dimension())
{ cout << "\nErreur : les dimensions du centre et du vecteur ne sont pas identique !";
cout <<"\ndim B = " <<B.Dimension() <<", dim centre =" << centre.Dimension();
cout << "\nCercle::change_donnees( const Coordonnee& B, const double& r)" << endl;
Sortie(1);
};
#endif
centre = B;
if (r < 0.)
{ cout << "\n erreur dans la construction d'un Cercle, le rayon " << r
<< " est negatif !! "
<< "\n Cercle::change_donnees( ) ";
};
rayon = r;
// maintenant cas de la normale
if ((ParaGlob::Dimension() != 3) && (No != NULL))
{ cout << "\nwarning : on n'est pas en 3D, la normale n'est pas utile pour le cercle !";
cout << "\nCercle::change_donnees( const Coordonnee& No)" << endl;
}
else if ((ParaGlob::Dimension() == 3) && (No == NULL))
{ cout << "\nwarning : la normale au cercle n'est pas defini !";
cout << "\nCercle::change_donnees()" << endl;
if (N != NULL)
delete N;
}
else if (ParaGlob::Dimension() == 3) // cas de la dimension 3 et No non null
{
#ifdef MISE_AU_POINT
if (No->Dimension() != 3)
{ cout << "\nErreur : la dimension de la normale devrait etre 3 alors qu'elle est : "
<<No->Dimension() <<" ";
cout << "\nCercle::change_donnees()" << endl;
Sortie(1);
};
#endif
if (N == NULL)
{ N = new Coordonnee (*No);}
else
{ *N = *No;};
double norme = N->Norme();
if (norme < ConstMath::petit)
{ cout << "\n erreur dans la construction d'un Cercle, le vecteur normal "
<< " a une norme (=" << norme <<") trop petite !! " ;
cout << "\nCercle::change_donnees(.." << endl;
Sortie(1);
};
(*N) /= norme; // on normalise le vecteur
};
};
// calcul les deux intercections M1 et M2 d'une droite avec le Cercle,
// 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 deux points d'intercection
int Cercle::Intersection( const Droite & dr,Coordonnee& M1, Coordonnee& M2)
{switch (ParaGlob::Dimension())
{case 2:
{// C le centre du cercle, H la projection de C sur la droite
// M1 et M2 les deux intersections potentielles
const Coordonnee& A = dr.PointDroite(); // pour simplifier
const Coordonnee& U = dr.VecDroite();
Coordonnee CA = A - centre;
double ah = - CA * U;
Coordonnee CH = CA + ah * U;
double ch_2 = CH * CH;
if (sqrt(ch_2) > rayon)
return 0; // cas pas d'intersection
// sinon
double hm = sqrt(rayon*rayon - ch_2);
M1 = centre + CH + hm * U;
M2 = centre + CH - hm * U;
return 1;
break;
}
case 3:
{
#ifdef MISE_AU_POINT
if (N==NULL)
{ cout << "\nErreur : la normale au cercle n'est pas defini, "
<< " on ne peut pas calculer la distance au cercle ";
cout << "\nCercle::Distance_au_Cercle()" << endl;
Sortie(1);
};
#endif
// on calcul le projeté sur le plan d'un point de la droite
const Coordonnee& A = dr.PointDroite(); // pour simplifier
const Coordonnee& U = dr.VecDroite();
Coordonnee CA = A-centre;
double ha = CA*(*N); // distance du point au plan du cercle
if (ha > ConstMath::pasmalpetit)
return 0; // cas pas d'intersection
// sinon le point appartient au plan: même procédure qu'en 2D
// C le centre du cercle, H la projection de C sur la droite
// M1 et M2 les deux intersections potentielles
double ah = - CA * U;
Coordonnee CH = CA + ah * U;
double ch_2 = CH * CH;
if (sqrt(ch_2) > rayon)
return 0; // cas pas d'intersection
// sinon
double hm = sqrt(rayon*rayon - ch_2);
M1 = centre + CH + hm * U;
M2 = centre + CH - hm * U;
return 1;
break;
}
case 1:
{ cout << "\nErreur : la dimension est 1 et on utilise un cercle ??? ";
cout << "\nCercle::Intersection(.." << endl;
Sortie(1);
};
};
// normalement on n'arrive jamais ici
return -1;
};
// calcul la distance d'un point au Cercle
double Cercle::Distance_au_Cercle(const Coordonnee& M) const
{ switch (ParaGlob::Dimension())
{case 2:
return Dabs((centre - M).Norme()-rayon);break;
case 3:
{
#ifdef MISE_AU_POINT
if (N==NULL)
{ cout << "\nErreur : la normale au cercle n'est pas defini, "
<< " on ne peut pas calculer la distance au cercle ";
cout << "\nCercle::Distance_au_Cercle()" << endl;
Sortie(1);
};
#endif
// on calcul le projeté sur le plan
Coordonnee CM = M-centre;
double cm = (CM*(*N));
Coordonnee CH = CM - cm * (*N);
// calcul du point intersection CH avec le cercle -> E
double ch = CH.Norme();
if (ch < ConstMath::trespetit)
{ // cela veut dire que le point M est sur l'axe du cercle
return (sqrt(cm*cm + rayon*rayon));
}
else
{ Coordonnee CE = CH * (rayon/ch);
return (CM-CE).Norme();
};
break;
}
case 1:
{ cout << "\nErreur : la dimension est 1 et on utilise un cercle ??? ";
cout << "\nCercle::Distance_au_Cercle()" << endl;
Sortie(1);
};
};
// normalement on n'arrive jamais ici
return ConstMath::tresgrand;
};
// calcul du projeté d'un point sur le plan
Coordonnee Cercle::Projete(const Coordonnee& M) const
{ switch (ParaGlob::Dimension())
{case 2: // le projeté est directement M
return M;break;
case 3:
{
#ifdef MISE_AU_POINT
if (N==NULL)
{ cout << "\nErreur : la normale au cercle n'est pas defini, "
<< " on ne peut pas calculer la distance au cercle ";
cout << "\nCercle::Projete(.." << endl;
Sortie(1);
};
#endif
// on calcul le projeté sur le plan
return (M - ((M-centre)* (*N))* (*N));
break;
}
case 1:
{ cout << "\nErreur : la dimension est 1 et on utilise un cercle ??? ";
cout << "\nCercle::Projete(.." << endl;
Sortie(1);
};
};
// normalement on n'arrive jamais ici
Coordonnee toto; // on fait en deux fois pour la version linux
return (toto);
};
// fonction valable uniquement en 2D, sinon erreur
// ramène true si le point est à l'intérieur du cercle, false sinon
bool Cercle::Dedans(const Coordonnee& M)const
{ if (ParaGlob::Dimension() == 2)
{ return ((M-centre).Norme() < rayon);}
else
{ cout << "\n erreur, la fonction Cercle::Dedans(... n'est utilisable "
<< " qu'en dimension 2D, alors que l'on est en " << ParaGlob::Dimension() << " ";
Sortie(1);
return false; // pour éviter un warning
};
};
// surcharge de l'operateur de lecture
istream & operator >> (istream & entree, Cercle & cer)
{ // vérification du type
string nom;
entree >> nom;
#ifdef MISE_AU_POINT
if (nom != "_cercle_")
{ cout << "\nErreur, en lecture d'une instance cercle "
<< " on attendait _cercle_ et on a lue: " << nom ;
cout << "istream & operator >> (istream & entree, Cercle & cer)\n";
Sortie(1);
};
#endif
// puis lecture des différents éléments
entree >> nom >> cer.centre >> nom >> cer.rayon;
bool existe_N;
entree >> existe_N;
if (existe_N)
{ if (cer.N == NULL)
cer.N = new Coordonnee(ParaGlob::Dimension());
entree >> (*(cer.N));
double norme = cer.N->Norme();
if (norme < ConstMath::petit)
{ cout << "\n erreur dans la construction d'un Cercle, le vecteur normal "
<< " a une norme (=" << norme <<") trop petite !! " ;
cout << "\nCercle::change_donnees(.." << endl;
Sortie(1);
};
(*cer.N) /= norme; // on normalise le vecteur
}
else
{ if (cer.N != NULL)
delete cer.N;
};
return entree;
};
// surcharge de l'operateur d'ecriture
ostream & operator << ( ostream & sort,const Cercle & cer)
{ // tout d'abord un indicateur donnant le type
sort << " _cercle_ " ;
// puis les différents éléments
sort << "\n C= " << cer.centre << " r= " << cer.rayon << " ";
if (cer.N == NULL)
{ sort << 0 ;}
else
{ sort << 1 << (*(cer.N)) << " " ;}
return sort;
};