rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity

ajout des fichiers du répertoire contact

Browse source
This commit is contained in:
Gérard Rio 2021-09-24 08:24:03 +02:00
parent d367a29e8c
commit 5da18302a8
17 changed files with 10910 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

416
contact/Cercle.cc Normal file
View file

@ -0,0 +1,416 @@
// 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;
};

137
contact/Cercle.h Normal file
View file

@ -0,0 +1,137 @@
// 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/>.
/************************************************************************
* DATE: 19/01/2007 *
* $ *
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
* $ *
* PROJET: Herezh++ *
* $ *
************************************************************************
* BUT: Def géométrie d'un cercle: un centre , un rayon *
* et une normale au plan du cercle dans le cas 3D *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
* VERIFICATION: *
* *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* ! ! ! ! *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* MODIFICATIONS: *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* $ *
************************************************************************/
#ifndef CERCLEE_H
#define CERCLEE_H
#include "Droite.h"
#include "Coordonnee.h"
/// @addtogroup Groupe_sur_les_contacts
/// @{
///
class Cercle
{ // surcharge de l'operator de lecture
friend istream & operator >> (istream &, Cercle &);
// surcharge de l'operator d'ecriture
friend ostream & operator << (ostream &, const Cercle &);
public :
// CONSTRUCTEURS :
// par defaut
Cercle ();
// avec les datas
Cercle ( const Coordonnee& B, const double& r, const Coordonnee* N);
// avec la dimension
Cercle (int dim);
// de copie
Cercle ( const Cercle& a);
// DESTRUCTEUR :
~Cercle ();
// surcharge des operator
Cercle& operator = ( const Cercle & P);
// METHODES PUBLIQUES :
// retourne le centre du Cercle
inline const Coordonnee& CentreCercle() const { return centre;};
// retourne le rayon du Cercle
inline double RayonCercle() const { return rayon;};
// retourne la normale au cercle en 3D sinon NULL
inline const Coordonnee* NormaleCercle() const { return N;};
// change le centre du Cercle
void Change_centre( const Coordonnee& B);
// change le rayon du Cercle
void Change_rayon( const double & r);
// change la normale au Cercle en 3D, sinon c'est inutile, warning
void Change_Normale( const Coordonnee& N);
// change toutes les donnees
// N n'est utilisé qu'en 3D, en 2D on met NULL
void change_donnees( const Coordonnee& B, const double& r,const Coordonnee* N);
// 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'intercection
// ne peut pas etre calculee
// et 1 s'il y a deux points d'intercection
int Intersection( const Droite & D,Coordonnee& M1, Coordonnee& M2);
// calcul la distance d'un point au Cercle
double Distance_au_Cercle(const Coordonnee& M) const;
// calcul du projeté d'un point sur le plan
Coordonnee Projete(const Coordonnee& M) const;
// indique si oui ou non le projeté du point est dans le cercle
bool ProjeteDedans(const Coordonnee& M) const
{ return ((Projete(M)-centre).Carre() < rayon*rayon);};
// fonction valable uniquement en 2D, sinon erreur
// ramène true si le point est à l'intérieur du cercle, false sinon
bool Dedans(const Coordonnee& M)const ;
protected :
// VARIABLES PROTEGEES :
Coordonnee centre; // centre du Cercle
double rayon; // rayon du Cercle
Coordonnee* N; // normal au cercle : utilisé qu'en 3D
// la grandeur stockée est normée
// METHODES PROTEGEES :
};
/// @} // end of group
#endif

189
contact/Cylindre.cc Normal file
View file

@ -0,0 +1,189 @@
// 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 "Cylindre.h"
#include "MathUtil.h"
#include "ParaGlob.h"
// CONSTRUCTEURS :
// par defaut
Cylindre::Cylindre () :
cercle(3)
{ if (ParaGlob::Dimension() != 3)
{ cout << "\nErreur : la dimension est differente de 3 et on veut utiliser un Cylindre ??? ";
cout << "\nCylindre::Cylindre ()" << endl;
Sortie(1);
};
};
// avec les datas
Cylindre::Cylindre ( const Cercle& cer):
cercle(cer)
{ if (ParaGlob::Dimension() != 3)
{ cout << "\nErreur : la dimension est differente de 3 et on veut utiliser un Cylindre ??? ";
cout << "\nCylindre::Cylindre (.." << endl;
Sortie(1);
};
};
// avec la dimension
Cylindre::Cylindre (int dim):
cercle(dim)
{ if (ParaGlob::Dimension() != 3)
{ cout << "\nErreur : la dimension est differente de 3 et on veut utiliser un Cylindre ??? ";
cout << "\nCylindre::Cylindre ()" << endl;
Sortie(1);
};
if (dim != 3)
{ cout << "\nErreur : la dimension demandee est differente de 3 et on veut utiliser un Cylindre ??? ";
cout << "\nCylindre::Cylindre ()" << endl;
Sortie(1);
};
};
// de copie
Cylindre::Cylindre ( const Cylindre& a):
cercle(a.cercle)
{ };
// DESTRUCTEUR :
Cylindre::~Cylindre ()
{ };
// surcharge des operator
Cylindre& Cylindre::operator = ( const Cylindre & cer)
{ cercle = cer.cercle;
return *this;
};
// METHODES PUBLIQUES :
// change le cercle du Cylindre
void Cylindre::Change_cercle( const Cercle& cer)
{
#ifdef MISE_AU_POINT
if (cer.CentreCercle().Dimension() != cercle.CentreCercle().Dimension())
{ cout << "\nErreur : les dimensions du cercle demande n'est pas identique a celle existante !";
cout <<"\ndim nouveau = " <<cer.CentreCercle().Dimension()
<<", dim ancien =" << cercle.CentreCercle().Dimension();
cout << "\nSphere::Cylindre::Change_cercle( const Cercle& cer)" << endl;
Sortie(1);
};
#endif
cercle = cer;
};
// change toutes les donnees
void Cylindre::change_donnees( const Cercle& cer)
{
#ifdef MISE_AU_POINT
if (cer.CentreCercle().Dimension() != cercle.CentreCercle().Dimension())
{ cout << "\nErreur : les dimensions du cercle demande n'est pas identique a celle existante !";
cout <<"\ndim nouveau = " <<cer.CentreCercle().Dimension()
<<", dim ancien =" << cercle.CentreCercle().Dimension();
cout << "\nCylindre::change_donnees( const Cercle& cer,const double& haut)" << endl;
Sortie(1);
};
#endif
cercle = cer;
};
// calcul la distance d'un point au Cylindre
double Cylindre::Distance_au_Cylindre(const Coordonnee& M) const
{ // on commence par calculer la projection du point sur l'axe
Coordonnee CM = M-cercle.CentreCercle();
const Coordonnee& U = (*(cercle.NormaleCercle()));
double ch = CM * U;
Coordonnee CB = ch * U; // B = projection de M sur l'axe
// def de la distance
return Dabs((CM-CB).Norme() - cercle.RayonCercle());
};
// ramène true si le point est à l'intérieur du cylindre, false sinon
bool Cylindre::Dedans(const Coordonnee& M)const
{ // on commence par calculer la projection du point sur l'axe
Coordonnee CM = M-cercle.CentreCercle();
const Coordonnee& U = (*(cercle.NormaleCercle()));
double ch = CM * U;
Coordonnee CB = ch * U; // B = projection de M sur l'axe
// calcul de BM
Coordonnee BM =(CM-CB);
double dist_axe = BM.Norme();
return (dist_axe < cercle.RayonCercle()) ;
};
// projection d'un point M sur la parois du cylindre
// dans le cas où M appartient à l'axe du cylindre, la projection n'est pas
// possible, dans ce cas projection_ok = false en retour, sinon true
Coordonnee Cylindre::Projete(const Coordonnee& M,bool& projection_ok) const
{ // on commence par calculer la projection du point sur l'axe
Coordonnee CM = M-cercle.CentreCercle();
const Coordonnee& U = (*(cercle.NormaleCercle()));
double ch = CM * U;
Coordonnee CB = ch * U; // B = projection de M sur l'axe
// calcul du projeté
Coordonnee BM =(CM-CB);
double dist_axe = BM.Norme();
if (dist_axe < ConstMath::pasmalpetit)
{projection_ok = false;// on signale que la projection n'est pas possible
Coordonnee toto;
return (toto);// on ramène un point par défaut
}
else
{projection_ok = true;
Coordonnee P = cercle.CentreCercle() + BM * (cercle.RayonCercle() / dist_axe );
return P;
};
};
// surcharge de l'operateur de lecture
istream & operator >> (istream & entree, Cylindre & cer)
{ // vérification du type
string nom;
entree >> nom;
#ifdef MISE_AU_POINT
if (nom != "_Cylindre_")
{ cout << "\nErreur, en lecture d'une instance Cylindre "
<< " on attendait _Cylindre_ et on a lue: " << nom ;
cout << "istream & operator >> (istream & entree, Cylindre & cer)\n";
Sortie(1);
};
#endif
// puis lecture des différents éléments
entree >> nom >> cer.cercle ;
return entree;
};
// surcharge de l'operateur d'ecriture
ostream & operator << ( ostream & sort,const Cylindre & cer)
{ // tout d'abord un indicateur donnant le type
sort << " _Cylindre_ " ;
// puis les différents éléments
sort << "\n Cercle= " << cer.cercle << " ";
return sort;
};

121
contact/Cylindre.h Normal file
View file

@ -0,0 +1,121 @@
// 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/>.
/************************************************************************
* DATE: 19/01/2007 *
* $ *
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
* $ *
* PROJET: Herezh++ *
* $ *
************************************************************************
* BUT: Def géométrie d'un cylindre: un cercle 3D, donc avec une *
* normale. Le cylindre n'est pas limité en hauteur. *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
* VERIFICATION: *
* *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* ! ! ! ! *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* MODIFICATIONS: *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* $ *
************************************************************************/
#ifndef CYLINDREE_H
#define CYLINDREE_H
#include "Droite.h"
#include "Coordonnee.h"
#include "Cercle.h"
/// @addtogroup Groupe_sur_les_contacts
/// @{
///
class Cylindre
{ // surcharge de l'operator de lecture
friend istream & operator >> (istream &, Cylindre &);
// surcharge de l'operator d'ecriture
friend ostream & operator << (ostream &, const Cylindre &);
public :
// CONSTRUCTEURS :
// par defaut
Cylindre ();
// avec les datas
Cylindre ( const Cercle& B);
// avec la dimension
Cylindre (int dim);
// de copie
Cylindre ( const Cylindre& a);
// DESTRUCTEUR :
~Cylindre ();
// surcharge des operator
Cylindre& operator = ( const Cylindre & P);
// METHODES PUBLIQUES :
// retourne le cercle d'embase du Cylindre
inline const Cercle& CercleCylindre() const { return cercle;};
// change le centre du Cylindre
void Change_cercle( const Cercle& B);
// change toutes les donnees
void change_donnees( const Cercle& cer);
// calcul les deux intercections M1 et M2 d'une droite avec le Cylindre,
// ramene 0 s'il n'y a pas d'intersection, ramene -1 si l'intercection
// ne peut pas etre calculee
// et 1 s'il y a deux points d'intercection
int Intersection( const Droite & D,Coordonnee& M1, Coordonnee& M2);
// calcul la distance d'un point au Cylindre
double Distance_au_Cylindre(const Coordonnee& M) const;
// ramène true si le point est à l'intérieur du cylindre, false sinon
bool Dedans(const Coordonnee& M)const ;
// projection d'un point M sur la parois du cylindre
// dans le cas où M appartient à l'axe du cylindre, la projection n'est pas
// possible, dans ce cas projection_ok = false en retour, sinon true
Coordonnee Projete(const Coordonnee& M,bool& projection_ok) const;
protected :
// VARIABLES PROTEGEES :
Cercle cercle; // cercle du Cylindre
// METHODES PROTEGEES :
};
/// @} // end of group
#endif

345
contact/Droite.cc Normal file
View file

@ -0,0 +1,345 @@
//#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;
};

135
contact/Droite.h Normal file
View file

@ -0,0 +1,135 @@
// 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/>.
/************************************************************************
* DATE: 23/01/97 *
* $ *
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
* $ *
* PROJET: Herezh++ *
* $ *
************************************************************************
* BUT: Def de la geometrie d'une droite: un point et un vecteur *
* directeur. *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
* VERIFICATION: *
* *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* ! ! ! ! *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* MODIFICATIONS: *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* $ *
************************************************************************/
#ifndef DROITE_H
#define DROITE_H
#include "Coordonnee.h"
/// @addtogroup Groupe_sur_les_contacts
/// @{
///
class Droite
{ // surcharge de l'operator de lecture
friend istream & operator >> (istream &, Droite &);
// surcharge de l'operator d'ecriture
friend ostream & operator << (ostream &, const Droite &);
public :
// CONSTRUCTEURS :
// par defaut définit une droit // à x et passant par 0, et de dimension celle de l'espace de travail
Droite ();
// avec les datas
// vec est quelconque, non norme
Droite ( const Coordonnee& B, const Coordonnee& vec);
// avec la dimension: la droite générée par défaut est une droite qui passe par 0 et est // à l'axe des x
Droite (int dim);
// de copie
Droite ( const Droite& a);
// DESTRUCTEUR :
~Droite ();
// surcharge des operator
Droite& operator = ( const Droite & P);
// METHODES PUBLIQUES :
// retour le point de ref la droite
inline const Coordonnee& PointDroite() const { return A;};
// retourne le vecteur directeur de la droite (de norme = 1 )
inline const Coordonnee& VecDroite() const { return U;};
// change la dimension de la droite
void Change_dim(int dima) {A.Change_dim(dima);U.Change_dim(dima);};
// change le point de ref de la droite
void Change_ptref( const Coordonnee& B);
// change le vecteur de la droite
void Change_vect( const Coordonnee& vec);
// change toutes les donnees
void change_donnees( const Coordonnee& B, const Coordonnee& vec);
// calcul l'intercection 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 Intersection( const Droite & D,Coordonnee& M);
// calcul si un point donné appartient ou non à la droite (à la précision donnée)
bool Appartient_a_la_droite(const Coordonnee& B);
// ramene une normale, en 2D il y a une solution, en 3D a chaque appel on ramene
// une nouvelle normale calculée aleatoirement
Coordonnee UneNormale();
// calcul la distance d'un point à la droite
double Distance_a_la_droite(const Coordonnee& M) const;
// calcul du projeté d'un point sur la doite
Coordonnee Projete(const Coordonnee& M) const
{return (A+((M-A)*U) * U);};
// 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 DuMemeCote(const Coordonnee& M1, const Coordonnee& M2) const;
protected :
// VARIABLES PROTEGEES :
Coordonnee A; // un point de la droite
Coordonnee U; // vecteur directeur de la droite
static int alea; // une variables aleatoire qui sert pour la sortie d'une normale
// aleatoire en 3D
// METHODES PROTEGEES :
};
/// @} // end of group
#endif

1759
contact/ElContact.cc Normal file
View file

File diff suppressed because it is too large Load diff

418
contact/ElContact.h Normal file
View file

@ -0,0 +1,418 @@
// 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/>.
/************************************************************************
* DATE: 23/01/97 *
* $ *
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
* $ *
* PROJET: Herezh++ *
* $ *
************************************************************************
* BUT: Element de contact. *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
* VERIFICATION: *
* *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* ! ! ! ! *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* MODIFICATIONS: *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* $ *
************************************************************************/
#ifndef ELCONTACT_H
#define ELCONTACT_H
#include "Front.h"
#include "Noeud.h"
#include "Condilineaire.h"
#include "LesFonctions_nD.h"
/** @defgroup Groupe_sur_les_contacts
*
* BUT: groupe relatif aux contacts
*
*
* \author Gérard Rio
* \version 1.0
* \date 23/01/97
* \brief groupe relatif aux contacts
*
*/
/// @addtogroup Groupe_sur_les_contacts
/// @{
///
class ElContact
{
public :
// une classe structure de transfert pour simplifier le passage de paramètres
class Fct_nD_contact
{ public:
Fct_nD_contact();
Fct_nD_contact(const Fct_nD_contact& a);
~Fct_nD_contact();
Fct_nD_contact& operator= (const Fct_nD_contact& a);
// utilisation de fct nD
Fonction_nD * fct_nD_penalisationPenetration; // fct nD dans le cas d'une valeur pilotée
Fonction_nD * fct_nD_penetration_contact_maxi; // fct nD dans le cas d'une valeur pilotée
Fonction_nD * fct_nD_penetration_borne_regularisation; // fct nD dans le cas d'une valeur pilotée
Fonction_nD * fct_nD_force_contact_noeud_maxi; // fct nD dans le cas d'une valeur pilotée
Fonction_nD * fct_nD_penalisationTangentielle;
Fonction_nD * fct_nD_tangentielle_contact_maxi; // fct nD dans le cas d'une valeur pilotée
Fonction_nD * fct_nD_tangentielle_borne_regularisation; // fct nD dans le cas d'une valeur pilotée
Fonction_nD * fct_nD_force_tangentielle_noeud_maxi; // fct nD dans le cas d'une valeur pilotée
};
// CONSTRUCTEURS :
// par defaut
ElContact();
// la version avec fonction de pilotage nD
ElContact(Fct_nD_contact & fct_contact);
// fonction d'un pointeur d'element frontiere et d'un pointeur de noeud
// du fait éventuel qu'il peut-être collant ou pas
ElContact ( const Front * elfront, const Noeud * noeud, Fct_nD_contact & fct_contact_, int collant = 0);
// de copie
ElContact ( const ElContact & a);
// DESTRUCTEUR :
~ElContact ();
// METHODES PUBLIQUES :
// init d'une fonction nD pour le pilotage du type de contact 4
static void Init_fct_pilotage_contact4(Fonction_nD * fct_pilo_contact4)
{fct_pilotage_contact4 = fct_pilo_contact4;};
// méthode statique de modification éventuelle du type de contact utilisé localement
// par exemple pour le type 4 en fonction des itérations
static int Recup_et_mise_a_jour_type_contact() ;
// affichage à l'écran des informations liées au contact
void Affiche() const ;
// calcul d'un contact eventuel entre le noeud esclave et la frontiere
// ramene true s'il y a contact
// si init = false, on recherche le contact a partir du precedent point sauvegarde
// sinon on commence a l'aide d'element de reference,
// et on calcule et sauvegarde les coordonnées
// initiale locales theta^i du point de contact
// si le contact existe et si l'algo le demande (cf. ParaAlgoControle) :
// le noeud pourrait-être ramené sur la surface mais:
// on ne fait pas de projection, sinon on ne peut pas tester plusieurs contacts pour
// choisir le meilleur, puisque les choses changent entre avant et après le test de contact
// donc ici la position du noeud esclave n'est pas modifiée
bool Contact( bool init = true);
// juste après l'utilisation de la méthode Contact(), ramène le point en contact
const Coordonnee& Point_intersection() const { return Mtdt; };
// un stockage utilisé par les méthodes appelantes
int& Num_zone_contact() {return num_zone_contact;};
// calcul de la trajectoire a prendre en compte pour le contact
// ---- a) cas ou à t il n'y avait pas de contact, ou que l'on n'a pas fait de projection sur la surface (cas du contact cinématique)
// 4 cas : 1) le noeud bouge, dans ce cas la trajectoire est determinee
// par la variation de la position du noeud
// 2) le noeud est immobile, mais la frontiere bouge, la trajectoire est determine
// par une parallele a la variation moyenne de la frontiere (variation de G)
// 4) est une variante du 2), cas où l'on a une rotation autour de G, dans ce cas on prend comme trajectoire
// le maxi des déplacements de noeud
// 3) rien ne bouge, on utilise la normale au point de reference de l'element de frontiere
// pour calculer la trajectoire .
// dans le cas 3 la variable test = 0 , sinon elle vaut 1 pour le cas 1 et 2 pour le cas 2 , 4 pour le cas 4
// ---- b) cas ou à t on était déjà en contact avec projection sur la surface
// la trajectoire est alors systématiquement la direction de la dernière normale,
// retour : test=5
Coordonnee Trajectoire(int & test);
// calcul de l'Intersection de la trajectoire du noeud definit par le vecteur V
// avec l'element frontiere
// ramene les coordonnees du noeud projete
// dans le cas où il n'y a pas d'intersection, on ramène un point de dimension nulle
Coordonnee Intersection( const Coordonnee& V,bool init);
// construction de la condition lineaire de contact
// nb_assemb : indique le numéro d'assemblage correspondant
Condilineaire ConditionLi(int nb_assemb);
// ramene le tableau de tous les noeuds, le premier est celui esclave
Tableau <Noeud*>& TabNoeud() {return tabNoeud;};
const Tableau <Noeud*>& Const_TabNoeud() const {return tabNoeud;};
// ramene le tableau pour assemblage: dépend du Cas_solide()
// est cohérent avec TableauDdlCont
Tableau <Noeud*>& TabNoeud_pour_assemblage() {return tabNoeud_pour_assemblage;};
// retourne les tableaux de ddl associés aux noeuds, gere par l'element
// qui sont actifs au moment de la demande
// Tableau de DdlElement pour l'assemblage uniquement
const DdlElement & TableauDdlCont() const {return *ddlElement_assemblage;};
// ramene le noeud esclave
inline Noeud* Esclave() { return noeud;};
// ramene la force de contact pour consultation
const Coordonnee& Force_contact() const {return force_contact;};
// ramène les maxis concernant le noeud esclave
double F_N_MAX() const {return F_N_max;};
double F_T_MAX() const {return F_T_max;};
// ramène la pénalisation actuelle éventuelle
const double& Penalisation() const {return penalisation;};
const double& Penalisation_tangentielle() const {return penalisation_tangentielle;};
// ramène le déplacement tangentiel actuel
const Coordonnee& Dep_tangentiel() const {return dep_tangentiel;};
// ramène la normale actuelle
const Coordonnee& Normale_actuelle() const {return N;};
// idem pour les réactions sur les noeuds de la facette
const Tableau <Coordonnee>& TabForce_cont() const {return tabForce_cont;};
// actualisation de la projection du noeud esclave en fonction de la position de l'element
// maitre frontiere. Lorsque le noeud change d'element fontiere, on change l'element
// frontiere de l'element de contact en consequence
// dans le cas ou on ne trouve pas d'intersection, cas d'un noeud qui sort d'une zone de
// contact, on retourne false, sinon retourne true
// en fonction de la méthode de contact, le noeud est ramené éventuellement sur la frontière
bool Actualisation();
// ramene un pointeur sur l'element frontiere
inline Front * Elfront() const { return elfront;};
// permet de changer le deplacement maxi de tous les noeuds, qui sert
// pour éviter de considérer des contact trop éloignés
static void Change_dep_max(const double & deplac_max)
{dep_max = deplac_max * ParaGlob::param->ParaAlgoControleActifs().FacPourRayonAccostage();};
// test et met à jour le compteur de décollage du noeud
// si le noeud decolle ou non en fonction de la force de reaction
// ramene 1: s'il decolle
// 0: s'il ne décolle pas
bool Decol();
// change force: permet de changer la valeur de la force
// utile quand la force est calculée en dehors de l'élément de contact
void Change_force(const Coordonnee& force);
// idem pour les forces réparties sur la facette
void Change_TabForce_cont(const Tableau <Coordonnee>& tab) {tabForce_cont=tab;};
// gestion de l'activité
int Actif() const {return actif;}; // ramène l'activité du contact
void Met_Inactif() { actif = 0;nb_decol_tdt=0;}; // met en inactif
void Met_actif() { actif++;nb_decol_tdt=0;}; // met en actif une fois de plus
// ramène le nombre actuel de décolement
int Nb_decol() const {return nb_decol_tdt;};
// ramène le nombre de pénétration actuel
int Nb_pene() const {return nb_pene_tdt;};
// --- calcul des puissances virtuelles développées par les efforts de contact ----------
// et eventuellement calcul de la raideur associé
// -> explicite à tdt
virtual Vecteur* SM_charge_contact();
// -> implicite,
virtual Element::ResRaid SM_K_charge_contact();
// récupération des énergies intégrées sur l'éléments, résultants d'un précédent calcul
// explicite, ou implicite:
// 1- il s'agit ici de l'énergie développée par le frottement glissant ou pas
const EnergieMeca& EnergieFrottement() const {return energie_frottement;};
// 2- énergie de pénalisation (élastique a priori)
const double& EnergiePenalisation() const {return energie_penalisation;};
// cas d'une méthode avec pénalisation: calcul éventuel d'un pas de temps idéal,
// permettant de limiter la pénétration
// si oui retour de la valeur delta_t proposé
// sinon dans tous les autres cas retour de 0.
// le calcul se fait en fonction du pas de temps courant et de la pénétration
// donc nécessite que le contact ait déjà été étudié
double Pas_de_temps_ideal()const;
// ramène l'info sur le fait que le contact est avec un solide ou pas
// retour = 0 : contact bi déformable,
// = 1 le noeud est libre et la frontière est bloqué (solide)
// = 2 le noeud est bloqué (solide) la frontière est libre
// = 3: tout est bloqué (solide)
int Cas_solide() const {return cas_solide;};
// permet de modifier le contact entre collant ou non suivant "change"
void Change_contact_collant(bool change);
// récup de l'information concernant le contact collant ou pas
int Collant() const {return cas_collant;};
// récup des gaps calculés
const double & Gaptdt() const {return gap_tdt;};
const double & Dep_T_tdt() const {return dep_T_tdt;}
// mise à jour du niveau de commentaire
static void Mise_a_jour_niveau_commentaire(int niveau)
{niveau_commentaire = niveau;};
// récupération des ddl ou des grandeurs actives de tdt vers t
void TdtversT();
// actualisation des ddl et des grandeurs actives de t vers tdt
void TversTdt();
//----- lecture écriture de restart -----
void Lec_base_info_ElContact(ifstream& ent);
void Ecri_base_info_ElContact(ofstream& sort);
protected :
// VARIABLES PROTEGEES :
int actif; // un indicateur, disant si le contact est actif ou pas
// conteneur plutôt utilisé par les classes appelantes
// =1 -> premier contact, > 1 contact qui suit un contact
Front* elfront; // un element frontiere
Noeud * noeud; // un pointeur de noeud
int num_zone_contact; // un stockage uniquement utilisé par les méthodes appelantes
// pour éviter de le construire à chaque demande on définit un tableau de tous les noeuds
Tableau <Noeud*> tabNoeud; // tableau de tous les noeud, le premier est noeud
// le tableau des positions successives du noeud esclave, ceci pour effectuer une moyenne glissante
Tableau <Coordonnee> tab_posi_esclave;
// le nombre courant de positions de noeuds esclaves actuellement stockées
int nb_posi_esclave_stocker_t,nb_posi_esclave_stocker_tdt;
// indice dans tab_posi_esclave, du prochain stockage
int indice_stockage_glissant_t,indice_stockage_glissant_tdt;
Coordonnee Mtdt,Mt; // sauvegarde du point d'intercection s'il est recevable
Coordonnee M_noeud_tdt_avant_projection; // dans le cas d'une projection du noeud sur la surface maître
// sauvegarde des coordonnées du noeuds esclave: utilisation avec le type de contact 4
// les fonctions d'interpolation au premier point en contact: sert pour le contact collant
Vecteur phi_theta_0 ;
EnergieMeca energie_frottement; // énergie développée par le frottement glissant ou pas
double energie_penalisation; // énergie due à la pénalisation (élastique a priori)
// cas_solide permet de simplifier le contact dans le cas ou le maître ou l'esclave est solide
// sert pour diminuer la taille de la raideur uniquement
int cas_solide; // =0 contact bi déformable, =1 le noeud est libre et la frontière est bloqué (solide)
// = 2 le noeud est bloqué (solide) la frontière est libre
// = 3 tout est solide
int cas_collant; // prise en compte éventuelle d'un contact collant, si 0: contact normal
// = 1 : contact collant
Vecteur * residu; // residu local
Mat_pleine * raideur; // raideur locale
DdlElement * ddlElement_assemblage; // le ddlElement qui correspond
Tableau <Noeud*> tabNoeud_pour_assemblage; // tableau des noeuds qui servent pour l'assemblage
Coordonnee force_contact,force_contact_t; // force de contact sur le noeud principal
double F_N_max,F_T_max,F_N_max_t,F_T_max_t; // les maxis constatés
Tableau <Coordonnee> tabForce_cont,tabForce_cont_t; // le tableau des forces sur les noeuds de la facette
Coordonnee N; // la dernière normale calculée
Coordonnee dep_tangentiel; // le dernier déplacement tangentiel calculé
int nb_decol_t,nb_decol_tdt; // nombre de fois où le noeud décolle de manière consécutive
double gap_t,gap_tdt; // les pénétrations d'un pas à l'autre
double dep_T_t, dep_T_tdt; // les valeurs absolue des déplacements tangentiels d'un pas à l'autre
double nb_pene_t,nb_pene_tdt; // le nombre de penetration positives
// cas TypeCalculPenalisationPenetration() = 4 ou -4
// -> cas d'un facteur multiplicatif évolutif, on fait une moyenne glissante de 2
// on mémorise donc les facteurs multiplicatifs successifs
double mult_pene_t,mult_pene_tdt;
// cas TypeCalculPenalisationTangentielle() = 4 ou -4
// -> cas d'un facteur multiplicatif évolutif, on fait une moyenne glissante de 2
// on mémorise donc les facteurs multiplicatifs successifs
double mult_tang_t,mult_tang_tdt;
double penalisation,penalisation_tangentielle; // on sauvegarde la pénalisation
// l'intérêt est de pouvoir la visualiser
double penalisation_t,penalisation_tangentielle_t; // les grandeurs à t
// utilisation de fct nD
Fct_nD_contact fct_nD_contact;
// pour le contact 4, pour le calcul de la pénalisation avec une moyenne glissante
Tableau <double > val_penal; // tableau de stockage intermédiaire
int pt_dans_val_penal; // indice du prochain elem du tableau a remplir
// stocke le dernier type de trajectoire du noeud par rapport à la facette
// c-a-d la valeur de la variable test dans la fonction Trajectoire(int & test)
int type_trajectoire_t,type_trajectoire_tdt;
// concernant les enum de ddl associées aux éléments de contact, on définit un tableau global
// qui est utilisé par tous les éléments
static list <DdlElement> list_Ddl_global; // liste de tous les DdlElements des éléments de contact
static list <Vecteur> list_SM; // list de seconds membres locals: sert pour tous les éléments de contact
static list <Mat_pleine> list_raideur; // list de raideurs locales: " " " "
// stockage du maximum de distance tolérée entre noeud à tdt et le projeté, sert pour éliminer les contacts aberrants
static double dep_max;
static int niveau_commentaire;
// stockage d'une fonction nD pour le pilotage du type de contact 4
static Fonction_nD * fct_pilotage_contact4;
// METHODES PROTEGEES :
// calcul la normal en fonction de differente conditions
Coordonnee Calcul_Normale (int dim, Plan & pl, Droite & dr,int indic);
void Libere();
// construction du tableau de tous les noeuds, le premier est celui esclave
// et mise à jour de ddlElement et de list_Ddl_global éventuellement
void Construction_TabNoeud();
// récupération d'informations des classes internes pour le calcul du résidu
// N: le vecteur normal
// M_impact: le point d'impact sur la surface (ou ligne)
// phii : les fonctions d'interpolation au point d'impact
// si avec_var est vrai: il y a retour du tableau de variation de la normale
Tableau <Coordonnee >* RecupInfo(Coordonnee& N,Coordonnee& M_impact,Vecteur& phii,bool avec_var );
// mise à jour de cas_solide et donc de ddlElement en fonction de l'activité des ddl
// mise à jour du tableau de noeud pour l'assemblage tabNoeud_pour_assemblage
void Mise_a_jour_ddlelement_cas_solide_assemblage();
// récup d'une place pour le résidu local et mise à jour de list_SM éventuellement
void RecupPlaceResidu(int nbddl);
// récup d'une place pour la raideur locale et mise à jour de list_raideur éventuellement
void RecupPlaceRaideur(int nbddl);
// calcul du facteur de pénalisation en pénétration, en fonction de la géométrie
// du module de compressibilité et des différents possibles
// éventuellement, calcul de la dérivée: d_beta_gapdu, du facteur par rapport au gap
// la sensibilité dépend du type de calcul du facteur de pénalisation
double CalFactPenal(const double& gap,double & d_beta_gap,int contact_type);
// calcul du facteur de pénalisation tangentielle, en fonction de la géométrie
// du module de compressibilité et des différents cas possibles
// éventuellement, calcul de la dérivée: d_beta_gapdu, du facteur par rapport au gap
// la sensibilité dépend du type de calcul du facteur de pénalisation
double CalFactPenalTangentiel(const double& gap,double & d_beta_gap);
// limitation éventuelle au niveau de la pénétration maxi
//void Limitation_penetration_maxi(
// calcul éventuel de la moyenne glissante des positions successive du noeud esclave
void ChangeEnMoyGlissante(Coordonnee& Noe_atdt);
// calcul de la moyenne glissante de la pénalisation
void Moyenne_glissante_penalisation(double& penalisation, double& essai_penalisation);
// calcul d'une fonction nD relative à des données de contact
double Valeur_fct_nD(Fonction_nD * fct_nD) const;
};
/// @} // end of group
#endif

2655
contact/ElContact_2.cc Executable file
View file

File diff suppressed because it is too large Load diff

2193
contact/LesContacts.cc Normal file
View file

File diff suppressed because it is too large Load diff

456
contact/LesContacts.h Normal file
View file

@ -0,0 +1,456 @@
// 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/>.
/************************************************************************
* DATE: 23/01/97 *
* $ *
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
* $ *
* PROJET: Herezh++ *
* $ *
************************************************************************
* BUT: Creation et gestion des contacts. *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
* VERIFICATION: *
* *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* ! ! ! ! *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* MODIFICATIONS: *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* $ *
************************************************************************/
#ifndef LESCONTACTS_H
#define LESCONTACTS_H
#include "Front.h"
#include <list>
#include "List_io.h"
#include "ElContact.h"
#include "Condilineaire.h"
#include "Nb_assemb.h"
#include "LesMaillages.h"
#include "LesReferences.h"
#include "Basiques.h"
#include "TypeQuelconque.h"
#include "Temps_CPU_HZpp.h"
/// @addtogroup Groupe_sur_les_contacts
/// @{
///
class LesContacts
{
public :
// CONSTRUCTEURS :
// constructeur par defaut
LesContacts ();
// constructeur de copie
LesContacts (const LesContacts& a);
// DESTRUCTEUR :
~LesContacts ();
// METHODES PUBLIQUES :
// ramène la liste des éléments de contact
LaLIST <ElContact>& LesElementsDeContact() {return listContact;};
// initialisation des zones de contacts éventuelles à partir des éléments de frontières et des noeuds esclaves
// sauf les frontières qui sont les mêmes pendant tout le calcul
// --- en entrée:
// les maillages, les ref et les fonctions nD
// -- en sortie:
// cas du contact type 4 : on renseigne éventuellement une fonction de pilotage --> fct_pilotage_contact4
// récup de: nb_mail_Esclave, nbmailMaitre,
// récup du tableau "indice" : = la liste pour chaque noeud, des éléments qui contient ce noeud
// création des conteneurs internes : tesctotal, tesN_collant, t_listFront, tesN_encontact
void Init_contact(LesMaillages& lesMail
,const LesReferences& lesRef
,LesFonctions_nD* lesFonctionsnD);
// mise à jour du tableau "indice": là on utilise les numéros de noeuds qui peuvent changer
// via une renumérotation. Le tableau indice, est un tableau utilisé pour les recherches
// mais le numéro de noeud n'est pas utilisé pour les stockages divers (on utilise un pointeur de noeud)
// du coup le tableau indice peut évoluer ex: après un remaillage avec chg de num de noeud
// il doit donc être remis à jour avant l'étude de nouveau contact
// la liste pour chaque noeud, des éléments qui contient ce noeud
// indice(i)(j) : = la liste des éléments qui contiennent le noeud j, pour le maillage i
void Mise_a_jour_indice(const Tableau < const Tableau <List_io < Element* > > *> ind)
{ indice = ind;};
// indique si l'initialisation a déjà été effectuée ou pas, car celle-ci ne doit-être faite
// qu'une fois
bool Init_contact_pas_effectue() const { return (tesctotal.Taille() == 0);};
// verification qu'il n'y a pas de contact avant le premier increment de charge
void Verification();
// definition des elements de contact eventuels
// - imposition des conditions de non penetration
// dep_max : déplacement maxi des noeuds du maillage
// , sert pour def des boites d'encombrement maxi des frontières
// ramène true s'il y a effectivement création d'élément de contact
bool DefElemCont(double dep_max);
// reexamen du contact pour voir s'il n'y a pas de nouveau element de
// contact
// ramene false s'il n'y a pas de nouveau element de contact
// true sinon
// dep_max : déplacement maxi des noeuds du maillage
// , sert pour def des boites d'encombrement maxi des frontières
bool Nouveau(double dep_max);
// suppression définitive, si le contact à disparu, des éléments inactifs
// ramène false si aucun élément n'est finalement supprimé
bool SuppressionDefinitiveElemInactif();
// relachement des noeuds collés
// ramène true s'il y a des noeuds qui ont été relachés
bool RelachementNoeudcolle();
// definition de conditions lineaires de contact
// marquage des ddl correspondant aux directions bloquees s'il s'agit d'un contact de type 1
// casAssemb : donne le cas d'assemblage en cours
//
const Tableau <Condilineaire>& ConditionLin(const Nb_assemb& casAssemb);
// création d'un tableau de condition linéaire, correspondant à tous les éléments de contact en cours
// qu'ils soient actifs ou pas (a prior cette méthode est conçu pour donner des infos relativement à la largeur
// de bandes en noeuds due aux CLL)
// chacune des condition ne contient "d'exploitable" que le tableau de noeuds associés à la CLL,
const Tableau <Condilineaire>& ConnectionCLL();
// effacement du marquage de ddl bloque du au conditions lineaire de contact
void EffMarque();
// indique si les surfaces des maillages maîtres ont des déplacements fixés
// c-a-d sont de type solide imposé
// retourne true si les déplacements des maîtres sont imposés
bool Maitres_avec_deplacement_imposer() const;
// def de la largeur de bande des elements contacts
// casAssemb : donne le cas d'assemblage a prendre en compte
// les condi linéaires ne donnent pas des largeurs de bande sup et inf égales !!!
// demi = la demi largeur de bande maxi ,
// total = le maxi = la largeur sup + la largeur inf +1
// cumule = la somme des maxis, ce qui donnera la largeur finale, due à des multiples multiplications: une par conditions linéaires
// ceci dans le cas de la prise en compte par rotation (et uniquement dans ce cas)
// en retour, ramène un booleen qui :
// = true : si la largeur de bande en noeud est supérieure à 1 (cas d'un contact avec une surface déformable)
// = false : si non, ce qui signifie dans ce cas qu'il n'y a pas d'augmentation de la largeur
// en noeud (cas d'un contact avec une surface rigide)
bool Largeur_Bande(int& demi,int& total,const Nb_assemb& casAssemb,int& cumule);
// actualisation du contact, on examine que les elements de contact, dont on
// actualise la projection du noeud esclave en fonction de la position de l'element
// maitre frontiere (mais la position finale du noeud n'est pas forcément changée, cela dépend du
// modèle de contact (cinématique, pénalisation etc.)
// dans le cas où la réaction est négative, en fonction de l'algorithme l'élément de contact est
// inactivé, cependant tous les éléments de contact sont passés en revue (actif ou pas)
// ramène true si quelque chose à changé, false sinon
bool Actualisation();
// ramène une liste de noeuds dont la position a été perturbé par le contact
// (dépend du type de contact : ex cas contact = 4)
// la liste passée en paramètre est supprimée et remplacée par la liste résultat
void Liste_noeuds_position_changer(list <Noeud * >& li_noe);
// calcul des reactions de contact et stockage des valeurs
// solution : le vecteur residu
// test d'un decollement eventuelle, pour un noeud en contact
// ramene true s'il y a decollement, sinon false
// casAssemb : donne le cas d'assemblage en cours
void CalculReaction(Vecteur& residu,bool& decol,const Nb_assemb& casAssemb
,bool affiche);
// récupération via les éléments de contact des forces maxis
// un : le maxi en effort normal, deux: le maxi en effort tangentiel
DeuxDoubles Forces_contact_maxi(bool affiche) const;
// récupération via les éléments de contact des gaps maxis
// un : le maxi en gap normal, deux: le maxi en gap tangentiel
DeuxDoubles Gap_contact_maxi(bool affiche) const;
// cas d'une méthode avec pénalisation: calcul éventuel d'un pas de temps idéal,
// si oui retour de la valeur delta_t proposé
// sinon dans tous les autres cas retour de 0.
// le calcul se fait en fonction du pas de temps courant, des forces de réaction et de la pénétration
// donc nécessite que le contact ait déjà été étudié et que les efforts de contact ait été calculé
double Pas_de_temps_ideal()const;
// calcul d'une estimation du pas de temps critique du aux éléments de contact
// affichage des reactions de contact sur la sortie
void Affiche(ofstream& sort) const ;
// affichage à l'écran des informations liées au contact
void Affiche() const ;
// affichage et definition interactive des commandes
void Info_commande_LesContacts(UtilLecture & entreePrinc);
// lecture éventuelle des zones où le contact doit être recherché, à l'exclusion de tout
// autre zone, ainsi que la définition de l'auto-contact éventuel
// ainsi que des contacts solide-deformables éventuel
void Lecture_zone_contact(UtilLecture & entreePrinc,const LesReferences& lesRef);
// récupération des ddl ou des grandeurs actives de tdt vers t
void TdtversT();
// actualisation des ddl et des grandeurs actives de t vers tdt
void TversTdt();
//------- temps cpu -----
// retourne temps cumulé pour imposer les CL imposées
const Temps_CPU_HZpp& Temps_cpu_Contact() const {return tempsContact;};
//----- lecture écriture de restart -----
// cas donne le niveau de sauvegarde
void Ecri_base_info_LesContacts(ofstream& sort);
// on utilise deux pointeurs de fonctions qui permettent de récupérer le pointeur de noeud esclave
// idem au niveau de l'élément
template <class T>
void Lec_base_info_LesContacts(ifstream& ent
,T& instance // l'instance qui permet d'appeler les pointeurs de fonctions
,Noeud& (T::*RecupNoeud)(int i, int j) const
,Element& (T::*RecupElement_LesMaille) (int i, int j) const);
// méthode générale: récupération des grandeurs particulière (hors ddl )
// correspondant à liTQ
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
// ===> n'est pas utilisée dans le cas du contact, c'est la méthode spécifique
// Mise_a_jour_Pour_Grandeur_particuliere qui la remplace (qui n'est pas en const car elle
// modifie les conteneurs des noeuds et éventuellement éléments)
void Grandeur_particuliere
(bool absolue,List_io<TypeQuelconque>& ,Loi_comp_abstraite::SaveResul * ,list<int>& decal) const
{};
// Il s'agit ici de mettre à jour les conteneurs stockés aux noeuds et/ou aux éléments
// qui servent à récupérer les infos liés aux contact correspondant à liTQ
// actuellement les conteneurs passés en paramètre ne servent que pour
// les énumérés, et les informations résultantes sont stockées au niveau des noeuds
// constituant les éléments de contact
//--> important : les conteneurs sont supposés initialisés avec l'appel
void Mise_a_jour_Pour_Grandeur_particuliere(
List_io < TypeQuelconque >& li_restreinte_TQ
);
// récupération de la liste de tous les grandeurs particulières disponibles avec le contact
// cette liste est identique quelque soit le maillage: il n'y a donc pas de tableau indicé du num de maillage
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
List_io<TypeQuelconque> ListeGrandeurs_particulieres(bool absolue) const;
// concernant les grandeurs gérées par le contact:
// ramène une liste d'iterator correspondant aux List_io<TypeQuelconque> passé en paramètre
// idem pour une List_io < Ddl _enum_etendu >
void List_reduite_aux_contact(const List_io<TypeQuelconque>& liTQ
,List_io < TypeQuelconque >& li_restreinte_TQ
);
// initialisation des listes de grandeurs qu'ils faudra transférérer aux niveaux des noeuds des élements
// de contact, on définit si besoin, les conteneurs ad hoc au niveau des noeuds
// ok, mais à revoir sans doute cf. pense bete 14 oct
void Init_Grandeur_particuliere (bool absolue,List_io<TypeQuelconque>& li1);
// ne sert plus (à virer car problématique !! )
// mise à jour du stockage inter, pour prendre en
// compte une nouvelle numérotation des noeuds
// void Prise_en_compte_nouvelle_numerotation_noeud();
// initialisation de la liste de grandeurs qui sont effectivement gérées par le contact
// ok, mais à revoir sans doute cf. pense bete 14 oct
// List_io<TypeQuelconque> Init_list_grandeur_contact_a_sortir(const Tableau <List_io < TypeQuelconque > >& li1);
class ReactCont
{ // surcharge de l'operator de lecture
friend istream & operator >> (istream &, ReactCont &);
// surcharge de l'operator d'ecriture
friend ostream & operator << (ostream &, const ReactCont &);
public :
Noeud* noe; // noeud esclave
Coordonnee force ; // force de reaction au noeud esclave
Tableau <Noeud *> tabNoeud; // les noeuds de la frontiere maitre
Tableau <Coordonnee> tabForce; // les reac "" ""
// constructeur par defaut
ReactCont() ;
// constructeur en fonction des datas du noeud esclave seul
ReactCont(Noeud* no,const Coordonnee& forc) ;
// constructeur en fonction des datas de tous les noeuds
ReactCont(Noeud* no,const Coordonnee& forc,Tableau <Noeud *> tN,const Tableau <Coordonnee>& tFor);
// constructeur de copie
ReactCont(const ReactCont& a);
// affectation
ReactCont& operator = (const ReactCont& a);
// test
bool operator == (const ReactCont& a);
bool operator != (const ReactCont& a);
};
//----------------------------------------------------------------------------------------------------
private :
// VARIABLES PROTEGEES de la classe LesContacts:
LesFonctions_nD* sauve_lesFonctionsnD ; // sauvegarde à l'initialisation (méthode Init_contact)
ElContact::Fct_nD_contact fct_nD_contact; // fonctions nD de pilotage: peuvent ne pas exister
LaLIST <ElContact> listContact; // la liste des elements en contact
LaLIST <LaLIST <ElContact>::iterator> listContact_nouveau_tatdt; // la liste des nouveaux contacts qui sont apparus sur l'incrément
LaLIST <ElContact> listContact_efface_tatdt; // la liste des contacts effacés sur l'incrément
int nb_contact_actif; // nombre de contact actif courant
// list <TroisEntiers> numtesN; // .un : maillage, .deux: num zone de contact, .trois: num propre du noeud
// // la liste des numéros dans tesN des noeuds en contact
Tableau <ReactCont> tabReacCont; // les reactions
Tableau <Condilineaire> tabCoLin; // tableau des conditions lineaires
static MotCle motCle; // liste des mots clés
// la liste des éléments qui contiennent des frontières, est reconstitué au démarrage avec Init_contact(..
list <Element *> liste_elemens_front;
// la liste pour chaque noeud, des éléments qui contient ce noeud : construite avec Init_contact
// indice(i)(j) : = la liste des éléments qui contiennent le noeud j, pour le maillage i
Tableau < const Tableau <List_io < Element* > > *> indice;
//--------- les tableaux de gestions pour la recherche de contact ----------------
list <quatre_string_un_entier> nom_ref_zone_contact; // liste des noms de références des zones de contact éventuelle
// cette liste peut être vide, dans ce cas cela signifie que toutes les frontières des pièces sont
// suceptible de rentrer en contact. Dans le cas où la liste n'est pas vide, seules les grandeurs
// référencées sont succeptibles de rentrer en contact
// nom1 -> nom_mail_ref_zone_contact pour les noeuds
// nom2 -> le nom de la référence de noeuds
// nom3 -> nom_mail_ref_zone_contact pour les frontières
// nom4 -> le nom de la référence de frontière
// n -> l'entier = 0 : pas d'utilisation
// = 1 : indique qu'il faut considérer un contact collant (glue)
// la liste des éléments frontières succeptibles d'entrer en contact
// ces Front (qui contiennent que des pointeurs sauf une boite d'encombrement)
// sont différents de ceux des maillages, et sont donc stocké en interne
Tableau < Tableau < LaLIST_io <Front> > > t_listFront;
// t_listFront(i)(j)(k) : maillage maître (i)
// zone de contact (j)
// (K) = kième frontière (dans l'ordre de la liste)
// ** pour le tableaux t_listFront, le numéros dit de maillage, n'est pas le numéro
// ** intrinsèque de maillage (telle que ceux associés aux noeuds et éléments)
// ** mais uniquement un numéro locale d'ordre
// ** mais on a: les éléments de frontière de t_listFront(i) font partie du maillage
// i + (nb_mail_Esclave-nbmailautocontact)
// les noeuds esclaves potentiels
Tableau < Tableau < Tableau < Noeud*> > > tesctotal;
// tesctotal(i)(j)(k) : maillage esclave (i), zone de contact (j)
// noeud esclave (k) : k= le numéro d'ordre dans tesctotal(i)(j),
// ==>> ce n'est pas le numéro du noeud dans le maillage !
// indicateur de noeuds collants
Tableau < Tableau < Tableau < int> > > tesN_collant;
// tesN_collant(i)(j)(k): maillage esclave (i), zone de contact (j)
// noeud esclave (k):k= le numéro d'ordre dans tesctotal(i)(j),
// ==>> ce n'est pas le numéro du noeud dans le maillage !
// indique pour si le noeud esclave doit être
// considéré en contact collant (=1) ou pas (=0)
// tesN_encontact: globalise tous les contacts sur un noeud (indépendamment des zones de contact)
// tesN_encontact(i)(num_noe) : nombre de contact avec le noeud
// ancien stockage: Tableau < Tableau < LaLIST < LaLIST<ElContact>::iterator > > > tesN_encontact;
// ancien stockage // tesN_encontact(i)(num_noe): maillage (i), noeud -> num_noe:
// indique pour le noeud esclave s'il est en contact ou non via la taille de la liste associée
// utilisation: pour chaque Noeud* = tesctotal(i)(k) -> la liste des éléments en contact
// contenant le noeud k
// on remplace par une map: l'avantage c'est que l'on utilise plus les numéros de noeuds pour
// retrouver la liste, mais on utilise à la place la valeur pointeur de noeud esclave qui
// doit être unique, du coup c'est indépendant de la numérotation des noeuds -> permet de la renumérotation
// sans changer la map
Tableau < std::map<Noeud*,LaLIST < LaLIST<ElContact>::iterator > > > tesN_encontact;
// tesN_encontact(numMail_esclave)[*pt_noeud] -> la liste des iterators d'élément en contact
// avec le noeud
//--------- fin tableaux de gestions pour la recherche de contact ----------------
int nb_mail_Esclave; // indique le nombre de maillages esclaves
int nbmailautocontact; // indique le nombre de maillages esclaves en auto contact, donc qui joue le rôle esclave et maître, ou maillages mixte: dépend de la manière dont les zones de contact ont été définit
int nbmailMaitre; // indique le nombre de maillages maitres
// = nombre total de maillage - (nb_mail_Esclave-nbmailautocontact)
// si l'on considère l'ensemble des maillages du calcul on a successivement:
// les nb_mail_Esclave-nbmailautocontact premier maillages: purement esclave
// puis les nbmailautocontact qui sont en auto contact: ils jouent le rôle d'esclaves et maîtres en même temps
// puis nbmailMaitre-nbmailautocontact : purement maître
list < Deux_String > cont_solide_defor; // liste des noms de maillages définissants les contacts
// solide-déformable: (*ie).nom1 -> le maillage solide , (*ie).nom2 -> le maillage deformable
// ---- pour le post traitement ---
List_io<TypeQuelconque> liQ_en_sortie; // liste de grandeurs quelconque qu'il faut sortir
// -------- une variable de travail pour la méthode LesContacts::ConnectionCLL()---
Tableau <Condilineaire> t_connectionCLL;
//------- temps cpu -----
// retourne temps cumulé pour imposer les CL imposées
Temps_CPU_HZpp tempsContact;
// METHODES PROTEGEES :
// mise à jour des boites d'encombrement pour les éléments qui contiennent une frontière
void Mise_a_jour_boite_encombrement_element_contenant_front();
// suppression du gap de contact pour les noeuds "collant avec suppression de gap"
void Suppression_gap_pour_noeud_collant();
// récupération de la zone de contact d'un élément de contact existant
// c'est un peu lourdinge, mais l'avantage c'est que cela s'adapte à une situation qui
// a par exemple changé
// si en retour le numéro de zone = 0, cela signifie que le contact ne peut plus exister
int Recup_ref( ElContact& al) ;
// récupère le nombre de contact actif et met à jour ce nombre
// normalement devrait toujours être correct mais ?? il y a quelque chose d'incorrecte quelque part
int Recalcul_Nb_contact_actif();
// création du conteneur Fct_nD_contact
void Creation_Fct_nD_contact();
};
/// @} // end of group
// pour faire de l'inline: nécessaire avec les templates
// on n'inclut que les méthodes templates
#include "LesContacts_2.cc"
#define LesContacts_2_deja_inclus
#endif

176
contact/LesContacts_2.cc Executable file
View file

@ -0,0 +1,176 @@
// 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/>.
#ifndef LesContacts_2_deja_inclus
#include "LesContacts.h"
#endif
#ifndef LesContacts_2_deja_inclus
#include <vector>
#include "ReferenceNE.h"
#include "ReferenceAF.h"
//----- lecture écriture base info -----
// lecture base info
// on utilise deux pointeurs de fonctions qui permettent de récupérer le pointeur de noeud esclave
// idem au niveau de l'élément
template <class T>
void LesContacts::Lec_base_info_LesContacts(ifstream& ent
,T& instance // l'instance qui permet d'appeler les pointeurs de fonctions
,Noeud& (T::*RecupNoeud)(int i, int j) const
,Element& (T::*RecupElement_LesMaille)(int i, int j) const)
{
int niveau_commentaire_contact = ParaGlob::param->ParaAlgoControleActifs().Niveau_commentaire_contact();
// tout d'abord on lit le type
string type_contact,no_taille; int taille=0; // taille par défaut
ent >> type_contact >> no_taille >> taille;
#ifdef MISE_AU_POINT
if (type_contact != "LesContacts:")
{ cout << "\n erreur dans la lecture des contacts, on attendait la chaine: LesContacts:"
<< " alors que l'on a lue : " << type_contact
<< "\n LesContacts::Lec_base_info_LesContacts(i...";
Sortie(1);
};
if (no_taille != "taille")
{ cout << "\n erreur dans la lecture des contacts, on attendait la chaine: taille"
<< " alors que l'on a lue : " << no_taille
<< "\n LesContacts::Lec_base_info_LesContacts(i...";
Sortie(1);
};
#endif
// on supprime tous les éléments de contact sauf... les collants qui restent collant
{
LaLIST <ElContact>::iterator iE ;
// on met .end(), car cette borne peut changer au gré des suppression
for (iE = listContact.begin(); iE != listContact.end(); iE++)
if (!((*iE).Collant()) // et si ce n'est pas collant: si collant, on ne change rien
)
{ // sinon on supprime
LaLIST <ElContact>::iterator iiE = iE;
iE--; // pointe sur le precedent element
LaLIST < LaLIST<ElContact>::iterator > & list_tesN
= tesN_encontact((*iiE).Esclave()->Num_Mail())[(*iiE).Esclave()];
list_tesN.remove(iiE);
listContact.erase(iiE); // efface l'element
};
Recalcul_Nb_contact_actif();
};
// nb_contact_actif = 0 ; // init
Tableau <int> compteur_inesc(nb_mail_Esclave,1); // compteur de compteur de noeuds esclaves
// on initialise la liste des éléments en contact
// listContact.clear(); // comme on redémarre, on efface la situation actuelle
// non, il faut garder les éléments collants
// --- on s'occupe du deuxième niveau de tesN_encontact
// une map vide par défaut
std::map<Noeud*,LaLIST < LaLIST<ElContact>::iterator > > map_vide;
tesN_encontact.Change_taille(nb_mail_Esclave,map_vide);
for (int i=1;i<= taille;i++)
{// on commence par créer un élément de contact de base en fonction des infos sauvegardées
// --- le noeud esclave
string toto; int numMail_esclave, numNoeud_esclave;
string type;
int numMail_element, num_element, num_front_element;
Enum_type_geom enu ; string nom_enu; //N_es 1 72 El 2 2 SURFACE 1
ent >> toto >> numMail_esclave >> numNoeud_esclave // le noeud esclave
>> type >> numMail_element >> num_element // partie relative à l'élément
>> enu >> num_front_element; // et la frontière
// on ne continue que si le noeud fait partie des maillages esclave ou en auto-contact
// le nombre de maillages esclave et en auto-contact, peut-être redéfinit lors d'un restart
// donc n'est pas sauvegardé
if (numMail_esclave <= nb_mail_Esclave)
{ Noeud& no = (instance.*RecupNoeud)(numMail_esclave, numNoeud_esclave); // ""
// --- création de la frontière
Element& elem = (instance.*RecupElement_LesMaille)(numMail_element,num_element);
ElFrontiere* elfr;
switch (enu)
{ case POINT_G : elfr = elem.Frontiere_points(num_front_element,true) ; break;
case LIGNE : elfr = elem.Frontiere_lineique(num_front_element,true) ; break;
case SURFACE : elfr = elem.Frontiere_surfacique(num_front_element,true) ; break;
default :
cout << "\nErreur : valeur incorrecte du type de frontiere " << Nom_type_geom(enu) << " !\n";
cout << "LesContacts::Lec_base_info_LesContacts(... \n";
Sortie(1);
};
// création d'un Front
Front froont(*elfr,&elem,num_front_element);
// création de l'élément de contact
ElContact courant(&froont,&no,fct_nD_contact);
// lecture des données particulières
courant.Lec_base_info_ElContact(ent);
// on regarde si c'était un contact actif ou pas, si oui on incrémente le nombre
if (courant.Actif())
nb_contact_actif++;
// sauvegarde
listContact.push_front(courant);
LaLIST<ElContact>::iterator ipp = listContact.begin();
// on met à jour le tableau tesN_encontact
//*** la liste peut très bien ne pas exister !
// #ifdef MISE_AU_POINT
// if (tesN_encontact(numMail_esclave).find((*ipp).Esclave())
// == tesN_encontact(numMail_esclave).end() )
// { cout << "\n*** Erreur : on ne trouve pas la liste d'element en contact avec le noeud esclave "
// << (*ipp).Esclave()->Num_noeud() << " du maillage " << numMail_esclave
// << " la suite n'est pas possible "
// << " LesContacts::Lec_base_info_LesContacts(.. \n";
// Sortie(1);
// };
// #endif
// ajout dans la liste associé au noeud esclave
tesN_encontact(numMail_esclave)[(*ipp).Esclave()].push_front(listContact.begin());
} // fin de l'ajout d'un élément de contact
else
{if (ParaGlob::NiveauImpression() > 0)
cout << "\n *** attention, le maillage "<<numMail_esclave<<" n'est plus esclave !"
<< " l'element de contact non pris en compte: "
<< "\n N_es " <<numMail_esclave << " "<<numNoeud_esclave
<< type << " "<< numMail_element << " "<< num_element << " "
<< Nom_type_geom(enu) <<" "<<num_front_element;
};
};
if (((ParaGlob::NiveauImpression() > 2) && (nb_contact_actif != 0))
|| (niveau_commentaire_contact > 0 )
)
cout << "\n >> restart: lecture elements contact : "<< nb_contact_actif << " elements lus " << endl ;
string mot;
ent >> mot >> tempsContact;
};
#endif // fin du test d'inclusion : ifndef LesContacts_2_deja_inclus

1247
contact/LesContacts_3.cc Executable file
View file

File diff suppressed because it is too large Load diff

211
contact/Plan.cc Normal file
View file

@ -0,0 +1,211 @@
//#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 "Plan.h"
#include "ConstMath.h"
#include "MathUtil.h"
#include "Util.h"
// CONSTRUCTEURS :
// par defaut
Plan::Plan () :
A(), N()
{};
// avec les datas
Plan::Plan ( const Coordonnee& B, const Coordonnee& vec) :
A(B), N(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 << "\nPlan::Plan (Coordonnee& B,Coordonnee& vec)" << endl;
Sortie(1);
};
#endif
double d = N.Norme();
if (d <= ConstMath::petit)
{ cout << "\nErreur : la norme du vecteur est trop petite !";
cout <<"\nnorme = " << d;
cout << "\nPlan::Plan (Coordonnee& B,Coordonnee& vec)" << endl;
Sortie(1);
};
N /= d; // N est ainsi norme
};
// avec la dimension
Plan::Plan (int dim) :
A(dim), N(dim)
{ };
// de copie
Plan::Plan ( const Plan& a) :
A(a.A), N(a.N)
{};
// DESTRUCTEUR :
Plan::~Plan () {};
// surcharge des operator
Plan& Plan::operator = ( const Plan & P)
{ this->A = P.A; this->N = P.N;
return *this;
};
// METHODES PUBLIQUES :
// change le point de ref du plan
void Plan::Change_ptref( const Coordonnee& B)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != N.Dimension())
{ cout << "\nErreur : les dimensions du point et du vecteur ne sont pas identique !";
cout <<"\ndim point = " <<B.Dimension() <<", dim vecteur =" << N.Dimension();
cout << "\nPlan::Change_ptref(Coordonnee& B)" << endl;
Sortie(1);
};
#endif
A = B;
};
// change le vecteur normal du plan
void Plan::Change_normal( 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 << "\nPlan::Change_normal(Coordonnee& vec)" << endl;
Sortie(1);
};
#endif
double d = vec.Norme();
if (d <= ConstMath::petit)
{ cout << "\nErreur : la norme du vecteur est trop petite !";
cout <<"\nnorme = " << d;
cout << "\nPlan::Change_normal(Coordonnee& vec)" << endl;
Sortie(1);
};
N = vec / d; // N est ainsi norme
};
// change toutes les donnees
void Plan::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 << "\nPlan::change_donnees(Coordonnee& B,Coordonnee& vec)" << endl;
Sortie(1);
};
#endif
A = B;
double d = vec.Norme();
if (d <= ConstMath::petit)
{ cout << "\nErreur : la norme du vecteur est trop petite !";
cout <<"\nnorme = " << d;
cout << "\nPlan::change_donnees(Coordonnee& B,Coordonnee& vec)" << endl;
Sortie(1);
};
N = vec / d; // N est ainsi norme
};
// calcul l'intercection M d'une droite avec le plan, ramene 0 s'il n'y
// a pas d'intercection, ramene -1 si l'intercection ne peut pas etre calculee
// et 1 s'il y a un point d'intercection
int Plan::Intersection( const Droite & D,Coordonnee& M)const
{ const Coordonnee & B = D.PointDroite(); // par commodite d'ecriture
const Coordonnee & V = D.VecDroite(); // ""
Coordonnee AB = B - A;double Nab = AB.Norme();
// existance du point d'intercection et cas particulier
// - droite // au plan ?
double abN = AB * N;
double vn = V * N;
//cout << "\n **Plan::Intersection debug ** ";
//cout << "\n B ";B.Affiche();
//cout << "\n V ";V.Affiche();
//cout << "\n N ";N.Affiche();
//cout << "\n M ";M.Affiche();
//cout << "\n abN= "<<abN << " vn= "<<vn << endl;
//
if ( Dabs(vn) <= ConstMath::pasmalpetit)
// cas //
// - on regarde si la droite appartiend au plan
if ( Dabs(abN) <= ConstMath::pasmalpetit)
// droite dans le plan -> indetermination
return -1;
else
// il n'y a pas d'intercection
return 0;
else
// cas ou il y a intersection
{ M = B - (abN /vn) * V;
return 1;
}
};
// calcul la distance d'un point à la droite
double Plan::Distance_au_plan(const Coordonnee& M) const
{ int dima = ParaGlob::Dimension();
Coordonnee AM(M-A);
return Dabs(AM * N);
};
// ramène true si les deux points sont du même coté du plan, false sinon
bool Plan::DuMemeCote(const Coordonnee& M1, const Coordonnee& M2) const
{ Coordonnee AM1 = M1-A;
Coordonnee AM2 = M2-A;
return ((AM1 * N) * (AM2 * N) >= 0);
};
// surcharge de l'operateur de lecture
istream & operator >> (istream & entree, Plan & pl)
{ // vérification du type
string nom;
entree >> nom;
#ifdef MISE_AU_POINT
if (nom != "_plan_")
{ cout << "\nErreur, en lecture d'une instance Plan "
<< " on attendait _plan_ et on a lue: " << nom ;
cout << "istream & operator >> (istream & entree, Plan & pl)\n";
Sortie(1);
};
#endif
// puis lecture des différents éléments
entree >> nom >> pl.A >> nom >> pl.N;
return entree;
};
// surcharge de l'operateur d'ecriture
ostream & operator << ( ostream & sort,const Plan & pl)
{ // tout d'abord un indicateur donnant le type
sort << " _plan_ " ;
// puis les différents éléments
sort << "\n A= " << pl.A << " U= " << pl.N << " ";
return sort;
};

124
contact/Plan.h Normal file
View file

@ -0,0 +1,124 @@
// 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/>.
/************************************************************************
* DATE: 23/01/97 *
* $ *
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
* $ *
* PROJET: Herezh++ *
* $ *
************************************************************************
* BUT: Def de la geometrie d'un plan: un point et une normale. *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
* VERIFICATION: *
* *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* ! ! ! ! *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* MODIFICATIONS: *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* $ *
************************************************************************/
#ifndef PLAN_H
#define PLAN_H
#include "Droite.h"
#include "Coordonnee.h"
/// @addtogroup Groupe_sur_les_contacts
/// @{
///
class Plan
{ // surcharge de l'operator de lecture
friend istream & operator >> (istream &, Plan &);
// surcharge de l'operator d'ecriture
friend ostream & operator << (ostream &, const Plan &);
public :
// CONSTRUCTEURS :
// par defaut
Plan ();
// avec les datas
Plan ( const Coordonnee& B, const Coordonnee& vec);
// avec la dimension
Plan (int dim);
// de copie
Plan ( const Plan& a);
// DESTRUCTEUR :
~Plan ();
// surcharge des operator
Plan& operator = ( const Plan & P);
// METHODES PUBLIQUES :
// retourne le point de ref du plan
inline const Coordonnee& PointPlan() const { return A;};
// retourne la normale au plan
inline const Coordonnee& Vecplan() const { return N;};
// change la dimension de la droite
void Change_dim(int dima) {A.Change_dim(dima);N.Change_dim(dima);};
// change le point de ref du plan
void Change_ptref( const Coordonnee& B);
// change le vecteur normal du plan
void Change_normal( const Coordonnee& vec);
// change toutes les donnees
void change_donnees( const Coordonnee& B, const Coordonnee& vec);
// calcul l'intercection M d'une droite avec le plan, ramene 0 s'il n'y
// a pas d'intersection, ramene -1 si l'intercection ne peut pas etre calculee
// et 1 s'il y a un point d'intercection
int Intersection( const Droite & D,Coordonnee& M)const;
// calcul la distance d'un point à la droite
double Distance_au_plan(const Coordonnee& M) const;
// calcul du projeté d'un point sur le plan
Coordonnee Projete(const Coordonnee& M) const
{return (M - ((M-A)*N)*N);};
// ramène true si les deux points sont du même coté du plan, false sinon
bool DuMemeCote(const Coordonnee& M1, const Coordonnee& M2) const;
protected :
// VARIABLES PROTEGEES :
Coordonnee A; // un point du plan
Coordonnee N; // normale au plan
// METHODES PROTEGEES :
};
/// @} // end of group
#endif

202
contact/Sphere.cc Normal file
View file

@ -0,0 +1,202 @@
// 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 "Sphere.h"
#include "MathUtil.h"
#include "ParaGlob.h"
// CONSTRUCTEURS :
// par defaut
Sphere::Sphere () :
centre(3),rayon(0.)
{ if (ParaGlob::Dimension()!= 3)
{ cout << "\nErreur : la dimension n'est pas 3 et on veut utiliser une sphere ??? ";
cout << "\nSphere::Sphere ()" << endl;
Sortie(1);
};
};
// avec les datas
Sphere::Sphere ( const Coordonnee& B, const double& r):
centre(B),rayon(r)
{ if (ParaGlob::Dimension()!= 3)
{ cout << "\nErreur : la dimension n'est pas 3 et on veut utiliser une sphere ??? ";
cout << "\nSphere::Sphere (.." << endl;
Sortie(1);
};
if (rayon < 0.)
{ cout << "\n erreur dans la construction d'une sphere, le rayon " << rayon
<< " est negatif !! ";
Sortie(1);
};
};
// avec la dimension
Sphere::Sphere (int dim):
centre(dim),rayon(0.)
{ if (ParaGlob::Dimension()!= 3)
{ cout << "\nErreur : la dimension n'est pas 3 et on veut utiliser une sphere ??? ";
cout << "\nSphere::Sphere (.." << endl;
Sortie(1);
};
if (dim != 3)
{ cout << "\nErreur : la dimension demandee n'est pas 3 et on veut utiliser une sphere ??? ";
cout << "\nSphere::Sphere (.." << endl;
Sortie(1);
};
};
// de copie
Sphere::Sphere ( const Sphere& a):
centre(a.centre),rayon(a.rayon)
{};
// DESTRUCTEUR :
Sphere::~Sphere () {};
// surcharge des operator
Sphere& Sphere::operator = ( const Sphere & sph)
{ centre = sph.centre; rayon = sph.rayon; return *this; };
// METHODES PUBLIQUES :
// change le centre de la sphere
void Sphere::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 << "\nSphere::Change_centre( const Coordonnee& B)" << endl;
Sortie(1);
};
#endif
centre = B;
};
// change le rayon de la sphere
void Sphere::Change_rayon( const double & r)
{if (r < 0.)
{ cout << "\n erreur dans la construction d'une sphere, le rayon " << r
<< " est negatif !! "
<< "\n Sphere::Change_rayon( const double & r) ";
};
rayon = r;
};
// change toutes les donnees
void Sphere::change_donnees( const Coordonnee& B, const double& r)
{
#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 << "\nSphere::change_donnees( const Coordonnee& B, const double& r)" << endl;
Sortie(1);
};
#endif
centre = B;
if (r < 0.)
{ cout << "\n erreur dans la construction d'une sphere, le rayon " << r
<< " est negatif !! "
<< "\n Sphere::Change_rayon( const double & r) ";
Sortie(1);
};
rayon = r;
};
// calcul les deux intercections M1 et M2 d'une droite avec la sphere,
// 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 Sphere::Intersection( const Droite & dr,Coordonnee& M1, Coordonnee& M2)
{ // C le centre de la sphère, 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 = dr.PointDroite() - 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;
};
// calcul la distance d'un point à la sphere
double Sphere::Distance_a_sphere(const Coordonnee& M) const
{ return Dabs((centre - M).Norme()-rayon);
};
// projection d'un point M sur la parois de la sphère
// dans le cas où M = le centre de la sphère, la projection n'est pas
// possible, dans ce cas projection_ok = false en retour, sinon true
Coordonnee Sphere::Projete(const Coordonnee& M,bool& projection_ok) const
{ // on commence par calculer la projection du point sur l'axe
Coordonnee CM = M-centre;
double dist_centre = CM.Norme();
if (dist_centre < ConstMath::pasmalpetit)
{projection_ok = false;// on signale que la projection n'est pas possible
Coordonnee toto; // on fait l'opération en deux temps pour linux
return (toto);// on ramène un point par défaut
}
else
{projection_ok = true;
Coordonnee P = centre + CM * (rayon / dist_centre );
return P;
};
};
// surcharge de l'operateur de lecture
istream & operator >> (istream & entree, Sphere & sph)
{ // vérification du type
string nom;
entree >> nom;
#ifdef MISE_AU_POINT
if (nom != "_sphere_")
{ cout << "\nErreur, en lecture d'une instance sphere "
<< " on attendait _sphere_ et on a lue: " << nom ;
cout << "istream & operator >> (istream & entree, Sphere & sph)\n";
Sortie(1);
};
#endif
// puis lecture des différents éléments
entree >> nom >> sph.centre >> nom >> sph.rayon;
return entree;
};
// surcharge de l'operateur d'ecriture
ostream & operator << ( ostream & sort,const Sphere & sph)
{ // tout d'abord un indicateur donnant le type
sort << " _sphere_ " ;
// puis les différents éléments
sort << "\n C= " << sph.centre << " r= " << sph.rayon << " ";
return sort;
};

126
contact/Sphere.h Normal file
View file

@ -0,0 +1,126 @@
// 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/>.
/************************************************************************
* DATE: 19/01/2007 *
* $ *
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
* $ *
* PROJET: Herezh++ *
* $ *
************************************************************************
* BUT: Def géométrie d'une sphere: un centre et un rayon *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
* VERIFICATION: *
* *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* ! ! ! ! *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* MODIFICATIONS: *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* $ *
************************************************************************/
#ifndef SPHEREE_H
#define SPHEREE_H
#include "Droite.h"
#include "Coordonnee.h"
/// @addtogroup Groupe_sur_les_contacts
/// @{
///
class Sphere
{ // surcharge de l'operator de lecture
friend istream & operator >> (istream &, Sphere &);
// surcharge de l'operator d'ecriture
friend ostream & operator << (ostream &, const Sphere &);
public :
// CONSTRUCTEURS :
// par defaut
Sphere ();
// avec les datas
Sphere ( const Coordonnee& B, const double& r);
// avec la dimension
Sphere (int dim);
// de copie
Sphere ( const Sphere& a);
// DESTRUCTEUR :
~Sphere ();
// surcharge des operator
Sphere& operator = ( const Sphere & P);
// METHODES PUBLIQUES :
// retourne le centre de la Sphere
inline const Coordonnee& CentreSphere() const { return centre;};
// retourne le rayon de la sphere
inline double RayonSphere() const { return rayon;};
// change le centre de la sphere
void Change_centre( const Coordonnee& B);
// change le rayon de la sphere
void Change_rayon( const double & r);
// change toutes les donnees
void change_donnees( const Coordonnee& B, const double& r);
// calcul les deux intercections M1 et M2 d'une droite avec la sphere,
// ramene 0 s'il n'y a pas d'intersection, ramene -1 si l'intercection
// ne peut pas etre calculee
// et 1 s'il y a deux points d'intercection
int Intersection( const Droite & D,Coordonnee& M1, Coordonnee& M2);
// calcul la distance d'un point à la sphere
double Distance_a_sphere(const Coordonnee& M) const;
// ramène true si le point est à l'intérieur de la sphère, false sinon
bool Dedans(const Coordonnee& M)const
{ return ((centre - M).Norme() < rayon);};
// projection d'un point M sur la parois de la sphère
// dans le cas où M = le centre de la sphère, la projection n'est pas
// possible, dans ce cas projection_ok = false en retour, sinon true
Coordonnee Projete(const Coordonnee& M,bool& projection_ok) const;
protected :
// VARIABLES PROTEGEES :
Coordonnee centre; // centre de la sphere
double rayon; // rayon de la sphère
// METHODES PROTEGEES :
};
/// @} // end of group
#endif