// 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) .
//
// 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 .
//
// For more information, please consult: .
#include"ElemGeomC0.h"
#include "MathUtil.h"
// CONSTRUCTEURS :
// pardefaut
ElemGeomC0::ElemGeomC0() :
ptelem(),ptInteg() ,tabPhi() , tabDPhi () , WI()
,NONF(),NONS(),face(), seg(),NONFt(),NONSs(),INVCONNEC(),extrapol(1), IND(0)
,id_type_pt_integ(PTI_GAUSS)
{ NBNE = 0; // pas de noeud
dimension = 0; // pas de dimension
NBFE = 0; // pas de face
NBSE = 0; // pas de segment
id_interpol=RIEN_INTERPOL;
id_geom=RIEN_GEOM;
};
// cas ou l'on connait nbi et nbne
ElemGeomC0::ElemGeomC0(int dim,int nbi,int nbne,int nbfe,int nbse
,Enum_geom geom, Enum_interpol interpol
,Enum_type_pt_integ type_pti ) :
dimension(dim),NBNE(nbne),NBFE(nbfe), NBSE(nbse),NONFt(nbfe),NONSs(nbse)
,ptelem(nbne),ptInteg(nbi) ,tabPhi(nbi) , tabDPhi (nbi) , WI(nbi)
,NONF(nbfe),NONS(nbse),face(nbfe), seg(nbse),id_interpol(interpol),id_geom(geom)
,INVCONNEC(nbne),extrapol(1),IND(0),id_type_pt_integ(type_pti)
{ for (int i= 1; i<=nbi; i++)
{ tabPhi(i).Change_taille(NBNE);
tabDPhi(i).Initialise (dimension,NBNE,0.);
ptInteg(i).Change_dim(dimension);
};
for (int i1 =1; i1<= NBFE; i1++)
face(i1) = NULL;
for (int i2 =1; i2<= NBSE; i2++)
seg(i2) = NULL;
extrapol(1).tab.Change_taille(NBNE);
extrapol(1).indir.Change_taille(NBNE);
for (int ne=1;ne<=NBNE;ne++)
{extrapol(1).tab(ne).Change_taille(nbi);
};
};
// de copie
ElemGeomC0::ElemGeomC0(const ElemGeomC0& a) :
dimension(a.dimension),NBNE(a.NBNE),NBFE(a.NBFE),NBSE(a.NBSE)
,ptelem(a.ptelem) ,ptInteg(a.ptInteg) ,tabPhi(a.tabPhi)
, tabDPhi (a.tabDPhi) , WI(a.WI)
,NONF(a.NONF),NONS(a.NONS),face((a.face).Taille()),seg((a.seg).Taille())
,id_interpol(a.id_interpol),id_geom(a.id_geom),NONFt(a.NONFt),NONSs(a.NONSs)
,INVCONNEC(a.INVCONNEC),extrapol(a.extrapol),IND(a.IND),id_type_pt_integ(a.id_type_pt_integ)
{ // concernant les faces et les segments éventuelles il faut créer ces éléments
// int facetaille= face.Taille();
// for (int i=1;i<=facetaille;i++)
// if (a.face(i) != NULL) face(i)= (a.face(i))->newElemGeomC0(face(i));
// int segtaille= seg.Taille();
// for (int i=1;i<=segtaille;i++)
// if (a.seg(i) != NULL) seg(i)=(a.seg(i))->newElemGeomC0(seg(i));
};
// DESTRUCTEUR :
ElemGeomC0::~ElemGeomC0()
{ // il faut supprimer les faces ou segments pointés, eventuelle
// --- en fait due au type de construction des faces et segments, ---
// ceux-ci sont construit dans l'élément dérivé spécifique et donc supprimé dans
// également dans les éléments spécifiques
// int facetaille= face.Taille();
// for (int i=1;i<=facetaille;i++)
// if (face(i) != NULL) delete face(i);
// int segtaille= seg.Taille();
// for (int i=1;i<=segtaille;i++)
// if (seg(i) != NULL) delete seg(i);
};
// 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
bool ElemGeomC0::Bases_naturelles_duales(const Tableau & tab_M
,Tableau & gi_B,Tableau & gi_H )
{// on récupère la dimension de la métrique
int dime = tab_M.Taille();
// traitement suivant la taille
switch (dime)
{ case 1:
{// calcul de la base naturelle
gi_B.Change_taille(1);
// gi_B(1)= tab_M(2)-tab_M(1); // modif 29 oct 2020
gi_B(1)= 0.5*(tab_M(2)-tab_M(1));
// calcul de la métrique
double g11BB=gi_B(1)*gi_B(1);
// inverse de la métrique
double det = g11BB;
if (Dabs(det) <= ConstMath::trespetit)
{ cout << "\nAttention : le determinant du systeme est nulle !"
<< " il n'est pas possible de calculer la base duale ";
cout << "\n ElemGeomC0::Bases_naturelles_duales(..";
return false;
};
double unsurdet = 1./det;
double g11HH = unsurdet * g11BB;
// base duale
gi_H(1)= g11HH * gi_B(1);
}
break;
case 2:
{// calcul de la base naturelle
gi_B.Change_taille(2);
gi_B(1)= tab_M(2)-tab_M(1);
gi_B(2)= tab_M(3)-tab_M(1);
// calcul de la métrique
double g11BB=gi_B(1)*gi_B(1);
double g12BB=gi_B(1)*gi_B(2);
double g22BB=gi_B(2)*gi_B(2);
// inverse de la métrique
double det = g11BB*g22BB - g12BB*g12BB;
if (Dabs(det) <= ConstMath::trespetit)
{ cout << "\nAttention : le determinant du systeme est nulle !"
<< " il n'est pas possible de calculer la base duale ";
cout << "\n ElemGeomC0::Bases_naturelles_duales(..";
return false;
};
double unsurdet = 1./det;
double g11HH = unsurdet * g22BB;
double g22HH = unsurdet * g11BB;
double g12HH = - unsurdet * g12BB;
// base duale
gi_H(1)= g11HH * gi_B(1) + g12HH * gi_B(2);
gi_H(2)= g12HH * gi_B(1) + g22HH * gi_B(2);
}
break;
case 3:
{// calcul de la base naturelle
gi_B.Change_taille(3);
gi_B(1)= tab_M(2)-tab_M(1);
gi_B(2)= tab_M(3)-tab_M(1);
gi_B(3)= tab_M(4)-tab_M(1);
// calcul de la métrique
Mat_pleine gijBB(3,3);
for (int i=1;i<=3;i++)
for (int j=1;j<=i;j++)
gijBB(i,j)=gijBB(j,i)=gi_B(i)*gi_B(j);
// inverse de la métrique
Mat_pleine gijHH(3,3); // def
try
{gijHH = gijBB.Inverse();}
catch (ErrResolve_system_lineaire excep)
// cas d'une erreur pendant l'inversion
{ cout << "\nAttention : l'extrapolation via un tetraedre lineaire !"
<< " n'est pas possible ";
cout << "\n ElemGeomC0::Bases_naturelles_duales(..";
return false;
};
// base duale
gi_H(1)= gijHH(1,1) * gi_B(1) + gijHH(1,2) * gi_B(2) + gijHH(1,3) * gi_B(3);
gi_H(2)= gijHH(2,1) * gi_B(1) + gijHH(2,2) * gi_B(2) + gijHH(2,3) * gi_B(3);
gi_H(3)= gijHH(3,1) * gi_B(1) + gijHH(3,2) * gi_B(2) + gijHH(3,3) * gi_B(3);
}
break;
default:
cout << "\n erreur !!! cas dime = " << dime << " non prevu "
<< "\n ElemGeomC0::Bases_naturel_duales(...";
Sortie(1);
};
// retour ok
return true;
};
// 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
void ElemGeomC0::Coor_phi(const Coordonnee& O
,const Tableau & giH_, const Coordonnee& A
,Vecteur& phi_, Coordonnee& theta)
{// on récupère la dimension de la métrique
int dime = giH_.Taille();
// vecteur local
Coordonnee OA = A - O;
// traitement suivant la taille
switch (dime)
{ case 1:
{ // calcul des coordonnées locales
theta.Change_dim(1);
theta(1) = OA * giH_(1);
// les fonctions d'interpolation
phi_.Change_taille(2);
phi_(1)=0.5*(1.-theta(1));
phi_(2)=0.5*(1.+theta(1));;
}
break;
case 2:
{ // calcul des coordonnées locales
theta.Change_dim(2);
theta(1) = OA * giH_(1); theta(2) = OA * giH_(2);
// les fonctions d'interpolation
phi_.Change_taille(3);
phi_(1)=1.-theta(1)-theta(2);
phi_(2)=theta(1);
phi_(3)=theta(2);
}
break;
case 3:
{ // calcul des coordonnées locales
theta.Change_dim(3);
theta(1) = OA * giH_(1); theta(2) = OA * giH_(2);theta(3) = OA * giH_(3);
// les fonctions d'interpolation
phi_.Change_taille(4);
phi_(1)=1.-theta(1)-theta(2)-theta(3);
phi_(2)=theta(1);
phi_(3)=theta(2);
phi_(4)=theta(3);
}
break;
default:
cout << "\n erreur !!! cas dime = " << dime << " non prevu "
<< "\n ElemGeomC0::Bases_naturel_duales(...";
Sortie(1);
};
};
// -- 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 ElemGeomC0::Bases_naturel_duales(const Tableau & indirec
,Tableau & gi_B,Tableau & gi_H
,int cas) const
{// on récupère la dimension de la métrique
int dime = indirec.Taille();
// traitement suivant la taille
switch (dime)
{ case 2: // on interpole linéairement entre deux points définissant une droite
{// calcul de la base naturelle
gi_B.Change_taille(1);
// gi_B(1)= ptInteg(indirec(2))-ptInteg(indirec(1)); // modif 29 oct 2020
gi_B(1)= 0.5*(ptInteg(indirec(2))-ptInteg(indirec(1)));
// calcul de la métrique
double g11BB=gi_B(1)*gi_B(1);
// inverse de la métrique
double det = g11BB;
if (Dabs(det) <= ConstMath::trespetit)
{ cout << "\nAttention : le determinant du systeme est nulle !"
<< " il n'est pas possible de calculer la base duale ";
cout << "\n ElemGeomC0::Bases_naturelles_duales(..";
return false;
};
// double unsurdet = 1./det;
double g11HH = 1./ g11BB;
// base duale
gi_H.Change_taille(1);
gi_H(1)= g11HH * gi_B(1);
}
break;
case 3: // on interpole linéairement entre 3 points définissant un triangle
{// calcul de la base naturelle
gi_B.Change_taille(2);
gi_B(1)= ptInteg(indirec(2))-ptInteg(indirec(1));
gi_B(2)= ptInteg(indirec(3))-ptInteg(indirec(1));
// calcul de la métrique
double g11BB=gi_B(1)*gi_B(1);
double g12BB=gi_B(1)*gi_B(2);
double g22BB=gi_B(2)*gi_B(2);
// inverse de la métrique
double det = g11BB*g22BB - g12BB*g12BB;
if (Dabs(det) <= ConstMath::trespetit)
{ cout << "\nAttention : le determinant du systeme est nulle !"
<< " il n'est pas possible de calculer la base duale ";
cout << "\n ElemGeomC0::Bases_naturelles_duales(..";
return false;
};
double unsurdet = 1./det;
double g11HH = unsurdet * g22BB;
double g22HH = unsurdet * g11BB;
double g12HH = - unsurdet * g12BB;
// base duale
gi_H.Change_taille(2);
gi_H(1)= g11HH * gi_B(1) + g12HH * gi_B(2);
gi_H(2)= g12HH * gi_B(1) + g22HH * gi_B(2);
}
break;
case 4:
if (cas == 1 )// on interpole linéairement entre 4 points définissant un tétraèdre
{// calcul de la base naturelle
gi_B.Change_taille(3);
gi_B(1)= ptInteg(indirec(2))-ptInteg(indirec(1));
gi_B(2)= ptInteg(indirec(3))-ptInteg(indirec(1));
gi_B(3)= ptInteg(indirec(4))-ptInteg(indirec(1));
// calcul de la métrique
Mat_pleine gijBB(3,3);
for (int i=1;i<=3;i++)
for (int j=1;j<=i;j++)
gijBB(i,j)=gijBB(j,i)=gi_B(i)*gi_B(j);
// inverse de la métrique
Mat_pleine gijHH(3,3); // def
{try
{gijHH = gijBB.Inverse();
}
catch (ErrResolve_system_lineaire excep)
// cas d'une erreur pendant l'inversion
{ cout << "\nAttention : l'extrapolation via un tetraedre lineaire !"
<< " n'est pas possible ";
cout << "\n ElemGeomC0::Bases_naturelles_duales(..";
return false;
};
}
// base duale
gi_H.Change_taille(3);
gi_H(1)= gijHH(1,1) * gi_B(1) + gijHH(1,2) * gi_B(2) + gijHH(1,3) * gi_B(3);
gi_H(2)= gijHH(2,1) * gi_B(1) + gijHH(2,2) * gi_B(2) + gijHH(2,3) * gi_B(3);
gi_H(3)= gijHH(3,1) * gi_B(1) + gijHH(3,2) * gi_B(2) + gijHH(3,3) * gi_B(3);
}
else // on interpole bi-linéairement en métrique 2D
{// calcul de la base naturelle
gi_B.Change_taille(2);
gi_B(1)= ptInteg(indirec(2))-ptInteg(indirec(1));
gi_B(2)= ptInteg(indirec(3))-ptInteg(indirec(1));
// calcul de la métrique
double g11BB=gi_B(1)*gi_B(1);
double g12BB=gi_B(1)*gi_B(2);
double g22BB=gi_B(2)*gi_B(2);
// inverse de la métrique
double det = g11BB*g22BB - g12BB*g12BB;
if (Dabs(det) <= ConstMath::trespetit)
{ cout << "\nAttention : le determinant du systeme est nulle !"
<< " il n'est pas possible de calculer la base duale ";
cout << "\n ElemGeomC0::Bases_naturelles_duales(..";
return false;
};
double unsurdet = 1./det;
double g11HH = unsurdet * g22BB;
double g22HH = unsurdet * g11BB;
double g12HH = - unsurdet * g12BB;
// base duale
gi_H.Change_taille(2);
gi_H(1)= g11HH * gi_B(1) + g12HH * gi_B(2);
gi_H(2)= g12HH * gi_B(1) + g22HH * gi_B(2);
}
break;
default:
cout << "\n erreur !!! cas dime = " << dime << " non prevu "
<< "\n ElemGeomC0::Bases_naturel_duales(...";
Sortie(1);
};
// retour ok
return true;
};