357 lines
20 KiB
C++
357 lines
20 KiB
C++
|
|
// 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-2022 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 geometrique generique pour une discrétisation C0. *
|
|
* $ *
|
|
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
|
|
* VERIFICATION: *
|
|
* *
|
|
* ! date ! auteur ! but ! *
|
|
* ------------------------------------------------------------ *
|
|
* ! ! ! ! *
|
|
* $ *
|
|
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
|
|
* MODIFICATIONS: *
|
|
* ! date ! auteur ! but ! *
|
|
* ------------------------------------------------------------ *
|
|
* $ *
|
|
************************************************************************/
|
|
#ifndef ELEMGEOMC0_H
|
|
#define ELEMGEOMC0_H
|
|
|
|
// #include <bool.h>
|
|
#include "Mat_pleine.h"
|
|
#include "Tableau_T.h"
|
|
#include "Vecteur.h"
|
|
#include "Coordonnee.h"
|
|
#include "Enum_interpol.h"
|
|
#include "Enum_geom.h"
|
|
#include "Enum_type_pt_integ.h"
|
|
|
|
/** @defgroup Les_Elements_de_geometrie
|
|
*
|
|
* BUT: groupe concernant les éléments de géométrie 1D, 2D, 3D
|
|
*
|
|
* \author Gérard Rio
|
|
* \version 1.0
|
|
* \date 23/01/97
|
|
* \brief groupe concernant les éléments de géométrie 1D, 2D, 3D
|
|
*
|
|
*/
|
|
|
|
|
|
/// @addtogroup Les_Elements_de_geometrie
|
|
/// @{
|
|
///
|
|
|
|
class ElemGeomC0
|
|
{
|
|
public :
|
|
// CONSTRUCTEURS :
|
|
ElemGeomC0(); // pardefaut
|
|
// cas ou l'on connait la dimension = dim,
|
|
// le nombre de point d'integration nbi et le nombre
|
|
// de noeud de l'element nbne
|
|
// le nombre de face nbfe, le nombre de segment nbse
|
|
|
|
// IMPORTANT : par defaut le nombre de pt d'integ pour les faces ou segents
|
|
// s'il y en a, est calcule a partir du nombre total de point d'integ
|
|
// se referer a la description des elements
|
|
// par défaut le type de point d'intégration est de gauss
|
|
|
|
ElemGeomC0(int dim,int nbi,int nbne,int nbfe,int nbse
|
|
,Enum_geom geom, Enum_interpol interpol
|
|
,Enum_type_pt_integ type_pti = PTI_GAUSS);
|
|
// de copie
|
|
ElemGeomC0(const ElemGeomC0& a);
|
|
// DESTRUCTEUR :
|
|
virtual ~ElemGeomC0();
|
|
// METHODES PUBLIQUES :
|
|
// retourne la dimension de l'element
|
|
inline int Dimension() const { return dimension;};
|
|
// retourne le nombre total de point d'integration
|
|
inline int Nbi() const { return tabPhi.Taille();};
|
|
// retourne le nombre total de noeud
|
|
inline int Nbne() const { return NBNE;};
|
|
// retourne le nombre de face
|
|
inline int NbFe() const { return NBFE;};
|
|
// retourne le nombre de segment
|
|
inline int NbSe() const { return NBSE;};
|
|
// retourne les coordonnees du point d'integration i
|
|
inline Coordonnee const & CoorPtInteg(int i) const { return ptInteg(i);};
|
|
// retourne les fonctions d'interpolation au point d'integration i
|
|
inline Vecteur const & Phi(int i) const { return tabPhi(i);};
|
|
// retourne les tableau de fonctions d'interpolation
|
|
inline Tableau <Vecteur> const & TaPhi() const { return tabPhi;};
|
|
// retourne les derivees des fonctions d'interpolation au point d'integration i
|
|
inline Mat_pleine const& Dphi(int i) { return tabDPhi(i);};
|
|
// retourne le tableau des derivees des fonctions d'interpolation
|
|
inline Tableau < Mat_pleine > const& TaDphi() const { return tabDPhi;};
|
|
// retourne les poids d'integration du point d'integration i
|
|
inline double Wi(int i) const { return WI(i);};
|
|
// retourne le vecteur des poids d'integration
|
|
inline Vecteur const & TaWi() const { return WI;};
|
|
// retourne l'"element face de référence" correspondant a la face i
|
|
inline ElemGeomC0 const & ElemFace(int i) const { return *face(i);};
|
|
// retourne l'"element segment de référence" correspondant a l'arrète i
|
|
inline ElemGeomC0 const & ElemSeg(int i) const { return *seg(i);};
|
|
// retourne true si l'element est complet sinon false
|
|
bool Complet() const ;
|
|
// retourne la connection des noeuds des faces par rapport a ceux de l'element
|
|
inline Tableau<Tableau<int> > const & Nonf() const { return NONF;};
|
|
// retourne la connection des noeuds des arêtes par rapport a ceux de l'element
|
|
inline Tableau<Tableau<int> > const & NonS() const { return NONS;};
|
|
// retourne les coordonnées des noeuds locaux
|
|
inline Tableau <Coordonnee > const & PtelemRef() const {return ptelem;};
|
|
// retourne le type de géométrie
|
|
inline Enum_geom TypeGeometrie() const {return id_geom;};
|
|
// retourne le type d'interpolation
|
|
inline Enum_interpol TypeInterpolation() const {return id_interpol;};
|
|
// retourne le type de points d'intégration
|
|
inline Enum_type_pt_integ TypePointIntegration() const {return id_type_pt_integ;};
|
|
|
|
// cas de l'extrapolation de grandeur des points d'intégrations aux noeuds
|
|
// ramène un tableau de pondération tab(i,j) qu'il faut appliquer
|
|
// aux noeuds pour avoir la valeur aux noeuds
|
|
// val_au_noeud(i) = somme_(de j=indir(i)(1) à indir(i)(taille(indir(i)) )) {tab(i)(j) * val_pt_integ(j) }
|
|
// cas = 1: la valeur au noeud = la valeur au pt d'integ le plus près ou une moyenne des
|
|
// pt les plus près (si le nb de pt d'integ < nb noeud)
|
|
class ConteneurExtrapolation
|
|
{ // surcharge de l'operator de lecture avec le type
|
|
friend istream & operator >> (istream &, ConteneurExtrapolation &);
|
|
// surcharge de l'operator d'ecriture
|
|
friend ostream & operator << (ostream &, const ConteneurExtrapolation &);
|
|
public:
|
|
ConteneurExtrapolation() : tab(),indir() {};
|
|
ConteneurExtrapolation(const ConteneurExtrapolation &a) :
|
|
tab(a.tab),indir(a.indir) {};
|
|
ConteneurExtrapolation& operator= ( const ConteneurExtrapolation& a)
|
|
{tab=a.tab;indir=a.indir; return *this;};
|
|
// les tableaux publiques
|
|
Tableau<Tableau<double> > tab;
|
|
Tableau<Tableau<int> > indir;
|
|
};
|
|
inline ConteneurExtrapolation const & ExtrapolationNoeud(int cas) const { return extrapol(cas);};
|
|
|
|
// création d'élément identiques : cette fonction est analogue à la fonction new
|
|
// elle y fait d'ailleurs appel. l'implantation est spécifique dans chaque classe
|
|
// dérivée
|
|
// pt est le pointeur qui est affecté par la fonction
|
|
virtual ElemGeomC0 * newElemGeomC0(ElemGeomC0 * pt) = 0;
|
|
|
|
// creation d'une numérotation correspondant à un élément d'orientation inverse
|
|
// la numérotation est donnée par rapport à la numérotation normale de l'élément (de 1 à NBNE)
|
|
inline Tableau<int> const & InvConnec() const { return INVCONNEC;};
|
|
|
|
// ramène le tableau des tranches de numérotation
|
|
// c-a-d: ex: un quadrangle quadratique à 9 noeuds à 3 tranches: 4, 4, 1
|
|
// - la première tranche correspond aux 4 noeuds sommets, la numérotation peut être cyclique (permutation circulaire)
|
|
// parmi ses 4 noeuds
|
|
// - la second tranche correspond aux 4 noeuds intermédiaires, là aussi, la numérotation peut être une
|
|
// permutation circulaire parmi ses 4 noeuds
|
|
// - la troisième tranche correspond au noeud central
|
|
// NB: la somme de toutes les tranches = le nombre total de noeud de l'élément
|
|
inline Tableau <int> const & Ind() const { return IND;};
|
|
|
|
// ramène les tableaux de permutations permettant de calculer le tableau
|
|
// de connection permuté, ceci par rapport à la numérotation normale
|
|
// le second tableau correspond aux noeuds sommets (angle par exemple), les autres tableaux sont
|
|
// les noeuds internes aux arrêtes ou à l'élément
|
|
// chaque tableau contient le double du nombre de noeud associé,
|
|
// Utilisation: soit un noeud de numéro local nl, soit nl <= permut(2).Taille()= t1, dans ce cas il s'agit du
|
|
// premier tableau : la numérotation sera pour i=1 à t1: permut(2)(nl+i-1)
|
|
// pour les numéros au-dessus de t1, on décale de permut(1)(1) = d1, c-a-d si permut(3).Taille()=t2
|
|
// la numérotation qui suit est : pour j=1 à t2 : permut(3)(
|
|
// en chantier inline Tableau<Tableau<int> > const & Permutation() const { return permut;};
|
|
|
|
//--------- cas de coordonnees locales quelconques ----------------
|
|
// retourne les fonctions d'interpolation au point M (en coordonnees locales)
|
|
virtual const Vecteur& Phi_point(const Coordonnee& M) = 0;
|
|
// retourne les derivees des fonctions d'interpolation au point M (en coordonnees locales)
|
|
virtual const Mat_pleine& Dphi_point(const Coordonnee& M) = 0;
|
|
// en fonction de coordonnees locales, retourne true si le point est a l'interieur
|
|
// de l'element, false sinon
|
|
virtual bool Interieur(const Coordonnee& M) = 0;
|
|
// en fonction de coordonnees locales, retourne le point local P, maximum intérieur à l'élément, donc sur la frontière
|
|
// dont les coordonnées sont sur la droite GM: c-a-d GP = alpha GM, avec apha maxi et P appartenant à la frontière
|
|
// de l'élément, G étant le centre de gravité, sauf si GGM est nul, dans ce cas retour de M
|
|
virtual Coordonnee Maxi_Coor_dans_directionGM(const Coordonnee& M)
|
|
{ cout << "\n Pas implante pour l'instant, Maxi_Coor_dans_direction(..."; Sortie(1);
|
|
Coordonnee P=M;
|
|
return P;}; // =0; pour l'instant essaie
|
|
|
|
//--------- triangulation des faces lorsqu'elles existes -------
|
|
// le retour est la liste de la connection par rapport à la numérotation locale
|
|
// de tous les facettes triangulaire linéaire qui composent les faces de l'élément
|
|
// l'indice de premier niveau indique le numéro de la face
|
|
// l'indice de deuxième niveau correspond à la numérotation du triangle
|
|
// le troisième indice qui varie de 1 à 3 indique la connection
|
|
// tab(i)(j)(k) : pour la face i, le triangle j , le numéro du noeud dans la
|
|
// numérotation locale des noeuds de l'élément
|
|
const Tableau<Tableau<Tableau<int> > >& Trian_lin() const { return NONFt; };
|
|
|
|
//--------- segmentation des arêtes -------
|
|
// le retour est la liste de la connection par rapport à la numérotation locale
|
|
// de tous les segments linéaire qui composent les arêtes de l'élément
|
|
// l'indice de premier niveau indique le numéro de l'arête
|
|
// l'indice de deuxième niveau correspond à la numérotation du segment
|
|
// le troisième indice qui varie de 1 à 2 indique la connection
|
|
// tab(i)(j)(k) : pour l'arête i, le segment j , le numéro du noeud dans la
|
|
// numérotation locale des noeuds de l'élément
|
|
const Tableau<Tableau<Tableau<int> > >& Trian_seg() const { return NONSs; };
|
|
|
|
// fonctions utilitaires spécifiques à l'interpolation linéaire:
|
|
// méthode utilisée pour extrapoler une grandeur à partir:
|
|
// soit de 2 points -> ex: pour une extrapolation linéaire (métrique en 1D)
|
|
// soit de 3 points -> ex: pour une extrapolation linéaire (métrique en 2D )
|
|
// soit de 4 points (non coplanaires) -> ex: pour une extrapolation linéaire (métrique en 3D )
|
|
// ces nipt points sont définit par le tableau tab_M,
|
|
//et en sortie on récupère les dim vecteurs de la base naturelle et les dim vecteurs de la base duale
|
|
// dans le cas où le calcul n'est pas possible (points trop près) on ramène false, sinon true
|
|
static bool Bases_naturelles_duales(const Tableau <Coordonnee >& tab_M
|
|
,Tableau <Coordonnee> & giB,Tableau <Coordonnee> & giH );
|
|
// on donne les vecteurs de la base naturelle, et les vecteurs de la base duale
|
|
// l'origine de la base O,
|
|
// et un point A,
|
|
// en sortie: les coordonnées locale de A dans la base naturelle, et les fonctions d'interpolation
|
|
// linéaire (1,2 ou 3D ce qui correspond à ligne, triangle, tetraedre) au point A
|
|
static void Coor_phi(const Coordonnee& O
|
|
,const Tableau <Coordonnee> & giH_, const Coordonnee& A
|
|
,Vecteur& phi, Coordonnee& theta) ;
|
|
|
|
protected :
|
|
// VARIABLES PROTEGEES :
|
|
int dimension; // dimension de l'element
|
|
int NBNE; // nombre de noeuds de l'element
|
|
int NBFE; // nombre de face par element
|
|
int NBSE; // nombre de segment par element
|
|
Tableau<Tableau<int> > NONF; // connexion des noeuds des faces par rapport a ceux des elements
|
|
// NONF(j)(i) noeud i de la face j = noeud de l'element
|
|
Tableau<Tableau<int> > NONS; // connexion des noeuds des segments par rapport a ceux des elements
|
|
// NONS(j)(i) noeud i du segment j = noeud de l'element
|
|
Tableau<int> INVCONNEC; // conection d'un élément d'orientation inverse du sens normal,
|
|
// la numérotation est donnée par rapport à la numérotation normale (de 1 à NBNE)
|
|
Tableau <int> IND; //tableau donnant les tranches de numérotation qui ont une numérotation cyclique
|
|
// c-a-d: ex: un quadrangle quadratique à 9 noeuds à 3 tranches: 4, 4, 1
|
|
// - la première tranche correspond aux 4 noeuds sommets, la numérotation peut être cyclique (permutation circulaire)
|
|
// parmi ses 4 noeuds
|
|
// - la second tranche correspond aux 4 noeuds intermédiaires, là aussi, la numérotation peut être une
|
|
// permutation circulaire parmi ses 4 noeuds
|
|
// - la troisième tranche correspond au noeud central
|
|
// NB: la somme de toutes les tranches = le nombre total de noeud de l'élément
|
|
|
|
|
|
// Tableau<Tableau<int> > permut; // définit les tableaux de permutations permettant de calculer le tableau
|
|
// de connection permuté, ceci par rapport à la numérotation normale
|
|
// le premier tableau correspond aux noeuds sommets (angle par exemple), les autres tableaux sont
|
|
// les noeuds internes aux arrêtes ou à l'élément
|
|
// en chantier !!!!!
|
|
// ramène les tableaux de permutations permettant de calculer le tableau
|
|
// de connection permuté, ceci par rapport à la numérotation normale
|
|
// le second tableau correspond aux noeuds sommets (angle par exemple), les autres tableaux sont
|
|
// les noeuds internes aux arrêtes ou à l'élément
|
|
// chaque tableau contient le double du nombre de noeud associé,
|
|
// Utilisation: soit un noeud de numéro local nl, soit nl <= permut(2).Taille()= t1, dans ce cas il s'agit du
|
|
// premier tableau : la numérotation sera pour i=1 à t1: permut(2)(nl+i-1)
|
|
// pour les numéros au-dessus de t1, on décale de permut(1)(1) = d1, c-a-d si permut(3).Taille()=t2
|
|
// la numérotation qui suit est : pour j=1 à t2 : permut(3)(
|
|
|
|
// cas de l'extrapolation de grandeur des points d'intégrations aux noeuds
|
|
// def du tableau de pondération tab(i)(j) qu'il faut appliquer
|
|
// aux noeuds pour avoir la valeur aux noeuds
|
|
// val_au_noeud(i) = somme_(de j=borne_inf(i) à borne_sup(i)) {tab(i)(j) * val_pt_integ(j) }
|
|
// cas = 1: la valeur au noeud = la valeur au pt d'integ le plus près ou une moyenne des
|
|
// pt les plus près (si le nb de pt d'integ < nb noeud)
|
|
Tableau <ConteneurExtrapolation > extrapol;
|
|
|
|
// pour l'element
|
|
Tableau <Coordonnee > ptelem ; // coordonnees des points de l'élément
|
|
Tableau <Coordonnee > ptInteg ; // coordonnees des points d'integration
|
|
Tableau <Vecteur> tabPhi ; // tabPhi(ni) = phi = fonctions d'interpolation
|
|
// au point d'interpolation ni, phi a la dimension de nbne
|
|
Tableau < Mat_pleine > tabDPhi; // tabDPhi(ni) = Dphi tel que, Dphi(i,r) =
|
|
// valeur de la derivee de la fonction phi(r) par rapport a la coordonnee`
|
|
// locale i (= 1 ou 2 ou 3, ceci dependant de la dimension de l'element)
|
|
Vecteur WI; // poids d'integration
|
|
// pour les faces
|
|
Tableau < ElemGeomC0 * > face ; // pointe sur les elements faces
|
|
// pour les segments
|
|
Tableau < ElemGeomC0 * > seg ; // pointe sur les elements segments
|
|
// type de géométrie
|
|
Enum_geom id_geom;
|
|
// type d'interpolation
|
|
Enum_interpol id_interpol;
|
|
// type de point d'intégration
|
|
Enum_type_pt_integ id_type_pt_integ;
|
|
|
|
// pour la triangulation linéaires des faces lorsqu'elles existent
|
|
// connexion des noeuds des faces par rapport a ceux des elements
|
|
// l'indice de premier niveau indique le numéro de la face
|
|
// l'indice de deuxième niveau correspond à la numérotation du triangle
|
|
// le troisième indice qui varie de 1 à 3 indique la connection
|
|
Tableau<Tableau<Tableau<int> > > NONFt;
|
|
|
|
// pour la segmentation linéaires des arêtes lorsqu'elles existent
|
|
// connexion des noeuds des arêtes par rapport a ceux des elements
|
|
// l'indice de premier niveau indique le numéro de l'arête
|
|
// l'indice de deuxième niveau correspond à la numérotation du segment
|
|
// le troisième indice qui varie de 1 à 2 indique la connection
|
|
Tableau<Tableau<Tableau<int> > > NONSs;
|
|
|
|
// METHODES PROTEGEES :
|
|
|
|
// -- méthode identique à Bases_naturelles_duales sauf que l'on utilise une numérotation
|
|
// indirecte ce qui permet d'éviter de contruire le tableau de coordonnées lorsque celui-ci
|
|
// existe globalement
|
|
|
|
// méthode utilisée pour extrapoler une grandeur à partir:
|
|
// soit de 2 points -> ex: pour une extrapolation linéaire (métrique en 1D)
|
|
// soit de 3 points -> ex: pour une extrapolation linéaire (métrique en 2D )
|
|
// soit de 4 points cas = 2 : extrapolation bi-linéaire (métrique en 2D)
|
|
// soit de 4 points (non coplanaires) cas = 1 -> ex: pour une extrapolation linéaire (métrique en 3D )
|
|
// soit de 8 points cas = 1 : pour une extrapolation bi-linéaire (métrique en 3D)
|
|
// ces nipt points d'intégration, dont les numéros dans ptInteg sont défini par le tableau indirec(i),
|
|
// i=1 à nipt, avec ptInteg(indirec(i)) le point à considérer, et en sortie on récupère les dim vecteurs
|
|
// de la base naturelle et les dim vecteurs de la base duale
|
|
// dans le cas où le calcul n'est pas possible (points trop près) on ramène false, sinon true
|
|
bool Bases_naturel_duales(const Tableau <int>& indirec
|
|
,Tableau <Coordonnee> & giB,Tableau <Coordonnee> & giH_
|
|
, int cas = 1) const;
|
|
};
|
|
/// @} // end of group
|
|
|
|
#endif
|