// 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 "DeformationPP.h"
# include "Met_PiPoCo.h"
#include "ConstMath.h"
#include "ParaGlob.h"
#include "Util.h"
// il faut regarder et mettre à jour le calcul des différentes déformations ********
// il faut regarder dans déformation, en particulier il faut mettre à jour, ***********
// si nécessaire différents cas de déformation **********************
// il faut également mettre à jour l'influence de l'indicateur premier_calcul ????????????
// dans le cas du premier passage, pour l'instant pas fait !!!!!!!!!!!!!!!!!!!!!!!!!!
// -------------------------- variables static -------------------
int DeformationPP::numInteg_ep = 0;
int DeformationPP::nbtotalep = 0;
int DeformationPP::nbtotalaxpl = 0;
Vecteur DeformationPP::theta_z(1);
// on utilise deux types d'interpolation
Tableau <Mat_pleine const *> DeformationPP::taDphi(2);
Tableau <Vecteur const *> DeformationPP::taPhi(2);
// ---------- constructeur----------------------------------------
DeformationPP::DeformationPP () // constructeur ne devant pas etre utilise
{
#ifdef MISE_AU_POINT
{ cout << "\nErreur : le constructeur par defaut ne doit pa etre utilise !\n";
cout << "DeformationPP::DeformationPP () \n";
Sortie(1);
};
#endif
};
// constructeur normal dans le cas d'un ou de plusieurs pt d'integration
DeformationPP::DeformationPP (Met_abstraite & a,Tableau<Noeud *>& tabnoeud
,Tableau <Mat_pleine> const & tDphH,Tableau <Vecteur> const & tPhH
,Tableau <Mat_pleine> const & tDphS,Tableau <Vecteur> const & tPhS ):
Deformation (a,tabnoeud,tDphS,tPhS )
{ type_deformation = DEFORMATION_POUTRE_PLAQUE_STANDART;
tabDphiH = &(tDphH);
tabPhiH = &(tPhH);
// a priori pas de numero d'integration affectes
numInteg = 0;
numInteg_ep = 0;
};
// constructeur de copie
DeformationPP::DeformationPP (const DeformationPP& a) :
Deformation(a)
{ numInteg_ep = a.numInteg_ep;
tabDphiH = a.tabDphiH;
tabPhiH = a.tabPhiH;
};
DeformationPP::~DeformationPP ()
{ };
// ============ METHODES PUBLIQUES : ==================
// Surcharge de l'operateur = : realise l'affectation
// fonction virtuelle
Deformation& DeformationPP::operator= (const Deformation& )
{ // normalement ne devrait pas être utilisé
cout << "\n erreur , l operateur d affectation a utiliser doit etre celui explicite "
<< " de la classe DeformationPP \n"
<< " Deformation& DeformationPP::operator= (const Deformation& def) \n";
Sortie(1);
return (*this);
};
// fonction normale
DeformationPP& DeformationPP::operator= (const DeformationPP& def)
{ tabDphiH = def.tabDphiH;
tabPhiH = def.tabPhiH;
(*this) = Deformation::operator=(def);
return (*this);
};
// change les numeros d'integration de surface et d'epaisseur courant
void DeformationPP::ChangeNumIntegSH
(int nisurf, int niepaiss)
{ //PiPoCo & elem = *((PiPoCo*) element);
numInteg_ep = niepaiss; numInteg = nisurf;
// on fait les changements dus aux nouveaux pts integ
AppliquePtInteg();
} ;
// calcul explicit à t :tous les parametres sont de resultats
const Met_abstraite::Expli& DeformationPP::Cal_explicit_t
( const Tableau <double>& def_equi_t,TenseurBB & epsBB_t,Tableau <TenseurBB *> & d_epsBB
,Tableau <double>& def_equi,TenseurBB& DepsBB,TenseurBB& DeltaEpsBB,bool premier_calcul)
{ bool gradV_instantane = false; // ************ pour l'instant figé
// appel de la metrique
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
const Met_abstraite::Expli& ex = met_pipoco->Cal_explicit_t
(*tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul);
/// voir Deformation car ça a changé !!
/* // epsBB_t est dimensionne dans l'element
epsBB_t = 0.5 * ( *(ex.gijBB_t) - *(ex.gijBB_0));
// ici on considère que le delta_epsBB=epsBB_t ! c-a-d le delta de 0 à t
DeltaEpsBB = epsBB_t;
// calcul simplifié éventuelle de la vitesse de déformation si les ddl vitesses ne sont pas présent
if (pas_de_gradV)
{ // on choisit une vitesse de déformation pour l'instant identique à deltaeps/t
const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
// dans le cas où l'incrément de temps est nul on garde la vitesse actuelle
double deltat=vartemps.IncreTempsCourant();
if (deltat >= ConstMath::trespetit)
DepsBB = DeltaEpsBB/deltat;
}
else // cas où le gradient a été calculé
DepsBB = 0.5 * ((*ex.gradVBB_t) + ex.gradVBB_t->Transpose());
// variation de la déformation / au ddl
int d_epsBBTaille = d_epsBB.Taille();
for (int i=1; i<= d_epsBBTaille; i++)
*(d_epsBB(i)) = 0.5 * (*((*(ex.d_gijBB_t))(i)));*/
return ex;
};
// calcul explicit à tdt :tous les parametres sont de resultats
const Met_abstraite::Expli_t_tdt& DeformationPP::Cal_explicit_tdt
( const Tableau <double>& def_equi_t,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB
,Tableau <double>& def_equi,TenseurBB& DepsBB,TenseurBB& delta_epsBB_tdt,bool premier_calcul)
{ bool gradV_instantane = false; // ************ pour l'instant figé
// appel de la metrique
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
const Met_abstraite::Expli_t_tdt& ex= met_pipoco->Cal_explicit_tdt
(*tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul);
/// voir Deformation car ça a changé !!
/* // epsBB_tdt est dimensionne auparavant
epsBB_tdt = 0.5 * ( *(ex.gijBB_t) - *(ex.gijBB_0));
// epsBB_tdt est dimensionne auparavant
delta_epsBB_tdt = 0.5 * ( *(ex.gijBB_tdt) - *(ex.gijBB_t));
if (pas_de_gradV)
{ // on choisit une vitesse de déformation pour l'instant identique à deltaeps/t
const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
// dans le cas où l'incrément de temps est nul on garde la vitesse actuelle
double deltat=vartemps.IncreTempsCourant();
if (deltat >= ConstMath::trespetit)
DepsBB = delta_epsBB_tdt/deltat;
}
else // cas où le gradient a été calculé
DepsBB = 0.5 * ((*ex.gradVBB_tdt) + ex.gradVBB_tdt->Transpose());
// variation de la déformation / au ddl
int d_epsBBTaille = d_epsBB.Taille();
for (int i=1; i<= d_epsBBTaille; i++)
*(d_epsBB(i)) = 0.5 * (*((*(ex.d_gijBB_tdt))(i)));*/
return ex;
};
// cas implicite : tous les parametres sont de resultats
const Met_abstraite::Impli& DeformationPP::Cal_implicit
( const Tableau <double>& def_equi_t,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB_tdt
,Tableau <double>& def_equi,Tableau2 <TenseurBB *>& d2_epsBB_tdt
,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul)
{ bool gradV_instantane = false; // ************ pour l'instant figé
// recup d'un pointeur sur la metrique sfe1
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
// toutes les variables de passage a metrique apres l'appel
// pointeront sur des variables deja dimensionnees
// pour les Tableau <> il y a dimensionnement auto a l'affectation
const Met_abstraite::Impli& ex =met_pipoco->Cal_implicit ( *tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul);
/// voir Deformation car ça a changé !!
/* epsBB_tdt = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_0));
delta_epsBB = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_t));
if (pas_de_gradV)
{ // on choisit une vitesse de déformation pour l'instant identique à deltaeps/t
const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
// dans le cas où l'incrément de temps est nul on garde la vitesse actuelle
double deltat=vartemps.IncreTempsCourant();
if (deltat >= ConstMath::trespetit)
DepsBB = delta_epsBB/deltat;
}
else // cas où le gradient a été calculé
DepsBB = 0.5 * ((*ex.gradVBB_tdt) + ex.gradVBB_tdt->Transpose());
// variation de la déformation / au ddl
int d_epsBB_tdtTaille = d_epsBB_tdt.Taille();
for (int i=1; i<= d_epsBB_tdtTaille; i++)
*(d_epsBB_tdt(i)) = 0.5 * (*((*(ex.d_gijBB_tdt))(i)));
// calcul éventuelle de la dérivée seconde de la déformation
if (pa.Var_D())
{// recup des variations secondes de la déformation
int d2_epsBB_tdtTaille1 = d2_epsBB_tdt.Taille1();
for (int i=1; i<= d2_epsBB_tdtTaille1; i++)
for (int j=1; j<= i; j++) // symétrie du tableau et du tenseur
{ *(d2_epsBB_tdt(i,j)) = 0.5 * (*((*(ex.d2_gijBB_tdt))(i,j))) ;
*(d2_epsBB_tdt(j,i)) = *(d2_epsBB_tdt(i,j)) ;
}
} */
return ex;
};
// ---------------- calcul des variables primaires autre que pour la mécanique --------
// ------------ donc par de retour relatif aux déformations
// calcul explicit à t :tous les parametres sont de resultats
const Met_abstraite::Expli& DeformationPP::Cal_explicit_t(bool premier_calcul )
{ bool gradV_instantane = false; // ************ pour l'instant figé
// appel de la metrique
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
// appel de la metrique
const Met_abstraite::Expli& ex = met_pipoco->Cal_explicit_t(*tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul);
return ex;
};
// calcul explicit à tdt :tous les parametres sont de resultats
const Met_abstraite::Expli_t_tdt& DeformationPP::Cal_explicit_tdt(bool premier_calcul )
{ bool gradV_instantane = false; // ************ pour l'instant figé
// appel de la metrique
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
// appel de la metrique
const Met_abstraite::Expli_t_tdt& ex = met_pipoco->Cal_explicit_tdt(*tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul);
return ex;
};
// cas implicite : tous les parametres sont de resultats
const Met_abstraite::Impli& DeformationPP::Cal_implicit( bool premier_calcul)
{ bool gradV_instantane = false; // ************ pour l'instant figé
// recup d'un pointeur sur la metrique sfe1
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
// toutes les variables de passage a metrique apres l'appel
// pointeront sur des variables deja dimensionnees
// pour les Tableau <> il y a dimensionnement auto a l'affectation
// appel de la metrique
const Met_abstraite::Impli& ex =met_pipoco->Cal_implicit ( *tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul);
return ex;
};
// ========== remontee aux informations =========================
// cas sortie d'un calcul implicit
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * I^a
// tout ce passe comme si I^a est la nouvelle base vers laquelle on veut évoluer
const Met_abstraite::InfoImp DeformationPP::RemontImp(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
{ // recup d'un pointeur sur la metrique sfe1
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
Met_abstraite::InfoImp ex = met_pipoco->Cal_InfoImp
( *tabnoeud,taDphi,taPhi,nbNoeud);
// determination des matrices de transformation de base
const BaseB & giB0 = *(ex.giB_0);
const BaseB & giB = *(ex.giB_tdt);
const BaseH & giH0 = *(ex.giH_0);
const BaseH & giH = *(ex.giH_tdt);
// on différencie le cas des poutres et celui des plaques
int dim = giB0.NbVecteur();
if (dim == 1) // cas de la dimension 1 c-a-dire poutre
Deformation::BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin); // générale
else
BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin); // particulier plaques
return ex;
};
// cas sortie d'un calcul implicit
// sans les matrices de passage
const Met_abstraite::InfoImp DeformationPP::RemontImp()
{ // recup d'un pointeur sur la metrique sfe1
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
Met_abstraite::InfoImp ex = met_pipoco->Cal_InfoImp
( *tabnoeud,taDphi,taPhi,nbNoeud);
// determination des matrices de transformation de base
const BaseB & giB0 = *(ex.giB_0);
const BaseB & giB = *(ex.giB_tdt);
const BaseH & giH0 = *(ex.giH_0);
const BaseH & giH = *(ex.giH_tdt);
return ex;
};
// cas sortie d'un calcul explicit à t
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * I^a
// tout ce passe comme si I^a est la nouvelle base vers laquelle on veut évoluer
const Met_abstraite::InfoExp_t DeformationPP::RemontExp_t(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
{ // recup d'un pointeur sur la metrique sfe1
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
Met_abstraite::InfoExp_t ex = met_pipoco->Cal_InfoExp_t
( *tabnoeud,taDphi,taPhi,nbNoeud);
// determination des matrices de transformation de base
const BaseB & giB0 = *(ex.giB_0);
const BaseB & giB = *(ex.giB_t);
const BaseH & giH0 = *(ex.giH_0);
const BaseH & giH = *(ex.giH_t);
BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);
return ex;
};
// cas sortie d'un calcul explicit à t
// sans les matrices de passage
const Met_abstraite::InfoExp_t DeformationPP::RemontExp_t()
{ // recup d'un pointeur sur la metrique sfe1
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
Met_abstraite::InfoExp_t ex = met_pipoco->Cal_InfoExp_t
( *tabnoeud,taDphi,taPhi,nbNoeud);
// determination des matrices de transformation de base
const BaseB & giB0 = *(ex.giB_0);
const BaseB & giB = *(ex.giB_t);
const BaseH & giH0 = *(ex.giH_0);
const BaseH & giH = *(ex.giH_t);
return ex;
};
// cas sortie d'un calcul explicit à tdt
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * I^a
// tout ce passe comme si I^a est la nouvelle base vers laquelle on veut évoluer
const Met_abstraite::InfoExp_tdt DeformationPP::RemontExp_tdt(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
{ // recup d'un pointeur sur la metrique sfe1
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
Met_abstraite::InfoExp_tdt ex = met_pipoco->Cal_InfoExp_tdt
( *tabnoeud,taDphi,taPhi,nbNoeud);
// determination des matrices de transformation de base
const BaseB & giB0 = *(ex.giB_0);
const BaseB & giB = *(ex.giB_tdt);
const BaseH & giH0 = *(ex.giH_0);
const BaseH & giH = *(ex.giH_tdt);
BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);
return ex;
};
// cas sortie d'un calcul explicit à tdt
// sans les matrices de passage
const Met_abstraite::InfoExp_tdt DeformationPP::RemontExp_tdt()
{ // recup d'un pointeur sur la metrique sfe1
Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
Met_abstraite::InfoExp_tdt ex = met_pipoco->Cal_InfoExp_tdt
( *tabnoeud,taDphi,taPhi,nbNoeud);
// determination des matrices de transformation de base
const BaseB & giB0 = *(ex.giB_0);
const BaseB & giB = *(ex.giB_tdt);
const BaseH & giH0 = *(ex.giH_0);
const BaseH & giH = *(ex.giH_tdt);
return ex;
};
// gestion du parcours de tous les points d'integration
// bien voir que numInteg est relatif aux pt d'intégration de l'axe
void DeformationPP::PremierPtInteg()
{ nbtotalaxpl = tabPhi->Taille();
nbtotalep = tabPhiH->Taille();
numInteg = 1 ; numInteg_ep = 1;
// on fait les changements dus aux nouveaux pts integ
AppliquePtInteg();
};
bool DeformationPP::DernierPtInteg()
{ if (numInteg <= nbtotalaxpl)
return true;
else
return false;
};
void DeformationPP::NevezPtInteg()
{ if (numInteg_ep < nbtotalep)
numInteg_ep++;
else
{ numInteg++;
numInteg_ep = 1;
}
// on fait les changements dus aux nouveaux pts integ
// dans le cas ou ils sont valides
if (numInteg <= nbtotalaxpl)
AppliquePtInteg();
};
// ------------------------ METHODES PROTEGEES : -------------------------------
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * I^a
// tout ce passe comme si I^a est la nouvelle base vers laquelle on veut évoluer
void DeformationPP::BasePassage
(bool absolue,const BaseB & giB0,const BaseB & giB,const BaseH & giH0,const BaseH & giH,
Mat_pleine& Aa0,Mat_pleine& Aafin)
{
// determination des matrices de transformation de base
// l'objectif est de determiner un repere pertinant dans le cas de plaque , coque .
// Dans le cas d'une dimension globale 3D avec des plaques ou coque ->
// choix : un repere qui appartiend a la facette, obtenu apres projection
// du repere global
//------ cas de la config initiale, on regarde si la projection de I1 n'est pas trop petite
// def de la normale a la facette
Coordonnee N = (Util::ProdVec_coorBN(giB0(1),giB0(2))).Normer();
Coordonnee ve(0.,-N(3),-N(2)); // = ProdVec(N,I1)
double norve = ve.Norme();
Tableau <Coordonnee> vi(2); // les vecteurs plans orthonormes de la facette
if (norve >= ConstMath::petit)
{ vi(2) = ve.Normer();
vi(1) = Util::ProdVec_coor(vi(2),N);
}
else
{ ve.Change_Coordonnee(3,N(3),0.,-N(1)); // = ProdVec(I2,N)
vi(1) = ve.Normer();
vi(2) = Util::ProdVec_coor(N,vi(1));
}
for (int al=1 ;al<=2;al++)
{Aa0(1,al) = giH0.Coordo(1) * vi(al) ;
Aa0(2,al) = giH0.Coordo(2) * vi(al) ;
}
//------ cas de la config finale,
N = (Util::ProdVec_coorBN(giB(1),giB(2))).Normer();
ve.Change_Coordonnee(3,0.,N(3),-N(2)); // = ProdVec(N,I1)
norve = ve.Norme();
Tableau <Coordonnee> Aa(2);
if (norve >= ConstMath::petit)
{ Aa(2) = ve.Normer();
Aa(1) = Util::ProdVec_coor(Aa(2),N);
}
else
{ ve.Change_Coordonnee(3,N(3),0.,-N(1)); // = ProdVec(I2,N)
Aa(1) = ve.Normer();
Aa(2) = Util::ProdVec_coor(N,Aa(1));
}
for (int be=1;be<=2;be++)
{ Aafin(1,be) = giH.Coordo(1) * Aa(be);
Aafin(2,be) = giH.Coordo(2) * Aa(be);
}
return;
};
// applique les conséquences d'un changement de numéro de point d'intégration
// par exemple effectué via NevezPtInteg() ou ChangeNumIntegSH
void DeformationPP::AppliquePtInteg()
{ //PiPoCo & elem = *((PiPoCo*) element);
// on est obligé de préciser le caractère réel du tableau ?
//------- pour l'instant c'est inopérationnelle: les 3 lignes qui suivent donne une erreur à la compilation
// taPhi(1) = & ((Vecteur const)((*tabPhi)(numInteg)));
// taDphi(1) = &((Mat_pleine const)((*tabDphi)(numInteg)));
// taDphi(2) = &((Mat_pleine const)((*tabDphiH)(numInteg_ep)));
// calcul de la cote en epaisseur
// les trois lignes suivantes ont été mise en commentaire car on ne peut plus
// accèder à l'élément directement, a moins par exemple de le passer en
// paramètre ou de trouver une autre technique pour avoir les infos a doc
cout << "\n pour l'instant la méthode n'est pas disponible cf le listing";
cout << " void DeformationPP::AppliquePtInteg()";
Sortie(1);
// double epaisseur = elem.H();
// theta_z(1) = 0.5 * epaisseur * (elem.KSIepais(numInteg_ep));
// taPhi(2) = & theta_z;
} ;