248 lines
10 KiB
C++
248 lines
10 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/>.
|
|
|
|
//#include "Debug.h"
|
|
|
|
#include "GeomPentaCom.h"
|
|
#include <math.h>
|
|
#include "GeomSeg.h"
|
|
#include "GeomQuadrangle.h"
|
|
#include "GeomTriangle.h"
|
|
#include "MathUtil.h"
|
|
|
|
// constructeur
|
|
// le constructeur par défaut ne doit pas être utilisé
|
|
GeomPentaCom::GeomPentaCom()
|
|
{ cout << "\n erreur le constructeur par defaut ne doit pas être utilise " ;
|
|
cout << "\nGeomPentaCom::GeomPentaCom() " << endl;
|
|
Sortie(1);
|
|
};
|
|
|
|
// la dimension est 3, on a nbi pt d'integration, nbn noeuds et 5 faces, 9 aretes
|
|
GeomPentaCom::GeomPentaCom(int nbi , int nbn, Enum_interpol interpol) :
|
|
ElemGeomC0(3,nbi,nbn,5,9,PENTAEDRE,interpol)
|
|
{
|
|
// coordonnees des points d'integration
|
|
// ceux-ci sont obtenues par assemblage des points de surface et ceux d'épaisseur
|
|
int nbil,nbis;
|
|
switch (nbi)
|
|
{ case 1 : nbil = 1; nbis = 1; break;
|
|
case 2 : nbil = 2; nbis = 1; break;
|
|
case 3 : nbil = 3; nbis = 1; break;
|
|
case 6 : nbil = 2; nbis = 3; break;
|
|
case 8 : nbil = 2; nbis = 4; break;
|
|
case 9 : nbil = 3; nbis = 3; break;
|
|
case 12 : nbil = 3; nbis = 4; break;
|
|
case 18 : nbil = 3; nbis = 6; break;
|
|
default :
|
|
{cout << "\n erreur le pentaedre de nombre de point d'integration = " << nbi
|
|
<< "\n n\'est pas implante !! ";
|
|
cout << "\nGeomPentaCom::GeomPentaCom(int nbi) " << endl;
|
|
Sortie(1);
|
|
}
|
|
};
|
|
|
|
// un segment de base avec deux noeuds
|
|
GeomSeg segbase(nbil,2);
|
|
// un triangle de base avec 3 noeuds
|
|
// dans le cas de 3 pt d'integ de surface, on choisit des pt d'integ près des noeuds
|
|
// et non au milieu des arrêtes d'où nbis_tri pour ce choix
|
|
int nbis_tri(nbis); if (nbis == 3) nbis_tri = 1000+nbis;
|
|
GeomTriangle triabase(nbis_tri,3);
|
|
// on reconstitue les coordonnées en trois dimension
|
|
// et on calcul les poids d'intégration
|
|
int npi=1;
|
|
for (int niil=1;niil<=nbil;niil++)
|
|
for (int nis=1;nis<=nbis;nis++,npi++)
|
|
{ ptInteg(npi) = Coordonnee
|
|
((triabase.CoorPtInteg(nis))(1)
|
|
,(triabase.CoorPtInteg(nis))(2)
|
|
,(segbase.CoorPtInteg(niil))(1));
|
|
// les poids d'integration sont obtenus par produit de ceux de surface
|
|
// et ceux de hauteur
|
|
WI(npi)=triabase.Wi(nis) * segbase.Wi(niil);
|
|
};
|
|
};
|
|
|
|
// destructeur
|
|
GeomPentaCom::~GeomPentaCom()
|
|
{
|
|
};
|
|
// constructeur de copie
|
|
GeomPentaCom::GeomPentaCom(const GeomPentaCom& a) :
|
|
ElemGeomC0(a)
|
|
{
|
|
};
|
|
|
|
// en fonction de coordonnees locales, retourne true si le point est a l'interieur
|
|
// de l'element, false sinon
|
|
bool GeomPentaCom::Interieur(const Coordonnee& M)
|
|
{ if ((seg(1)->Interieur(M(3))) &&(face(1)->Interieur(Coordonnee(M(1),M(2)))))
|
|
return true;
|
|
else
|
|
return false;
|
|
};
|
|
|
|
// 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 GM est nul, dans ce cas retour de M
|
|
Coordonnee GeomPentaCom::Maxi_Coor_dans_directionGM(const Coordonnee& M)
|
|
{ // on décompose la recherche en segment et facette
|
|
double x=M(1); double y=M(2); double z=M(3);
|
|
// tout d'abord on calcul le point maxi dans le plan 1 2 et dans z indépendemment
|
|
Coordonnee XY(x,y); Coordonnee Z(z);
|
|
Coordonnee Pxy= face(1)->Maxi_Coor_dans_directionGM(XY);
|
|
Coordonnee Pz = seg(1)->Maxi_Coor_dans_directionGM(Z);
|
|
Coordonnee P(0.,0.,0.); // le retour
|
|
|
|
// G le centre de gravité: il a la cote nulle et il est centre de gravité du triangle médian
|
|
double untiers = 1./3.;
|
|
Coordonnee G(untiers,untiers,0.);
|
|
Coordonnee Gxy(untiers,untiers); // en 2D xy
|
|
|
|
// le pb est de savoir si on sort sur une des faces triangulaires, ou sur une des faces rectangulaires
|
|
// Si le point maxi P sort par une face quadrangulaire, cela signifie que z_P est compris entre 1 et -1
|
|
// --- on va donc commencer par calculer le z_P
|
|
// on considère les points: G (z_G=0), Pxy dont la côte est également nulle, M, et z_P ?
|
|
// soit un plan vertical (normale à Oxy) passant par G,M,Pxy. Il passe également par P. On note Mxy la projection
|
|
// de M sur le plan Oxy. Dans le plan verticale on peut faire une régle de trois:
|
|
// on a z_P/z_M = GPxy/GMxy d'où z_P = z_M * (GPxy/GMxy)
|
|
// NB: si ||GMxy|| = 0, là c'est sur ça ne sort pas par un quadrangle
|
|
|
|
|
|
// 1)-----
|
|
Coordonnee GMxy(x-untiers,y-untiers); // en 2D xy
|
|
int fin_traitement = false; // par défaut
|
|
double gmxy = GMxy.Norme();
|
|
if (gmxy > ConstMath::trespetit)
|
|
{ double z_P = z * ( (Pxy-Gxy).Norme()) / gmxy;
|
|
if (Abs(z_P) <= 1.) // là c'est ok on sort par un rectangle
|
|
{P(1)=Pxy(1); P(2)=Pxy(2); P(3)=z_P;
|
|
fin_traitement = true;// pour indiquer que c'est ok
|
|
};
|
|
};
|
|
// 2)----
|
|
// dans le cas où fin_traitement est toujours faux, cela signifie que le cas sortie par un quadrangle n'a pas marché
|
|
// donc cela signifie que le point P est sur la face sup ou inf. En fait cela dépend de z_M si c'est positif alors
|
|
// c'est la face sup donc z_P=1 sinon c'est la face inf donc z_P=-1
|
|
// on considère un plan vertical qui passe par G, M et P . On a GMxy/z_M = GPxy/z_P avec z_P connu
|
|
// d'où GPxy = z_P * GMxy/z_M = e (en module) d'où \vec{GPxy} = e * \vec{GMxy)/||GMxy|| = z_P / z_M * \vec{GMxy) d'où P
|
|
if (! fin_traitement)
|
|
{ if (z > ConstMath::trespetit)
|
|
{ // cas d'une sortie par la facette supérieure, z_P = 1.
|
|
Coordonnee GPxy = GMxy / z;
|
|
P(1) = untiers + GPxy(1); P(2) = untiers + GPxy(2); P(3) = 1.;
|
|
fin_traitement = true;// pour indiquer que c'est ok
|
|
}
|
|
else if (z < -ConstMath::trespetit)
|
|
{ // cas d'une sortie par la facette inférieure, z_P = -1.
|
|
Coordonnee GPxy = (-1./z) * GMxy ;
|
|
P(1) = untiers + GPxy(1); P(2) = untiers + GPxy(2); P(3) = -1.;
|
|
fin_traitement = true;// pour indiquer que c'est ok
|
|
}
|
|
else
|
|
// arrivée ici, cela veut dire que z est quasiment nulle, et pourtant due au 1)---
|
|
// - soit Abs(z_P) est > 1 , ce qui n'est pas possible car on vient de trouver z quasiment nulle
|
|
// - soit gmxy est également quasiment nulle : dans ce cas on choisit arbitrairement P(0,0,0)
|
|
{ P.Zero();};
|
|
};
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if ((P(1) + P(2)) > (1. + ConstMath::petit))
|
|
{ cout << "\n petit soucis, le point maxi calcule n'est pas correcte ";
|
|
P.Affiche(cout);
|
|
cout << "\n GeomPentaCom::Maxi_Coor_dans_directionGM(.." << endl;
|
|
};
|
|
#endif
|
|
/*
|
|
|
|
// en fait il y a un pb, car la suite à l'air correcte si on cherche l'intersection de OM avec la frontière
|
|
// mais ce qu'il nous faut c'est GM avec la frontière, il faut donc quelques modifs !!
|
|
cout << "\n ***************** GeomPentaCom::Maxi_Coor_dans_directionGM( a revoir ******" << endl;
|
|
|
|
if ((Dabs(x) <= ConstMath::petit) && (Dabs(y) <= ConstMath::petit))
|
|
{ // cas où x et y transformées sont quasi nulles, il ne reste plus qu'à examiner le cas de z
|
|
if (Dabs(z) <= ConstMath::petit)
|
|
{ // cas où toutes les coordonnées sont nulles, retour de M
|
|
// on met les coordonnées en x, y strictement à 0. pour éviter d'être légèrement à l'extérieur
|
|
return Coordonnee (0.,0.,z);
|
|
}
|
|
else // cas où l'on modifie suivant z (le point est sur le triangle sup ou inf)
|
|
{ return Coordonnee (0.,0.,Pz(1));
|
|
}
|
|
}
|
|
else // cas où x ou/et y sont non nulles
|
|
{// on calcul de facteur de réduction suivant xy
|
|
double facteurxy=0.;
|
|
if (Dabs(x) <= ConstMath::petit) // donc y est non nulle
|
|
{// cas donc où seule y a peut-être variée
|
|
facteurxy = Pxy(2)/y;
|
|
}
|
|
else if (Dabs(y) <= ConstMath::petit) // donc x est non nulle
|
|
{// cas donc où seule x a peut-être variée
|
|
facteurxy = Pxy(1)/x;
|
|
}
|
|
else // cas où x et y sont non nulles
|
|
{// cas donc où x et y x ont peut-être variées
|
|
facteurxy = MiN(Pxy(1)/x,Pxy(2)/y);
|
|
}
|
|
// on regarde maintenant du coté de z
|
|
if (Dabs(z) <= ConstMath::petit)
|
|
{ // z ne participe pas car nulle
|
|
return Coordonnee (Pxy(1),Pxy(2),z);
|
|
}
|
|
else // cas où la modification est possible suivant z
|
|
{ double facteurz= Pz(1)/z;
|
|
// on prend le plus petit des deux facteurs
|
|
if (facteurxy <= facteurz)
|
|
{ // c'est le facteur dans le plan qui agit
|
|
return Coordonnee (Pxy(1),Pxy(2),facteurxy*z);
|
|
}
|
|
else
|
|
{ // c'est le facteur suivant z qui agit
|
|
return Coordonnee (facteurz*x,facteurz*y,Pz(1));
|
|
}
|
|
}
|
|
} //-- fin du cas où x ou/et y sont non nulles
|
|
// cas qui ne doit jamais arriver car on se situe forcément dans un des trois secteurs
|
|
// précédemment étudiés
|
|
cout << "\n erreur inconnue !!"
|
|
<< "\n GeomPentaCom::Maxi_Coor_dans_directionGM(..";
|
|
Sortie(1); */
|
|
return M; // pour éviter un warning
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|