rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/Elements/Mecanique/Deformation_gene/PiPoCo.cc
Gérard Rio 0460624e94 mise à jour 7.010
2023-05-03 17:23:49 +02:00

345 lines
19 KiB
C++
Raw Permalink Blame History

// 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 "PiPoCo.h"
#include <iomanip>
#include "ConstMath.h"
#include "Loi_Umat.h"
#include "ExceptionsElemMeca.h"
#include "ExceptionsLoiComp.h"
// =========================== variables communes ===============================
// a priori les pointeurs de fonction pointent sur rien au debut
TenseurHH* PiPoCo::sig_bulk_pourPiPoCo_HH = NULL; // variable de travail pour le bulk
// =========================== fin variables communes ===========================
// =========================== methodes publiques ===============================
// calcul si un point est a l'interieur de l'element ou non
// il faut que M est la dimension globale
// les trois fonctions sont pour l'etude a t=0, t et tdt
// retour : =0 le point est externe, =1 le point est interne ,
// = 2 le point est sur la frontière à la précision près
// coor_locales : s'il est différent de NULL, est affecté des coordonnées locales calculées,
// uniquement précises si le point est interne
int PiPoCo::Interne_0(const Coordonnee& M,Coordonnee* )
{ cout << "\n fonction non implante, bool PiPoCo::Interne_0(Coordonnee& M)";
Sortie(1);
return Interne(TEMPS_0,M);
};
int PiPoCo::Interne_t(const Coordonnee& M,Coordonnee* )
{ cout << "\n fonction non implante, bool PiPoCo::Interne_t(Coordonnee& M)";
Sortie(1);
return Interne(TEMPS_t,M);
};
int PiPoCo::Interne_tdt(const Coordonnee& M,Coordonnee* )
{ cout << "\n fonction non implante, bool PiPoCo::Interne_tdt(Coordonnee& M)";
Sortie(1);
return Interne(TEMPS_tdt,M);
};
// METHODES VIRTUELLES:
// --------- calculs utils dans le cadre de la recherche du flambement linéaire
// dans un premier temps uniquement virtuelles, ensuite se sera virtuelle pure pour éviter
// les oubli de définition ----> AMODIFIER !!!!!!!!
// Calcul de la matrice géométrique et initiale
ElemMeca::MatGeomInit PiPoCo::MatricesGeometrique_Et_Initiale (const ParaAlgoControle & )
{ cout << "\n attention le calcul de la matrice g\'eom\'etrique et de la matrice "
<< " initiale ne sont pas implant\'ees \n" ;
this->Affiche(1);
cout << "\nvoid PiPoCo::MatricesGeometrique_Et_Initiale (Mat_pleine &,Mat_pleine &)"
<< endl;
Sortie(1);
return MatGeomInit(NULL,NULL); // pour supprimer le warning
};
// =========================== methodes protegees ===============================
// Calcul du residu local et de la raideur locale,
// pour le schema implicite
void PiPoCo::Cal_implicitPiPoCo (DdlElement & tab_ddl,
Tableau <TenseurBB *> & d_epsBB,Tableau < Tableau2 <TenseurBB *> > d2_epsBB_,
Tableau <TenseurHH *> & d_sigHH,int nbintS,Vecteur& poidsS,int nbintH
,Vecteur& poidsH,const ParaAlgoControle & pa,bool cald_Dvirtuelle)
{ int nbddl = tab_ddl.NbDdl();
double jacobien; // def du jacobien
Vecteur d_jacobien_tdt;
energie_totale.Inita(0.); // init de l'énergie totale sur l'élément
TenseurHH* sigHH_t_1 = (*lesPtIntegMecaInterne)(1).SigHH_t(); // simplification d'écriture
if (ElemMeca::Bulk_visco_actif())
{ if (sig_bulk_pourPiPoCo_HH == NULL) {sig_bulk_pourPiPoCo_HH = NevezTenseurHH(*sigHH_t_1);}
else if (sig_bulk_pourPiPoCo_HH->Dimension() != sigHH_t_1->Dimension())
{delete sig_bulk_pourPiPoCo_HH; sig_bulk_pourPiPoCo_HH = NevezTenseurHH(*sigHH_t_1);};
};
int indice = 1;
for (int nisur =1; nisur<= nbintS; nisur++) // boucle sur les pt d'integ de surface
for (int niepais =1; niepais<= nbintH; niepais++, indice++) // pt d'integ d'epaisseur
{ ((DeformationPP*) def)->ChangeNumIntegSH(nisur,niepais);
((DeformationPP*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(indice);
TenseurHH & sigHH = *(ptIntegMeca.SigHH());
TenseurBB & DepsBB_tdt = *(ptIntegMeca.DepsBB());
EnergieMeca& energ = tab_energ(indice);
EnergieMeca& energ_t = tab_energ_t(indice);
double poids = poidsS(nisur) * poidsH(niepais) ;
CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques
if (loiTP != NULL) {sTP = tabSaveTP(indice);}; // au pt d'integ si elles existes
if ((loiComp->Id_comport()==LOI_VIA_UMAT) // cas de l'utilisation d'une loi umat
|| (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),indice);
// -- on gère les exceptions éventuelles en mettant le bloc sous surveillance
const Met_abstraite::Impli* ex_pointe = NULL; // c'est pour construire ex
try // ce qui permet de donner les numéros d'éléments et de pti
{ ex_pointe = &loiComp->Cal_implicit(tabSaveDon(indice), *def,tab_ddl,ptIntegMeca
,d_epsBB,jacobien,d_jacobien_tdt,d_sigHH
,pa,sTP,loiTP,dilatation,energ,energ_t
,premier_calcul_meca_impli_expli);
}
catch (ErrNonConvergence_loiDeComportement excep)
// cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
// d'une loi de comportement incrémentale
{ cout << "\n erreur de loi de comportement element= " << this->Num_elt()
<< " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
throw excep;
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... )
{ cout << "\n erreur inconnue de loi de comportement, element= " << this->Num_elt()
<< " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
throw (Err_inconnue_ElemMeca());
};
const Met_abstraite::Impli& ex=(*ex_pointe);
// on calcul éventuellement la dérivée de la vitesse de déformation virtuelle
Tableau2 <TenseurBB *>& d2_epsBB = d2_epsBB_(indice); // pour simplifier
if (cald_Dvirtuelle)
def->Cal_var_def_virtuelle(false,d2_epsBB);
// examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
if (Bulk_visco_actif())
ModifContrainteBulk(TEMPS_tdt,*ex.gijHH_tdt,sigHH,DepsBB_tdt,(*sig_bulk_pourPiPoCo_HH) );
if (jacobien <= 0.) // vérif du jacobien
{ cout << "\n ********** attention on a trouve un jacobien negatif !! ********* " << endl;
cout << "\n on annulle sa contribution. Element nb: " << Num_elt() << " nb d'integ : " << indice;
}
else
{ // cas avec toutes les variations
{energie_totale += energ * (poids * jacobien);
for (int j =1; j<= nbddl; j++) // 1ere boucle sur les ddl
{(*residu)(j) -= ((*(d_epsBB(j))) && sigHH ) * poids * jacobien;
for (int i =1; i<= nbddl; i++) // 2ere boucle sur les ddl
(*raideur)(j,i) += ( ((*d_sigHH(i)) && (*(d_epsBB(j)))) * jacobien
+ (sigHH && (*(d2_epsBB(j,i)))) * jacobien
+ (sigHH && (*(d_epsBB(j)))) * d_jacobien_tdt(i)
) * poids;
}
}
}//-- fin des différents cas
if (Bulk_visco_actif()) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
{ sigHH -= (*sig_bulk_pourPiPoCo_HH); } // ceci pour la sauvegarde ou autre utilisation
}; // fin boucle sur les points d'intégration d'épaisseur
// liberation des tenseurs intermediaires
LibereTenseur();
};
// Calcul uniquement du residu local a l'instant t ou pas suivant la variable atdt
// pour le schema explicit par exemple
void PiPoCo::Cal_explicitPiPoCo (DdlElement & tab_ddl,Tableau <TenseurBB *> & d_epsBB,int nbintS,
Vecteur& poidsS,int nbintH,Vecteur& poidsH,const ParaAlgoControle & ,bool atdt)
{ double jacobien; // def du jacobien
int nbddl = tab_ddl.NbDdl();
int indice = 1;
energie_totale.Inita(0.); // init de l'énergie totale sur l'élément
TenseurHH* sigHH_t_1 = (*lesPtIntegMecaInterne)(1).SigHH_t(); // simplification d'écriture
if (ElemMeca::Bulk_visco_actif())
{ if (sig_bulk_pourPiPoCo_HH == NULL) {sig_bulk_pourPiPoCo_HH = NevezTenseurHH(*sigHH_t_1);}
else if (sig_bulk_pourPiPoCo_HH->Dimension() != sigHH_t_1->Dimension())
{delete sig_bulk_pourPiPoCo_HH; sig_bulk_pourPiPoCo_HH = NevezTenseurHH(*sigHH_t_1);};
};
for (int nisur =1; nisur<= nbintS; nisur++) // boucle sur les pt d'integ de surface
for (int niepais =1; niepais<= nbintH; niepais++, indice++) // pt d'integ d'epaisseur
{ ((DeformationPP*) def)->ChangeNumIntegSH(nisur,niepais);
((DeformationPP*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(indice);
TenseurHH & sigHH = *(ptIntegMeca.SigHH());
TenseurBB & DepsBB_ = *(ptIntegMeca.DepsBB());
EnergieMeca& energ = tab_energ(indice);
EnergieMeca& energ_t = tab_energ_t(indice);
double poids = poidsS(nisur) * poidsH(niepais) ;
CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques
if (loiTP != NULL) {sTP = tabSaveTP(indice);}; // au pt d'integ si elles existes
if ((loiComp->Id_comport()==LOI_VIA_UMAT) // cas de l'utilisation d'une loi umat
|| (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),indice);
// -- on gère les exceptions éventuelles en mettant le bloc sous surveillance
try // ce qui permet de donner les numéros d'éléments et de pti
{if (atdt)
{ const Met_abstraite::Expli_t_tdt& ex = loiComp->Cal_explicit_tdt
(tabSaveDon(indice), *def,tab_ddl,ptIntegMeca,d_epsBB,jacobien
,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_meca_impli_expli);
// examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
if (Bulk_visco_actif())
ModifContrainteBulk(TEMPS_tdt,*ex.gijHH_tdt,sigHH,DepsBB_ ,(*sig_bulk_pourPiPoCo_HH));
}
else
{ const Met_abstraite::Expli& ex = loiComp->Cal_explicit_t
(tabSaveDon(indice), *def,tab_ddl,ptIntegMeca,d_epsBB,jacobien
,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_meca_impli_expli);
// examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
if (Bulk_visco_actif())
ModifContrainteBulk(TEMPS_t,*ex.gijHH_t,sigHH,DepsBB_ ,(*sig_bulk_pourPiPoCo_HH));
};
}
catch (ErrNonConvergence_loiDeComportement excep)
// cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
// d'une loi de comportement incrémentale
{ cout << "\n erreur de loi de comportement element= " << this->Num_elt()
<< " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
throw excep;
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... )
{ cout << "\n erreur inconnue de loi de comportement, element= " << this->Num_elt()
<< " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
throw (Err_inconnue_ElemMeca());
};
if (jacobien <= 0.) // vérif du jacobien
{ cout << "\n ********** attention on a trouve un jacobien negatif !! ********* " << endl;
cout << "\n on annulle sa contribution. Element nb: " << Num_elt() << " nb d'integ : " << indice;
}
else
{ energie_totale += energ * (poids * jacobien);
for (int j =1; j<= nbddl; j++)
(*residu)(j) -= (sigHH && (*(d_epsBB(j)))) * jacobien * poids;
}
if (Bulk_visco_actif()) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
{ sigHH -= (*sig_bulk_pourPiPoCo_HH); } // ceci pour la sauvegarde ou autre utilisation
}; // fin boucle sur les points d'intégration
// liberation des tenseurs intermediaires
LibereTenseur();
};
// Calcul de la matrice géométrique et de la matrice initiale
// cette fonction est éventuellement appelée par les classes dérivées
// ddl represente les degres de liberte specifiques a l'element
// epsBB = deformation, sigHH = contrainte, d_epsbb = variation des def
// nbint = nb maxi de pt d'integration , poids = poids d'integration
// S pour surface et H pour épaisseur -> nbint et poids
void PiPoCo::Cal_matGeom_InitPiPoCo (Mat_pleine & matGeom, Mat_pleine & matInit,
DdlElement & tab_ddl,Tableau <TenseurBB *>& d_epsBB, Tableau < Tableau2 <TenseurBB *> > d2_epsBB_,
Tableau <TenseurHH *> & d_sigHH,int nbintS
,Vecteur& poidsS,int nbintH,Vecteur& poidsH,const ParaAlgoControle & pa,bool cald_Dvirtuelle)
{ // le début du programme est semblable au calcul implicit
// seule la constitution des matrices finales est différente
int nbddl = tab_ddl.NbDdl();
double jacobien; // def du jacobien
Vecteur d_jacobien_tdt;
int indice = 1;
for (int nisur =1; nisur<= nbintS; nisur++) // boucle sur les pt d'integ de surface
for (int niepais =1; niepais<= nbintH; niepais++, indice++) // pt d'integ d'epaisseur
{ ((DeformationPP*) def)->ChangeNumIntegSH(nisur,niepais);
((DeformationPP*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(indice);
TenseurHH & sigHH = *(ptIntegMeca.SigHH());
TenseurBB & DepsBB_tdt = *(ptIntegMeca.DepsBB());
EnergieMeca& energ = tab_energ(indice);
EnergieMeca& energ_t = tab_energ_t(indice);
double poids = poidsS(nisur) * poidsH(niepais) ;
CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques
if (loiTP != NULL) {sTP = tabSaveTP(indice);}; // au pt d'integ si elles existes
if ((loiComp->Id_comport()==LOI_VIA_UMAT) // cas de l'utilisation d'une loi umat
|| (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),indice);
// -- on gère les exceptions éventuelles en mettant le bloc sous surveillance
try // ce qui permet de donner les numéros d'éléments et de pti
{loiComp->Cal_flamb_lin(tabSaveDon(indice), *def,tab_ddl,ptIntegMeca,d_epsBB
,jacobien,d_jacobien_tdt, d_sigHH,pa,sTP,loiTP,dilatation
,energ,energ_t,premier_calcul_meca_impli_expli);
}
catch (ErrNonConvergence_loiDeComportement excep)
// cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
// d'une loi de comportement incrémentale
{ cout << "\n erreur de loi de comportement element= " << this->Num_elt()
<< " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
throw excep;
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... )
{ cout << "\n erreur inconnue de loi de comportement, element= " << this->Num_elt()
<< " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
throw (Err_inconnue_ElemMeca());
};
// on calcul éventuellement la dérivée de la vitesse de déformation virtuelle
Tableau2 <TenseurBB *>& d2_epsBB = d2_epsBB_(indice); // pour simplifier
if (cald_Dvirtuelle)
def->Cal_var_def_virtuelle(false,d2_epsBB);
for (int j =1; j<= nbddl; j++) // 1ere boucle sur les ddl
for (int i =1; i<= nbddl; i++) // 2ere boucle sur les ddl
{matInit(j,i) += ((*d_sigHH(i)) && (*(d_epsBB(j)))) * jacobien * poids;
matGeom(j,i) += (sigHH && (*(d2_epsBB(j,i)))) * jacobien * poids;
}
}
// liberation des tenseurs intermediaires
LibereTenseur();
};
Reference in a new issue View git blame Copy permalink