rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/Algo/GalerkinContinu/AlgoDynaExplicite/Algori_relax_dyna.cc

3664 lines
207 KiB
C++
Raw Normal View History

ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// 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.
//
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// 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 "Algori_relax_dyna.h"
#include "MatDiag.h"
using namespace std; //introduces namespace std
#include <cmath>
#include "TypeQuelconqueParticulier.h"
// ---------- pour le calcul de la masse pour la méthode de relaxation dynamique
// masseRelax_1 : pour le premier type de calcul de masse : type_calcul_mass= 1
Ddl_enum_etendu AlgoriRelaxDyna::masseRelax_1; // init à rien, def à l'initialisation
// masseRelax_2 : pour le deuxième type de calcul de masse: type_calcul_mass= 2
TypeQuelconque AlgoriRelaxDyna::masseRelax_2;
// CONSTRUCTEURS :
AlgoriRelaxDyna::AlgoriRelaxDyna () : // par defaut
Algori()
,cL_a_chaque_iteration(0)
,delta_t(),typeCalRelaxation(1)
,lambda(0.605),alph(1.),beta(1.),gamma(1.),theta(1.),casMass_relax(3)
,lambda_initial(0.605)
// pilotage automatique éventuel
,lambda_min(0.605),lambda_max(50.),delta_lambda(0.1)
,pilotage_auto_lambda(false),list_iter_relax()
,niveauLambda_grandeurGlobale(NULL),niveauLambda_ddlEtendu(NULL),niveauLambda_temps(NULL)
// calcul masse
,option_recalcul_mass(1,0)
,nom_fct_nD_option_recalcul_mass(1,""),fct_nD_option_recalcul_mass(1)
,ncycle_calcul(100),type_calcul_mass(2)
,fac_epsilon(1),type_activation_contact(1),choix_mini_masse_nul(0)
,proportion_cinetique(0.1),visqueux_activer(false)
,et_recalcul_masse_a_la_transition(1)
,niveauF_grandeurGlobale(NULL),niveauF_ddlEtendu(NULL),niveauF_temps(NULL)
,Ass1_(NULL),Ass2_(NULL),Ass3_(NULL)
,cas_combi_ddl(),icas(),compteur(0),prepa_avec_remont(false)
,brestart(false),type_incre(OrdreVisu::PREMIER_INCRE)
,compteur_demarrage(-1)
,vglobin(),vglobex(),vglobaal(),vcontact()
,save_X_t()
,X_Bl(),V_Bl(),G_Bl(),forces_vis_num(0)
,li_gene_asso(),t_assemb(),tenuXVG(),mat_masse(NULL),mat_masse_sauve(NULL),mat_C_pt(NULL)
,v_mass(),v_mass1(),matglob(NULL)
{ //vglobaal = &vglobin;
fct_nD_option_recalcul_mass(1)=NULL; // init
// message d'erreur
cout << "\n constructeur par defaut de AlgoriRelaxDyna, ne doit pas etre utilise !!";
Sortie(1);
};
// constructeur en fonction du type de calcul et du sous type
// il y a ici lecture des parametres attaches au type
AlgoriRelaxDyna::AlgoriRelaxDyna (const bool avec_typeDeCal
,const list <EnumSousTypeCalcul>& soustype
,const list <bool>& avec_soustypeDeCal
,UtilLecture& entreePrinc) :
Algori(RELAX_DYNA,avec_typeDeCal,soustype,avec_soustypeDeCal,entreePrinc)
,cL_a_chaque_iteration(0)
,delta_t(0.),typeCalRelaxation(1)
,lambda(0.605),alph(1.),beta(1.),gamma(1.),theta(1.),casMass_relax(3)
,lambda_initial(0.605)
// pilotage automatique éventuel
,lambda_min(0.605),lambda_max(50.),delta_lambda(0.1)
,pilotage_auto_lambda(false),list_iter_relax()
,niveauLambda_grandeurGlobale(NULL),niveauLambda_ddlEtendu(NULL),niveauLambda_temps(NULL)
// calcul masse
,option_recalcul_mass(1,0)
,nom_fct_nD_option_recalcul_mass(1,""),fct_nD_option_recalcul_mass(1)
,ncycle_calcul(100),type_calcul_mass(2)
,fac_epsilon(1),type_activation_contact(1),choix_mini_masse_nul(0)
,proportion_cinetique(0.1),visqueux_activer(false)
,et_recalcul_masse_a_la_transition(1)
,niveauF_grandeurGlobale(NULL),niveauF_ddlEtendu(NULL),niveauF_temps(NULL)
,Ass1_(NULL),Ass2_(NULL),Ass3_(NULL)
,cas_combi_ddl(),icas(),compteur(0),prepa_avec_remont(false)
,brestart(false),type_incre(OrdreVisu::PREMIER_INCRE)
,compteur_demarrage(-1)
,vglobin(),vglobex(),vglobaal(),vcontact()
,save_X_t()
,X_Bl(),V_Bl(),G_Bl(),forces_vis_num(0)
,li_gene_asso(),t_assemb(),tenuXVG(),mat_masse(NULL),mat_masse_sauve(NULL),mat_C_pt(NULL)
,v_mass(),v_mass1(),matglob(NULL)
{ //vglobaal = &vglobin;
fct_nD_option_recalcul_mass(1)=NULL; // init
// intro para globaux
int toto=1; // on est en cinétique par défaut
Grandeur_scalaire_entier grand_courant(toto);
TypeQuelconque typQ1(GENERIQUE_UNE_GRANDEUR_GLOBALE,NU_DDL,grand_courant);
ParaGlob::param->Ajout_grandeur_consultable(&typQ1,AMOR_CINET_VISQUEUX);
// lecture des paramètres attachés au type de calcul
switch (entreePrinc.Lec_ent_info())
{ case 0 :
{lecture_Parametres(entreePrinc); break;}
case -11 : // cas de la création d'un fichier de commande
{ Info_commande_parametres(entreePrinc); break;}
case -12 : // cas de la création d'un schéma XML, on ne fait rien à ce niveau
{ break;}
default:
Sortie(1);
};
};
// constructeur de copie
AlgoriRelaxDyna::AlgoriRelaxDyna (const AlgoriRelaxDyna& algo):
Algori(algo)
,cL_a_chaque_iteration(algo.cL_a_chaque_iteration)
,Ass1_(NULL),Ass2_(NULL),Ass3_(NULL)
,delta_t(0.),typeCalRelaxation(algo.typeCalRelaxation)
,lambda(algo.lambda),alph(algo.alph),beta(algo.beta),gamma(algo.gamma),theta(algo.theta)
,lambda_initial(algo.lambda_initial)
// pilotage automatique éventuel
,lambda_min(algo.lambda_min),lambda_max(algo.lambda_max),delta_lambda(algo.delta_lambda)
,pilotage_auto_lambda(algo.pilotage_auto_lambda),list_iter_relax(algo.list_iter_relax)
,niveauLambda_grandeurGlobale(NULL),niveauLambda_ddlEtendu(NULL),niveauLambda_temps(NULL)
,casMass_relax(algo.casMass_relax)
,option_recalcul_mass(algo.option_recalcul_mass)
,nom_fct_nD_option_recalcul_mass(algo.nom_fct_nD_option_recalcul_mass)
,fct_nD_option_recalcul_mass(algo.fct_nD_option_recalcul_mass)
,ncycle_calcul(algo.ncycle_calcul),type_calcul_mass(algo.type_calcul_mass)
,fac_epsilon(algo.fac_epsilon),type_activation_contact(algo.type_activation_contact)
,proportion_cinetique(algo.proportion_cinetique),visqueux_activer(algo.visqueux_activer)
,et_recalcul_masse_a_la_transition(algo.et_recalcul_masse_a_la_transition)
,niveauF_grandeurGlobale(NULL),niveauF_ddlEtendu(NULL),niveauF_temps(NULL)
,choix_mini_masse_nul(algo.choix_mini_masse_nul)
,cas_combi_ddl(),icas()
,compteur(0),prepa_avec_remont(false)
,brestart(false),type_incre(OrdreVisu::PREMIER_INCRE)
,compteur_demarrage(algo.compteur_demarrage)
,vglobin(),vglobex(),vglobaal(),vcontact()
,save_X_t()
,X_Bl(),V_Bl(),G_Bl(),forces_vis_num(0)
,li_gene_asso(),t_assemb(),tenuXVG(),mat_masse(NULL),mat_masse_sauve(NULL),mat_C_pt(NULL)
,v_mass(algo.v_mass),v_mass1(algo.v_mass1),matglob(NULL)
{//vglobaal = &vglobin;
// -- cas d'un lambda piloté
if (algo.niveauLambda_grandeurGlobale != NULL)
niveauLambda_grandeurGlobale = new Ponderation_GGlobal(*algo.niveauLambda_grandeurGlobale);
if (algo.niveauLambda_ddlEtendu != NULL)
niveauLambda_ddlEtendu = new Ponderation(*algo.niveauLambda_ddlEtendu);
if (algo.niveauLambda_temps != NULL)
niveauLambda_temps = new Ponderation_temps(*algo.niveauLambda_temps);
// -- cas de l'algo mixte : pilotage du basculement entre cinétique et visqueux
if (algo.niveauF_grandeurGlobale != NULL)
niveauF_grandeurGlobale = new Ponderation_GGlobal(*algo.niveauF_grandeurGlobale);
if (algo.niveauF_ddlEtendu != NULL)
niveauF_ddlEtendu = new Ponderation(*algo.niveauF_ddlEtendu);
if (algo.niveauF_temps != NULL)
niveauF_temps = new Ponderation_temps(*algo.niveauF_temps);
};
// destructeur
AlgoriRelaxDyna::~AlgoriRelaxDyna ()
{ if (mat_masse != NULL) delete mat_masse;
if (mat_masse_sauve != NULL) delete mat_masse_sauve;
if (mat_C_pt != NULL) delete mat_C_pt;
if (Ass1_ != NULL) delete Ass1_;
if (Ass2_ != NULL) delete Ass2_;
if (Ass3_ != NULL) delete Ass3_;
if (matglob != NULL) delete matglob;
// -- cas d'un lambda piloté
if (niveauLambda_grandeurGlobale != NULL)
delete niveauLambda_grandeurGlobale;
if (niveauLambda_ddlEtendu != NULL)
delete niveauLambda_ddlEtendu;
if (niveauLambda_temps != NULL)
delete niveauLambda_temps;
// cas des pondérations
if (niveauF_grandeurGlobale != NULL)
delete niveauF_grandeurGlobale;
if (niveauF_ddlEtendu != NULL)
delete niveauF_ddlEtendu;
if (niveauF_temps != NULL)
delete niveauF_temps;
};
// execution de l'algorithme dans le cas dynamique explicite, sans contact
void AlgoriRelaxDyna::Execution(ParaGlob * paraGlob,LesMaillages * lesMail
,LesReferences* lesRef,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD
,VariablesExporter* varExpor,LesLoisDeComp* lesLoisDeComp, DiversStockage* divStock
,Charge* charge,LesCondLim* lesCondLim,LesContacts* lesContacts,Resultats* resultats)
{ Tableau < Fonction_nD* > * tb_combiner = NULL; // ici ne sert pas
// on définit le type de calcul a effectuer :
if ( soustypeDeCalcul->size()==0 )
// cas où il n'y a pas de sous type, on fait le calcul d'équilibre classique
// signifie que le type principal est forcément valide
{ // initialisation du calcul : deux cas, soit avec une lecture initiale du .info, soit une lecture secondaire
if (paraGlob->EtatDeLaLecturePointInfo() == 0)
{InitAlgorithme(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );}
else {MiseAJourAlgo(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
};
// on ne continue que si on n'a pas dépasser le nombre d'incréments maxi ou le temps maxi
// bref que l'on n'a pas fini, sinon on passe
if (! (charge->Fin(icharge,true) ) )
{ // calcul de l'équilibre
CalEquilibre(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats
,tb_combiner);
// fin du calcul, pour l'instant on ne considère pas les autres sous-types
FinCalcul(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
};
}
else
{if ( avec_typeDeCalcul )
// cas où le type principal est valide et qu'il y a des sous_types
{ // on regarde si le sous-type "commandeInteractive" existe, si oui on le met en place
// détermine si le sous type de calcul existe et s'il est actif
if (paraGlob->SousTypeCalcul(commandeInteractive))
{// -- cas avec commandes interactives
// initialisation du calcul : deux cas, soit avec une lecture initiale du .info, soit une lecture secondaire
if (paraGlob->EtatDeLaLecturePointInfo() == 0)
{InitAlgorithme(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
}
else {MiseAJourAlgo(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
};
// calcul de l'équilibre tant qu'il y a des commandes
while (ActionInteractiveAlgo())
{ // on ne continue que si on n'a pas dépasser le nombre d'incréments maxi ou le temps maxi
// bref que l'on n'a pas fini, sinon on passe
if (! (charge->Fin(icharge,true) ) )
CalEquilibre(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats
,tb_combiner);
};
// fin du calcul, pour l'instant on ne considère pas les autres sous-types
FinCalcul(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
}
else // cas sans commandes interactives
{// on fait le calcul d'équilibre
// initialisation du calcul : deux cas, soit avec une lecture initiale du .info, soit une lecture secondaire
if (paraGlob->EtatDeLaLecturePointInfo() == 0)
{InitAlgorithme(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
}
else {MiseAJourAlgo(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
};
// on ne continue que si on n'a pas dépasser le nombre d'incréments maxi ou le temps maxi
// bref que l'on n'a pas fini, sinon on passe
if (! (charge->Fin(icharge,true) ) )
{ // calcul de l'équilibre
CalEquilibre(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats
,tb_combiner);
// fin du calcul, pour l'instant on ne considère pas les autres sous-types
FinCalcul(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
};
// ensuite on teste en fonction des calculs complémentaires
// dépendant des sous_types. Pour l'instant ici uniquement la remontée
list <EnumSousTypeCalcul>::const_iterator ili,ili_fin = soustypeDeCalcul->end();
list <bool>::const_iterator ila;
for (ili = soustypeDeCalcul->begin(),ila = avec_soustypeDeCalcul->begin();
ili!=ili_fin;ili++,ila++)
if (*ila) // cas où le sous type est valide
{if (Remonte_in(*ili)) // on test la présence du calcul de remonté
{ // certaines initialisations sont nécessaires car c'est le premier calcul
Algori::InitRemontSigma(lesMail,lesRef,divStock,charge,lesCondLim,lesContacts,resultats);
Algori::InitErreur(lesMail,lesRef,divStock,charge,lesCondLim,lesContacts,resultats);
Algori::RemontSigma(lesMail);
Algori::RemontErreur(lesMail);
}
else if ( (*ili) == sauveMaillagesEnCours )
{ cout << "\n================================================================="
<< "\n| ecriture des maillages en cours en .her et .lis |"
<< "\n================================================================="
<< endl;
// ----- sort les informations sur fichiers
// Affichage des donnees des maillages dans des fichiers dont le nom est construit
// à partir du nom de chaque maillage au format ".her" et ".lis"
lesMail->Affiche_maillage_dans_her_lis(TEMPS_tdt,*lesRef);
};
};
};// fin du cas sans commandes interactives
}
else
// cas ou le type principal n'est pas valide
// on ne fait que le calcul complémentaire
{ list <EnumSousTypeCalcul>::const_iterator ili,ili_fin = soustypeDeCalcul->end();
list <bool>::const_iterator ila;
for (ili = soustypeDeCalcul->begin(),ila = avec_soustypeDeCalcul->begin();
ili!=ili_fin;ili++,ila++)
if (*ila) // cas où le sous type est valide
{if (Remonte_in(*ili)) // on test la présence du calcul de remonté
{ // certaines initialisations sont nécessaires car c'est le premier calcul
Algori::InitRemontSigma(lesMail,lesRef,divStock,charge,lesCondLim,lesContacts,resultats);
Algori::InitErreur(lesMail,lesRef,divStock,charge,lesCondLim,lesContacts,resultats);
Algori::RemontSigma(lesMail);
Algori::RemontErreur(lesMail);
}
else if ( (*ili) == sauveMaillagesEnCours )
{ cout << "\n================================================================="
<< "\n| ecriture des maillages en cours en .her et .lis |"
<< "\n================================================================="
<< endl;
// ----- sort les informations sur fichiers
// Affichage des donnees des maillages dans des fichiers dont le nom est construit
// à partir du nom de chaque maillage au format ".her" et ".lis"
lesMail->Affiche_maillage_dans_her_lis(TEMPS_0,*lesRef);
};
};
}
}
// si on a forcé la sortie des itérations et incréments, il faut réinitialiser l'indicateur
if (!(pa.EtatSortieEquilibreGlobal()))
pa.ChangeSortieEquilibreGlobal(false);
};
// écriture des paramètres dans la base info
// = 1 : on écrit tout
// = 2 : on écrot uniquement les données variables (supposées comme telles)
void AlgoriRelaxDyna::Ecrit_Base_info_Parametre(UtilLecture& entreePrinc,const int& cas)
{ // récup du flot
ofstream& sort = *entreePrinc.Sort_BI();
// sort << "\n parametres_algo_specifiques_ "<< Nom_TypeCalcul(this->TypeDeCalcul());
switch (cas)
{case 1 : // ------- on sauvegarde tout -------------------------
{ // ecriture des paramètres de contrôle
sort << "\n typeCalRelaxation= " << typeCalRelaxation;
// en fait on sort tous les paramètres de contrôle qu'ils servent ou non, ce qui permettra éventuellement de faire
// des changements de type de relaxation en cours de calcul
sort << "\n lambda= " << lambda << " type_calcul_mass= " << type_calcul_mass;
// option de recalcul de masse
{int nb_option = option_recalcul_mass.Taille(); // nombre d'options
sort << " nombre_option_option_recalcul_mass : " << nb_option;
for (int i=1;i<= nb_option;i++)
{ if (nom_fct_nD_option_recalcul_mass(i).length())
{sort << "\n nom_fct_nD_option_recalcul_mass("<<i<<")= "<< nom_fct_nD_option_recalcul_mass(i);}
else
{sort << "\n valeur_fixe_option_recalcul_mass("<<i<<")= "<< option_recalcul_mass(i); }
};
};
// suite des paramètres
sort << " alpha= " << alph << " beta= " << beta << " gamma= " << gamma << " theta= " << theta
<< " casMass_relax= " << casMass_relax
<< " fac_epsilon= " << fac_epsilon << " type_activation_contact= " << type_activation_contact
<< " choix_mini_masse_nul= "<< choix_mini_masse_nul
<< " proportion_cinetique= " << proportion_cinetique
<< " et_recalcul_masse_a_la_transition= " << et_recalcul_masse_a_la_transition;
// écriture d'un pilotage de lambda éventuelle
sort << "\n lambda_min= " << lambda_min << " lambda_max= " << lambda_max
<< " delta_lambda " << delta_lambda
<< " pilotage_auto_lambda= " << pilotage_auto_lambda ;
// pondération éventuelle
bool sans_courbe = true;
if (niveauLambda_temps == NULL)
{sort << "\n ponderationTemps_lambda 0 ";}
else
{sort << "\n ponderationTemps_lambda 1 " ;
niveauLambda_temps->Ecriture_base_info(sort,cas,sans_courbe);
};
if (niveauLambda_grandeurGlobale == NULL)
{sort << "\n ponderationGrandeurGlobale_lambda 0 ";}
else
{sort << "\n ponderationGrandeurGlobale_lambda 1 " ;
niveauLambda_grandeurGlobale->Ecriture_base_info(sort,cas,sans_courbe);
};
// écriture des pondérations éventuelles
if (niveauF_temps == NULL)
{sort << "\n ponderationTemps_prop_cinetique 0 ";}
else
{sort << "\n ponderation_prop_cinetique 1 " ;
niveauF_temps->Ecriture_base_info(sort,cas,sans_courbe);
};
if (niveauF_grandeurGlobale == NULL)
{sort << "\n ponderationGrandeurGlobale_prop_cinetique 0 ";}
else
{sort << "\n ponderationGrandeurGlobale_prop_cinetique 1 " ;
niveauF_grandeurGlobale->Ecriture_base_info(sort,cas,sans_courbe);
};
sort << " ncycle_calcul= " << ncycle_calcul << " "
<< "\n cL_a_chaque_iteration "<< cL_a_chaque_iteration<< " ";
sort << flush;
break;
}
case 2 : // ----------- on sauvegarde en minimaliste --------------------
{ sort << "\n lambda= "<<lambda << " ";
break;
}
default :
{ cout << "\nErreur : valeur incorrecte du type de sauvegarde !\n";
cout << "AlgoriRelaxDyna::Ecrit_Base_info_Parametre(UtilLecture& ,const int& )"
<< " cas= " << cas << endl;
Sortie(1);
}
};
// sort << "\n fin_parametres_algo_specifiques_ ";
};
// lecture des paramètres dans la base info
// = 1 : on récupère tout
// = 2 : on récupère uniquement les données variables (supposées comme telles)
// choix = true : fonctionememt normal
// choix = false : la méthode ne doit pas lire mais initialiser les données à leurs valeurs par défaut
// car la lecture est impossible
void AlgoriRelaxDyna::Lecture_Base_info_Parametre(UtilLecture& entreePrinc,const int& cas,bool choix)
{// récup du flot
ifstream& ent = *entreePrinc.Ent_BI();
string toto;
if (choix)
{switch (cas)
{ case 1 : // ------- on récupère tout -------------------------
{ // lecture du type de calcul de la relaxation
ent >> toto >> typeCalRelaxation;
switch (typeCalRelaxation)
{ case 1: visqueux_activer = false; break;
case 2: visqueux_activer = true; break;
case 4: visqueux_activer = false; break; // au début en cinétique
default:
cout <<"\n typeCalRelaxation = " << typeCalRelaxation
<< " cas non pris en compte dans AlgoriRelaxDyna::Lecture_Base_info_Parametre( ";
Sortie(1);
break;
};
Transfert_ParaGlob_AMOR_CINET_VISQUEUX();
// en fait on lit tous les paramètres de contrôle qu'ils servent ou non, ce qui permettra éventuellement de faire
// des changements de type de relaxation en cours de calcul
ent >> toto >> lambda >> toto >> type_calcul_mass;
// option de recalcul de masse
{int nb_option = 0;
ent >> toto >> nb_option;
option_recalcul_mass.Change_taille(nb_option);
nom_fct_nD_option_recalcul_mass.Change_taille(nb_option);
fct_nD_option_recalcul_mass.Change_taille(nb_option);
for (int i=1;i<= nb_option;i++)
{ ent >> toto; // on lit l'identificateur
string iden("nom_fct_nD_option_recalcul_mass");
std::size_t found = toto.find(iden);
if (found!=std::string::npos) // si oui on a trouvé la chaine
{ ent >> nom_fct_nD_option_recalcul_mass(i);}
else // sinon cela signifie qu'il s'agit d'une valeur fixe
{ ent >> option_recalcul_mass(i); };
// dans tous les cas on met à NULL pour l'instant le pointeur de fonction
fct_nD_option_recalcul_mass(i) = NULL;
};
};
// suite des paramètres
ent >> toto >> alph >> toto >> beta >> toto >> gamma >> toto >> theta
>> toto >> casMass_relax
>> toto >> fac_epsilon >> toto >> type_activation_contact
>> toto >> choix_mini_masse_nul >> toto >> proportion_cinetique
>> toto >> et_recalcul_masse_a_la_transition;
// cas d'un pilotage de lambda éventuelle
ent >> toto >> lambda_min >> toto >> lambda_max
>> toto >> delta_lambda
>> toto >> pilotage_auto_lambda ;
// avec pondération éventuelle
// -> le temps
bool ponder_lambda_temps = false;
ent >> toto >> ponder_lambda_temps;
if (ponder_lambda_temps)
{ if (niveauLambda_temps == NULL)
{niveauLambda_temps = new Ponderation_temps();
niveauLambda_temps->Lecture_base_info(ent,cas);
}
else {niveauLambda_temps->Lecture_base_info(ent,cas);};
};
// -> les grandeurs globales
bool ponder_lambda_glob = false;
ent >> toto >> ponder_lambda_glob;
if (ponder_lambda_glob)
{ if (niveauLambda_grandeurGlobale == NULL)
{niveauLambda_grandeurGlobale = new Ponderation_GGlobal();
niveauLambda_grandeurGlobale->Lecture_base_info(ent,cas);
}
else {niveauLambda_grandeurGlobale->Lecture_base_info(ent,cas);};
};
// cas des pondérations éventuelles
// -> le temps
bool ponder_temps = false;
ent >> toto >> ponder_temps;
if (ponder_temps)
{ if (niveauF_temps == NULL)
{niveauF_temps = new Ponderation_temps();
niveauF_temps->Lecture_base_info(ent,cas);
}
else {niveauF_temps->Lecture_base_info(ent,cas);};
};
// -> les grandeurs globales
bool ponder_glob = false;
ent >> toto >> ponder_glob;
if (ponder_glob)
{ if (niveauF_grandeurGlobale == NULL)
{niveauF_grandeurGlobale = new Ponderation_GGlobal();
niveauF_grandeurGlobale->Lecture_base_info(ent,cas);
}
else {niveauF_grandeurGlobale->Lecture_base_info(ent,cas);};
};
ent >> toto >> ncycle_calcul ;
// lecture de cL_a_chaque_iteration
ent >> toto ;
if (toto != "cL_a_chaque_iteration")
{ cout << "\n erreur en lecture du parametre cL_a_chaque_iteration"
<< "\n on attendait le mot : cL_a_chaque_iteration , au lieu de " << toto
<< "\n AlgoriRelaxDyna::Lecture_Base_info_Parametre( ... ";
Sortie(1);
}
ent >> cL_a_chaque_iteration ;
break;
}
case 2 : // ----------- on récupère en minimaliste --------------------
{ double titi;
ent >> toto >> titi; // on ne relit pas le lambda ceci pour pouvoir le modifier
// au moment d'un restart par exemple
break;
}
default :
{ cout << "\nErreur : valeur incorrecte du type de sauvegarde !\n";
cout << "AlgoriRelaxDyna::Lecture_Base_info_Parametre(UtilLecture& ,const int& )"
<< " cas= " << cas << endl;
Sortie(1);
}
};
}
//---- dans le cas où un restart en relaxation n'est pas prévu, car par exemple cela suit l'utilisation
// d'un autre algo, on ne fait rien, et on se contente de garder les paramètres qui viennent d'être lues
// ---> décision du 21 mars 2019 .... à voir si cela ne crée pas de pb
// else
// { // cas où la lecture est impossible
// typeCalRelaxation=1; lambda = 0.605; lambda_initial=lambda;
// type_calcul_mass= 2;
// option_recalcul_mass.Change_taille(1,0);
// nom_fct_nD_option_recalcul_mass.Change_taille(1,"");
// fct_nD_option_recalcul_mass.Change_taille(1,NULL);
// alph=1; beta = 1.;gamma=1.;theta=1.; casMass_relax=3;
// fac_epsilon= 1.; type_activation_contact=1;choix_mini_masse_nul=0;proportion_cinetique = 0.1;
// et_recalcul_masse_a_la_transition = 1;
// };
};
// gestion et vérification du pas de temps et modif en conséquence si nécessaire
// cas = 0: premier passage à blanc, delta t = 0
// cas = 1: initialisation du pas de temps et vérif / au pas de temps critique
// ceci pour le temps t=0
// cas = 2: initialisation du pas de temps et vérif / au pas de temps critique
// ceci pour le temps t
// en entrée: modif_pas_de_temps: indique qu'il y a eu par ailleurs (via Charge->Avance())
// une modification du pas de temps depuis le dernier appel
// retourne vrai s'il y a une modification du pas de temps, faux sinon
bool AlgoriRelaxDyna::Gestion_pas_de_temps(bool modif_pas_de_temps,LesMaillages * lesMail,int cas)
{ bool modif_deltat = modif_pas_de_temps; // booleen pour la prise en compte éventuelle de la modif du temps éventuelle
if (modif_pas_de_temps) // dans le cas où il y a eu une modification externe du pas de temps on modifie la variable interne
delta_t = pa.Deltat(); // sinon elle ne sera pas mise à jour dans l'algo
// double ancien_pas = delta_t;
switch (cas)
{ case 0 :
{ // cas d'un passage à blanc, rien ne bouge et delta t est mis à 0
delta_t = 0.;
break;
}
case 1:
{ // dans le cas où le pas de temps ou le pas de temps maxi dépendent d'un pas critique,
// on calcul le pas de temps critique et on met à jour les pas de temps
if (pa.DeltatOuDeltatmaxDependantTempsCritique())
// on calcul le pas de temps minimal pour cela on utilise
// les caractéristiques dynamiques d'une biellette de longueur
// valant le minimum d'un coté d'arrête
{ // ---- on calcul le pas de temps minimal
double l_sur_c = lesMail->Longueur_arrete_mini_sur_c(TEMPS_0);
double delta_t_essai = l_sur_c ;
// mise à jour éventuel du pas de temps et du pas de temps maxi et mini s'ils sont définit à partir du temps critique
modif_deltat = pa.Modif_Deltat_DeltatMaxi(delta_t_essai,l_sur_c); // mais pas de test vis-a-vis des bornes
delta_t = pa.Deltat(); // mise à jour éventuelle du pas de temps
};
break;
}
case 2:
{ // dans le cas où le pas de temps ou le pas de temps maxi dépendent d'un pas critique,
// on calcul le pas de temps critique et on met à jour les pas de temps
if (pa.DeltatOuDeltatmaxDependantTempsCritique())
// on calcul le pas de temps minimal pour cela on utilise
// les caractéristiques dynamiques d'une biellette de longueur
// valant le minimum d'un coté d'arrête
{ // ---- on calcul le pas de temps minimal
double l_sur_c = lesMail->Longueur_arrete_mini_sur_c(TEMPS_t);
double delta_t_essai = l_sur_c ;
// mise à jour éventuel du pas de temps et du pas de temps maxi et mini s'ils sont définit à partir du temps critique
modif_deltat = pa.Modif_Deltat_DeltatMaxi(delta_t_essai,l_sur_c); // mais pas de test vis-a-vis des bornes
delta_t = pa.Deltat(); // mise à jour éventuelle du pas de temps
// // on vérifie que le pas de temps est correcte: on propose le deltat actuel qui sera modifié s'il n'est pas bon
// double deltat_actuel = pa.Deltat();
// pa.Modif_Detat_dans_borne(deltat_actuel);
};
break;
} //-- fin du cas == 2
default :
{cout << "\nErreur : valeur incorrecte du cas = " << cas << "\n";
cout << "AlgoriRelaxDyna::Gestion_pas_de_temps(... \n";
Sortie(1);
}
}; // fin du switch
// if (modif_deltat)
// cout << "\n --->>>> modif increment de temps de " << ancien_pas << " a " << delta_t;
// retour
return modif_deltat;
};
// lecture des paramètres du calcul
void AlgoriRelaxDyna::lecture_Parametres(UtilLecture& entreePrinc)
{ MotCle motCle; // ref aux mots cle
Transfert_ParaGlob_ALGO_GLOBAL_ACTUEL(RELAX_DYNA); // transfert info
deja_lue_entete_parametre = 1; // a priori pas de lecture d'entête
// on se positionne sur le prochain mot clé
do
{ entreePrinc.NouvelleDonnee();
}
while ( !motCle.SimotCle(entreePrinc.tablcar)) ;
// si le mot clé est "PARA_TYPE_DE_CALCUL" cela signifie
// qu'il y a un paramètre à lire
bool lecture_effective = false; // indiquera si oui ou non on a lue quelque chose
// donc que l'on doit passer à la ligne qui suit
string nom_class_methode("AlgoriRelaxDyna::lecture_Parametres"); // pour la lecture
int min=1,max=0;string mot_cle;int int_defaut=0; double double_defaut=0;// variables inter qui servent pour la lecture
double double_min = 1.; double double_max = 0.;
if (strstr(entreePrinc.tablcar,"PARA_TYPE_DE_CALCUL")!=NULL)
{ //cas de la définition de paramètres
// on signale à Algori qu'il y a eu déjà une lecture de paramètre
deja_lue_entete_parametre=2;
// --- lecture facultative du premier paramètres de l'algorithme:
entreePrinc.NouvelleDonnee(); // ligne suivante
int_defaut=1; mot_cle="typeCalRelaxation=";min=1;max=4;
{bool lec_eff = entreePrinc.Lecture_un_parametre_int(int_defaut,nom_class_methode,min,max,mot_cle,typeCalRelaxation);
lecture_effective = lecture_effective || lec_eff;
// dans le cas où on va utiliser une relaxation avec amortissement visqueux on l'indique au niveau des paramètres généraux
if ((typeCalRelaxation ==2 )||(typeCalRelaxation ==3 )||(typeCalRelaxation ==4 ))
{pa.Change_amort_visco_artificielle(4);};
};
// --- lecture facultative du paramètres lambda :
double_defaut=0.605; mot_cle="lambda=";
min=1;max=0; // le max < min ==> pas de prise en compte de valeurs
{bool lec_eff = entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,min,max,mot_cle,lambda);
lecture_effective = lecture_effective || lec_eff;
lambda_initial=lambda;
};
// --- lecture facultative du paramètres type_calcul_mass :
int_defaut=2; mot_cle="type_calcul_mass="; min=1;max=2;
{bool lec_eff = entreePrinc.Lecture_un_parametre_int(int_defaut,nom_class_methode,min,max,mot_cle,type_calcul_mass);
lecture_effective = lecture_effective || lec_eff;
};
// --- lecture éventuelle de l'option_recalcul_mass
if (strstr(entreePrinc.tablcar,"option_recalcul_mass=")!=NULL)
{ // à ce niveau on ne lit a priori qu'une valeur
// 1- init par défaut
option_recalcul_mass.Change_taille(1,0);
nom_fct_nD_option_recalcul_mass.Change_taille(1,"");
fct_nD_option_recalcul_mass.Change_taille(1,NULL);
int_defaut=0; mot_cle="option_recalcul_mass="; min=-1;max=5;
// 2- on regarde s'il y a une présence d'une fonction nD
// pour cela on lit la première chaine de caractère
string nom_inter;
bool lec_eff = entreePrinc.Lecture_mot_cle_et_string(nom_class_methode,mot_cle,nom_inter);
lecture_effective = lecture_effective || lec_eff;
// maintenant on scrute le nom_inter
string iden("nom_fct_nD_option_recalcul_mass_");
std::size_t found = nom_inter.find(iden);
if (found!=std::string::npos) // si oui on a trouvé la chaine
// on lit le nom de la fonction
{ *(entreePrinc.entree) >> nom_fct_nD_option_recalcul_mass(1);}
else // sinon cela signifie qu'il s'agit d'une valeur fixe
{option_recalcul_mass(1) = ChangeEntier(nom_inter);};
// dans tous les cas on met à NULL pour l'instant le pointeur de fonction
fct_nD_option_recalcul_mass(1) = NULL;
};
// --- lecture éventuelle de l'option_recalcul_mass spécifiquement pour l'amortissement visqueux
int_defaut=0; mot_cle="opt_visqueux_recal_mass="; min=-1;max=5;
{// on est obligé de regardé si l'indicateur existe avant toute chose
if (strstr(entreePrinc.tablcar,"opt_visqueux_recal_mass=")!=NULL)
{// 1) là if faut redimentionner les tableaux
option_recalcul_mass.Change_taille(2);
nom_fct_nD_option_recalcul_mass.Change_taille(2);
fct_nD_option_recalcul_mass.Change_taille(2);
// 2- on regarde s'il y a une présence d'une fonction nD
// pour cela on lit la première chaine de caractère
string nom_inter;
bool lec_eff = entreePrinc.Lecture_mot_cle_et_string(nom_class_methode,mot_cle,nom_inter);
lecture_effective = lecture_effective || lec_eff;
// maintenant on scrute le nom_inter
string iden("nom_fct_nD_option_recalcul_mass_");
std::size_t found = nom_inter.find(iden);
if (found!=std::string::npos) // si oui on a trouvé la chaine
// on lit le nom de la fonction
{ *(entreePrinc.entree) >> nom_fct_nD_option_recalcul_mass(2);}
else // sinon cela signifie qu'il s'agit d'une valeur fixe
{ option_recalcul_mass(2) = ChangeEntier(nom_inter); };
// dans tous les cas on met à NULL pour l'instant le pointeur de fonction
fct_nD_option_recalcul_mass(2) = NULL;
};
};
// ----------- maintenant on va lire les paramètres de contrôle éventuelles ------------
if (lecture_effective)
{entreePrinc.NouvelleDonnee(); // si on a lue des para sur la première ligne on passe
lecture_effective=false;
};
while ( (strstr(entreePrinc.tablcar,"parametre_recalcul_de_la_masse_")!=NULL)
|| (strstr(entreePrinc.tablcar,"parametre_calcul_de_la_masse_")!=NULL)
|| (strstr(entreePrinc.tablcar,"parametre_activation_du_contact_")!=NULL)
|| (strstr(entreePrinc.tablcar,"choix_mini_masse_nul=")!= NULL)
|| (strstr(entreePrinc.tablcar,"proportion_cinetique=")!= NULL)
|| (strstr(entreePrinc.tablcar,"propCinetiqueAvecPonderation_temps_")!=NULL)
|| (strstr(entreePrinc.tablcar,"propCinetiqueAvecPonderation_Globale_")!=NULL)
|| (strstr(entreePrinc.tablcar,"lambdaAvecPonderation_temps_")!=NULL)
|| (strstr(entreePrinc.tablcar,"lambdaAvecPonderation_Globale_")!=NULL)
|| (strstr(entreePrinc.tablcar,"pilotageAutolambda_")!=NULL)
|| (strstr(entreePrinc.tablcar,"cL_a_chaque_iteration_")!=NULL)
)
{
//=== paramètres de controle du calcul de la masse =====
if (strstr(entreePrinc.tablcar,"parametre_calcul_de_la_masse_")!=NULL)
{ mot_cle="parametre_calcul_de_la_masse_";
lecture_effective = entreePrinc.Lecture_et_verif_mot_cle(nom_class_methode,mot_cle);
// si on arrive ici c'est que l'on a bien lue le mot clé "parametre_calcul_de_la_masse_"
// donc il faudra passer une ligne de toute manière, donc ce n'est plus la peine de tester les booleen
if (type_calcul_mass == 1)
{ // cas, du modèle d'han lee amélioré
// --- lecture éventuelle de alpha
double_defaut=1.; mot_cle="alpha="; min=1;max=0;
entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,min,max,mot_cle,alph);
// --- lecture éventuelle de beta
double_defaut=1.; mot_cle="beta="; min=1;max=0;
entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,min,max,mot_cle,beta);
// --- lecture éventuelle de gamma
double_defaut=1.; mot_cle="gamma="; min=1;max=0;
entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,min,max,mot_cle,gamma);
// --- lecture éventuelle de theta
double_defaut=1.; mot_cle="theta="; min=1;max=0;
entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,min,max,mot_cle,theta);
// --- lecture éventuelle de casMass_relax
int_defaut=3; mot_cle="casMass_relax="; min=1;max=5;
entreePrinc.Lecture_un_parametre_int(int_defaut,nom_class_methode,min,max,mot_cle,casMass_relax);
}
else if (type_calcul_mass == 2)
{ // cas du modèle de calcul utilisant la matrice de raideur
// --- lecture éventuelle de casMass_relax
int_defaut=3; mot_cle="casMass_relax="; min=0;max=9;
entreePrinc.Lecture_un_parametre_int(int_defaut,nom_class_methode,min,max,mot_cle,casMass_relax);
};
};
if (lecture_effective)
{entreePrinc.NouvelleDonnee(); // si on a lue des para sur la ère ligne on passe à la suivante
lecture_effective=false;
};
//=== paramètres de controle du re-calcul de la masse =====
if (strstr(entreePrinc.tablcar,"parametre_recalcul_de_la_masse_")!=NULL)
{ mot_cle="parametre_recalcul_de_la_masse_";
lecture_effective = entreePrinc.Lecture_et_verif_mot_cle(nom_class_methode,mot_cle);
// si on arrive ici c'est que l'on a bien lue le mot clé "parametre_calcul_de_la_masse_"
// donc il faudra passer une ligne de toute manière, donc ce n'est plus la peine de tester les booleen
// --- lecture éventuelle de fac_epsilon
double_defaut=1.; mot_cle="fac_epsilon="; min=1;max=0;
entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,min,max,mot_cle,fac_epsilon);
// --- lecture éventuelle de ncycle_calcul
double_defaut=100; mot_cle="ncycle_calcul="; min=1;max=ConstMath::grand;
entreePrinc.Lecture_un_parametre_int(double_defaut,nom_class_methode,min,max,mot_cle,ncycle_calcul);
// --- vérification de l'existence éventuelle du mot clé: et_recalcul_masse_a_la_transition_
if (strstr(entreePrinc.tablcar,"et_pas_recalcul_masse_a_la_transition_")!=NULL)
{// si oui, on lit le mot clé
mot_cle="et_pas_recalcul_masse_a_la_transition_";
lecture_effective = entreePrinc.Lecture_et_verif_mot_cle(nom_class_methode,mot_cle);
et_recalcul_masse_a_la_transition = 0; // on enregistre
};
};
if (lecture_effective)
{entreePrinc.NouvelleDonnee(); // si on a lue des para sur la ère ligne on passe à la suivante
lecture_effective=false;
};
//=== paramètres de controle automatique de l'évolution de lambda =====
if (strstr(entreePrinc.tablcar,"pilotageAutolambda_")!=NULL)
{ mot_cle="pilotageAutolambda_";
lecture_effective = entreePrinc.Lecture_et_verif_mot_cle(nom_class_methode,mot_cle);
pilotage_auto_lambda = true;
// si on arrive ici c'est que l'on a bien lue le mot clé "pilotageAutolambda_"
// donc il faudra passer une ligne de toute manière, donc ce n'est plus la peine de tester les booleen
// --- lecture éventuelle de lambda_min
double xmin = 0.5; double xmax = 50.;
double_defaut=0.605; mot_cle="lambda_min=";
entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,xmin,xmax,mot_cle,lambda_min);
// --- lecture éventuelle de lambda_max
xmin = 0.5; xmax = ConstMath::grand;
double_defaut=ConstMath::grand; mot_cle="lambda_max="; ;
entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,xmin,xmax,mot_cle,lambda_max);
// --- lecture éventuelle de delta_lambda
xmin = 0.; xmax = ConstMath::grand;
double_defaut=0.1; mot_cle="delta_lambda="; ;
entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,xmin,xmax,mot_cle,delta_lambda);
};
if (lecture_effective)
{entreePrinc.NouvelleDonnee(); // si on a lue des para sur la ère ligne on passe à la suivante
lecture_effective=false;
};
//=== paramètres de controle de l'activation du contact =====
if (strstr(entreePrinc.tablcar,"parametre_activation_du_contact_")!=NULL)
{ mot_cle="parametre_activation_du_contact_";
lecture_effective = entreePrinc.Lecture_et_verif_mot_cle(nom_class_methode,mot_cle);
// si on arrive ici c'est que l'on a bien lue le mot clé "parametre_activation_du_contact_"
// donc il faudra passer une ligne de toute manière, donc ce n'est plus la peine de tester les booleen
// --- lecture éventuelle de type_activation_contact
int_defaut=1; mot_cle="type_activation_contact="; min=0;max=2;
entreePrinc.Lecture_un_parametre_int(int_defaut,nom_class_methode,min,max,mot_cle,type_activation_contact);
};
//=== paramètres de controle d'une masse nulle =====
if (strstr(entreePrinc.tablcar,"choix_mini_masse_nul=")!=NULL)
{ // --- lecture éventuelle de choix_mini_masse_nul
int_defaut=1; mot_cle="choix_mini_masse_nul="; min=0;max=2;
lecture_effective = entreePrinc.Lecture_un_parametre_int(int_defaut,nom_class_methode,min,max,mot_cle,choix_mini_masse_nul);
};
//=== paramètres de controle de la proportion de cinétique / visqueux pour le type 4 =====
if (strstr(entreePrinc.tablcar,"proportion_cinetique=")!=NULL)
{ // --- lecture éventuelle de proportion_cinetique
double_defaut=0.1; mot_cle="proportion_cinetique="; double_min=0;double_max=1.;
lecture_effective = entreePrinc.Lecture_un_parametre_double
(double_defaut,nom_class_methode,double_min,double_max,mot_cle,proportion_cinetique);
};
//=== paramètres de controle de la proportion de cinétique / visqueux pour le type 4 =====
// utilisation d'une dépendance au temps éventuelle: on lit le nom de la fonction
if (strstr(entreePrinc.tablcar,"propCinetiqueAvecPonderation_temps_")!=NULL)
{ // --- lecture du nom de la fonction de dépendance au temps
mot_cle="propCinetiqueAvecPonderation_temps_";
niveauF_temps = new Ponderation_temps();
string nom_fonction; // init
lecture_effective = entreePrinc.Lecture_mot_cle_et_string
(nom_class_methode,mot_cle,nom_fonction);
niveauF_temps->Nom_fonction() = nom_fonction; // affectation du nom de la fonction
};
//=== paramètres de controle de la proportion de cinétique / visqueux pour le type 4 =====
// utilisation d'une dépendance à des variables globales éventuelles: on lit le nom de la fonction
if (strstr(entreePrinc.tablcar,"propCinetiqueAvecPonderation_Globale_")!=NULL)
{ // --- lecture du nom de la fonction de dépendance au temps
mot_cle="propCinetiqueAvecPonderation_Globale_";
niveauF_grandeurGlobale = new Ponderation_GGlobal();
string nom_fonction; // init
lecture_effective = entreePrinc.Lecture_mot_cle_et_string
(nom_class_methode,mot_cle,nom_fonction);
niveauF_grandeurGlobale->Nom_fonction() = nom_fonction; // affectation du nom de la fonction
};
//=== paramètres de controle de la proportion de cinétique / visqueux pour le type 4 =====
// utilisation d'une dépendance au temps éventuelle: on lit le nom de la fonction
if (strstr(entreePrinc.tablcar,"lambdaAvecPonderation_temps_")!=NULL)
{ // --- lecture du nom de la fonction de dépendance au temps
mot_cle="lambdaAvecPonderation_temps_";
niveauLambda_temps = new Ponderation_temps();
string nom_fonction; // init
lecture_effective = entreePrinc.Lecture_mot_cle_et_string
(nom_class_methode,mot_cle,nom_fonction);
niveauLambda_temps->Nom_fonction() = nom_fonction; // affectation du nom de la fonction
};
//=== paramètres de controle de la proportion de cinétique / visqueux pour le type 4 =====
// utilisation d'une dépendance à des variables globales éventuelles: on lit le nom de la fonction
if (strstr(entreePrinc.tablcar,"lambdaAvecPonderation_Globale_")!=NULL)
{ // --- lecture du nom de la fonction de dépendance au temps
mot_cle="lambdaAvecPonderation_Globale_";
niveauLambda_grandeurGlobale = new Ponderation_GGlobal();
string nom_fonction; // init
lecture_effective = entreePrinc.Lecture_mot_cle_et_string
(nom_class_methode,mot_cle,nom_fonction);
niveauLambda_grandeurGlobale->Nom_fonction() = nom_fonction; // affectation du nom de la fonction
};
if (mot_cle == "cL_a_chaque_iteration_")
{ // on lit le paramètre
*(entreePrinc.entree) >> cL_a_chaque_iteration;
lecture_effective = true;
};
if (lecture_effective)
{entreePrinc.NouvelleDonnee(); // si on a lue des para sur la ère ligne on passe à la suivante
lecture_effective=false;
};
}; // fin du while
}
else // sinon , cas où il n'y a pas de paramètres de calcul, on met les valeurs par défaut
{ typeCalRelaxation=1; lambda = 0.605;lambda_initial=lambda;
type_calcul_mass= 2;
option_recalcul_mass.Change_taille(1,0);
nom_fct_nD_option_recalcul_mass.Change_taille(1,"");
fct_nD_option_recalcul_mass.Change_taille(1,NULL);
alph=1; beta = 1.;gamma=1.;theta=1.; casMass_relax=3;
fac_epsilon= 1.;type_activation_contact= 1;choix_mini_masse_nul=0;
proportion_cinetique=0.1;
// pilotage auto éventuel
lambda_min=0.605;lambda_max=50.;delta_lambda=0.1;
pilotage_auto_lambda=false;
cL_a_chaque_iteration = 0;
};
// mise à jour variable globale
switch (typeCalRelaxation)
{ case 1: visqueux_activer = false; break;
case 2: visqueux_activer = true; break;
case 4: visqueux_activer = false; break; // au début en cinétique
default:
cout <<"\n typeCalRelaxation = " << typeCalRelaxation
<< " cas non pris en compte dans AlgoriRelaxDyna::Lecture_Base_info_Parametre( ";
Sortie(1);
break;
};
Transfert_ParaGlob_AMOR_CINET_VISQUEUX();
// on prépare la prochaine lecture si la lecture a été effective et que l'on n'est pas
// sur un mot clé
if ((lecture_effective) && ( !motCle.SimotCle(entreePrinc.tablcar)))
entreePrinc.NouvelleDonnee(); // ligne suivante
// puis appel de la méthode de la classe mère
Algori::lecture_Parametres(entreePrinc);
};
// création d'un fichier de commande: cas des paramètres spécifiques
void AlgoriRelaxDyna::Info_commande_parametres(UtilLecture& entreePrinc)
{ // écriture dans le fichier de commande
ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
sort << "\n# -----------------------------------------------------------------------------"
<< "\n#| 1) parametres (falcultatifs ) associes au calcul de relaxation dynamique |"
<< "\n#| avec en continu, modification et optimisation de la matrice masse |"
<< "\n#| 2) un parametre (falcultatif ) specifique au contact |"
<< "\n# -----------------------------------------------------------------------------"
<< "\n"
<< "\n # PARA_TYPE_DE_CALCUL"
<< "\n # ........................................................."
<< "\n # / type d'algorithme /"
<< "\n #........................................................."
<< "\n # plusieurs methodes d'amortissement dynamique sont implantees "
<< "\n # le parametre typeCalRelaxation= indique le choix de la methode : "
<< "\n # "
<< "\n # typeCalRelaxation= 1 : Il s'agit d'une relaxation avec amortissement cinetique, qui utilise differente "
<< "\n # techniques pour le calcul de la pseudo-masse. Ces techniques sont differenciee a l'aide"
<< "\n # de parametres particuliers explicites par la suite. L'amortissement etant cinetique "
<< "\n # on peut donc modifier egalement les parametres de l'amortissement cinetique"
<< "\n # "
<< "\n # typeCalRelaxation= 2 : Il s'agit d'une relaxation avec amortissement visqueux critique, qui utilise les memes "
<< "\n # technique de calcul de la pseudo-masse que le cas typeCalRelaxation= 1 par contre "
<< "\n # l'amortissement est different car visqueux critique. "
<< "\n # "
<< "\n # typeCalRelaxation= 4 : Il s'agit d'une relaxation avec amortissement cinetique au debut du calcul "
<< "\n # puis a partir d'un seuil, le calcul se poursuit avec un amortissement visqueux."
<< "\n # La premiere partie suit la methodologie du cas : typeCalRelaxation= 1 "
<< "\n # La seconde partie suit la methodologie du cas: typeCalRelaxation= 2"
<< "\n # Le seuil est cacule a partir d'une fonction multidimensionnelle des grandeures globales"
<< "\n # exemple syntaxe, a inserer a la fin des parametres de l'algo : "
<< "\n # les_grandeurs_de_controle_= NORME_CONVERGENCE fin_grandeurs_ "
<< "\n # fonctions_changement_amortissement= fct1 fin_fonct_ "
<< "\n # "
<< "\n # "
<< "\n # Pour chaque type de methode, des parametres specifiques sont associes pour controler son comportement "
<< "\n # "
<< "\n # le typeCalRelaxation= est par defaut = 1 "
<< "\n # "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # ------------- typeCalRelaxation= 1 ----------------------- "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # dans le cas ou typeCalRelaxation= 1, il s'agit d'une relaxation avec amortissement cinetique, "
<< "\n # associee avec un type particulier de calcul de la matrice masse "
<< "\n # a) le parametre qui permet de choisir parmis les differents cas est type_calcul_mass "
<< "\n # Son utilisation est decrite par la suite "
<< "\n # b) on peut egalement intervenir sur les parametres qui pilotes l'amortissement cinetique (cf. doc) "
<< "\n # c) Il y a egalement un parametre lambda permettant de ponderer le calcul de la masse "
<< "\n # d) enfin un parametre option_recalcul_mass= permet de choisir le test de re-calcul de la matrice masse "
<< "\n # Exemple de declaration: "
<< "\n ## typeCalRelaxation= 1 lambda= 0.7 type_calcul_mass= 1 option_recalcul_mass= 1 "
<< "\n # "
<< "\n # >>>>> PARTICULARITE DU PARAMETRE: option_recalcul_mass= "
<< "\n # ==> A la place d'une valeur numerique pour option_recalcul_mass= il est possible d'indiquer le "
<< "\n # nom d'une fonction nD precede du mot cle : nom_fct_nD_option_recalcul_mass_ "
<< "\n # on aura alors la syntaxe suivane: "
<< "\n # option_recalcul_mass= nom_fct_nD_option_recalcul_mass_ <un nom de fonction nD> "
<< "\n # dans ce cas le recalcul de la masse suis l'option indiquee par le retour de la fonction nD "
<< "\n # en particulier si elle retourne 0 pendant l'execution, il n'y pas de recalcul de masse "
<< "\n # NB: La fonction doit evidemment exister, elle doit egalement uniquement utiliser des grandeurs"
<< "\n # globales. "
<< "\n # La fonction est appelee en particulier a chaque iteration. "
<< "\n # L'utilisation d'une fonction nD est possible avec tous les types de relaxation "
<< "\n # cinetique, visqueux et egalement mixte "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # ------------- typeCalRelaxation= 2 ----------------------- "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # dans le cas ou typeCalRelaxation= 2, il s'agit d'une relaxation avec amortissement visqueux critique, "
<< "\n # associee avec des types de calcul de la matrice masse different du cas de Barnes "
<< "\n # a) le parametre qui permet de choisir parmis les differents cas de masse est type_calcul_mass "
<< "\n # Son utilisation est decrite par la suite "
<< "\n # b) on peut egalement intervenir sur les parametres qui pilotes l'amortissement visqueux (cf. doc) "
<< "\n # c) Il y a egalement un parametre lambda permetant de ponderer le calcul de la masse "
<< "\n # d) enfin un parametre option_recalcul_mass= permet de choisir le test de re-calcul de la matrice masse "
<< "\n # Exemple de declaration: "
<< "\n ## typeCalRelaxation= 2 lambda= 0.7 type_calcul_mass= 2 option_recalcul_mass= 4 "
<< "\n # "
<< "\n # tous les parametres sont facultatifs: les valeurs par defaut sont: "
<< "\n # typeCalRelaxation= 2 ; lambda= 0.605 ; type_calcul_mass= 2; option_recalcul_mass= 0; "
<< "\n # "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # ------------- typeCalRelaxation= 4 ----------------------- "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # il s'agit d'une relaxation mixte qui debute en amortissement cinetique et peut terminer en visqueux "
<< "\n # On indique les parametres habituels: lambda , le type de calcul de masse , et son recalcul "
<< "\n # suivent dans un ordre quelconque les parametres facultatifs suivants ( explicites ensuite) : "
<< "\n # parametre_recalcul_de_la_masse_ parametre_calcul_de_la_masse_ "
<< "\n # parametre_activation_du_contact_ , choix_mini_masse_nul= , proportion_cinetique= "
<< "\n # propCinetiqueAvecPonderation_temps_ , propCinetiqueAvecPonderation_Globale_"
<< "\n # lambdaAvecPonderation_temps_ , lambdaAvecPonderation_Globale_ , pilotageAutolambda_ "
<< "\n # "
<< "\n # parmi l'ensemble de ces parametres, il faut noter: "
<< "\n # 1) proportion_cinetique= suivi d'un reel: par defaut vaut 0.1 "
<< "\n # ce parametre indique la proportion (proportion_cinetique) qu'il faut atteindre relativement "
<< "\n # a la precision demandee, pour basculer de l'amortissement cinetique a l'amortissement visqueux "
<< "\n # 2) propCinetiqueAvecPonderation_temps_ suivi d'un nom de fonction 1D "
<< "\n # --> indique le nom d'une fonction 1D du temps, qui permet de faire evoluer le parametre "
<< "\n # proportion_cinetique en fonction du temps : "
<< "\n # le resultat = la valeur initiale de proportion_cinetique * la valeur de la fonction(temps) "
<< "\n # 3) propCinetiqueAvecPonderation_Globale_ suivi d'un nom de fonction nD "
<< "\n # --> indique le nom d'une fonction nD de grandeurs globales, qui permet de faire "
<< "\n # evoluer le parametre proportion_cinetique en fonction de grandeurs globales : "
<< "\n # le resultat = la valeur initiale de proportion_cinetique * la valeur de la fonction nD "
<< "\n # 4) il est possible d'indiquer de ne pas recalculer la masse au moment de la transition "
<< "\n # cinetique - visqueux, pour cela il faut indiquer le mot cle "
<< "\n # et_pas_recalcul_masse_a_la_transition_ a la fin de la declaration "
<< "\n # des parametres de controle du re-calcul de la matrice masse (voir les declarations"
<< "\n # associees a ce parametre (par defaut, a la transition il y a recalcul de la masse)) "
<< "\n # "
<< "\n # 5) il est possible d'indiquer un parametre de recalcul de la masse, different du "
<< "\n # cas de l'amortissement cinetique. Pour cela, en plus et a la suite du mot cle "
<< "\n # option_recalcul_mass= on peut indiquer opt_visqueux_recal_mass= suivi du cas "
<< "\n # qui sera utilise pour le calcul de la masse lorsque l'on sera en amortissement "
<< "\n # visqueux"
<< "\n # Comme pour le premier parametre de recalcul de masse, a la place d'une "
<< "\n # valeur numerique, on peut indiquer la presence d'une fonction nD a l'aide "
<< "\n # du mot cle nom_fct_nD_option_recalcul_mass_ suivi du nom d'une fonction nD"
<< "\n # (cf. les informations precedentes concernant option_recalcul_mass) "
<< "\n # "
<< "\n # Important: en fait le parametre opt_visqueux_recal_mass= est utilisable avec "
<< "\n # tous les types de relaxation, lorsqu'il est presant, suivant qu'il s'agit d'un "
<< "\n # d'un amortissement cinetique ou visqueux, ce sera le premier ou le second"
<< "\n # type de recalcul de masse qui sera utilise "
<< "\n # "
<< "\n # Exemple de declaration, : "
<< "\n # (ici apres la transition vers l'amortissement visqueux, option_recalcul_mass "
<< "\n # initialement a 0, passe à 3) "
<< "\n # "
mise à jour 7.010
2023-05-03 17:23:49 +02:00
<< "\n # typeCalRelaxation= 4 lambda= 0.6 type_calcul_mass= 2 option_recalcul_mass= 0 \\ "
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
<< "\n # opt_visqueux_recal_mass= 3 "
<< "\n # parametre_calcul_de_la_masse_ casMass_relax= 3 "
<< "\n # proportion_cinetique= 0.1 "
<< "\n # parametre_recalcul_de_la_masse_ ncycle_calcul= 100 et_pas_recalcul_masse_a_la_transition_"
mise à jour 7.010
2023-05-03 17:23:49 +02:00
<< "\n # parametre_calcul_de_la_viscosite_ type_calcul_visqu_critique= 2 \\ "
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
<< "\n # opt_cal_C_critique= 1 f_= 0.9 "
<< "\n # mode_debug_= 100 "
<< "\n # "
<< "\n # "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # une fois typeCalRelaxation et les parametres generaux definis on trouve sur les lignes qui suivent "
<< "\n # (une ligne par paquet de parametres) "
<< "\n # 1) eventuellement les parametres permettant de controler le calcul de la pseudo-masse "
<< "\n # 2) eventuellement les parametres permettant de controler le calcul de la matrice visqueuse "
<< "\n # 3) eventuellement les parametres de controle du re-calcul de la matrice masse "
<< "\n # 4) eventuellement le parametre de controle du mode debug "
<< "\n # "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # --- parametres de controle du calcul de la masse ----------- "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # ** PARAMETRES FACULTATIF ** "
<< "\n # on indique tout d'abord le mot cle : parametre_calcul_de_la_masse_ suivi des parametres associe "
<< "\n # "
<< "\n # =====> si type_calcul_mass= 1 : on a le fonctionnement suivant: "
<< "\n # la pseudo-masse elementaire est calculee avec la formule: "
<< "\n # mi=lambda*(ep/8)*(alpha*K+beta*mu+gamma*Isig/3+theta/2*Sig_mises "
<< "\n # Les parametres lambda, alpha, beta, gamma et theta permettent donc de ponderer le calcul de la masse "
<< "\n # ep: epaisseur, K module de compressibilite, mu: module de cisaillement, Isig trace de sigma, "
<< "\n # Sig_mises la contrainte de mises "
<< "\n # casMass_relax permet de choisir parmi differentes methodes de calcul de la matrice masse "
<< "\n # =1 on cumule la formule aux noeuds, en fonction de tous les elements qui entourent l'element"
<< "\n # =2 on prend le maximum de la formule aux noeuds, en fonction des elements l'entourant"
<< "\n # =3 on prend la moyenne de la formule aux noeuds, en fonction des elements l'entourant "
<< "\n # =4 on prend la moyenne et on divise par la section, comme dans la formule originale de Barnes "
<< "\n # =5 idem 1, et on divise par la section, comme dans la formule originale de Barnes "
<< "\n\n # l'algorithme a recours automatiquement a la relaxation cinetique. On peut donc "
<< "\n # aussi modifier les paramtetres de reglage de cette partie, sinon ce sont les parametres "
<< "\n # par defaut qui sont retenu (voir doc) "
<< "\n # Exemple de declaration: "
<< "\n ## parametre_calcul_de_la_masse_ alpha= 1. beta= 1. gamma= 1. theta= 1. casMass_relax= 1 "
<< "\n # "
<< "\n # tous les parametres sont facultatifs: les valeurs par defaut sont: "
<< "\n # alpha= 1 ; beta= 1. ; gamma= 1. ; theta= 1. ; casMass_relax= 3 ; "
<< "\n # "
<< "\n # =====> si type_calcul_mass= 2 : on a le fonctionnement suivant: "
<< "\n # la masse elementaire du noeud i est calculee a partir de la raideur reelle "
<< "\n # notation: ai c'est le a du noeud i, idem pour b "
<< "\n # m_(a_i) = lambda/2 * S(a_i) , a=1 a dim de l'espace "
<< "\n # casMass_relax= permet de choisir comment la matrice de masse est calculer. "
<< "\n # Sur l'ensemble de la matrice de raideur assemblee normalement "
<< "\n # casMass_relax= 0 : S(ai) = |K(ai,ai)| "
<< "\n # casMass_relax= 1 : S(bi) = sup_a ( |K(ai,ai)|), a=1 a dim , b=1 a dim "
<< "\n # pour chaque assemblage de matrice elementaire : "
<< "\n # casMass_relax= 2 : S(ai) = sum_{k,b} |K_loc(ai,bk)|, k=1,nbn; b=1,dim; (majoration du th de Gerschgorin) "
<< "\n # casMass_relax= 3 : S(ci) = sup_a (sum_{k,b} |K_loc(ai,bk)|), k=1,nbn; b=1,dim; c=1,dim; (majoration du th de Gerschgorin) "
<< "\n # casMass_relax= 4 : S(ai) = sup (2.*|K_loc(ai,ai)|, sum_{k,b} |K_loc(ai,bk)|, k=1,nbn; b=1,dim) "
<< "\n # casMass_relax= 5 : S(cj) = sup_a (sup (2.*|K_loc(ai,ai)|, sum_{k,b} |K_loc(ai,bk)|, k=1,nbn; b=1,dim)), c=1,dim;"
<< "\n # "
<< "\n # de nouveau sur l'ensemble de la matrice de raideur assemblee normalement K = K_globale "
<< "\n # casMass_relax= 6 : S(ai) = sum_{k,b} |K(ai,bk)|, k=1,nbn; b=1,dim; (th de Gerschgorin) "
<< "\n # casMass_relax= 7 : S(ci) = sup_a (sum_{k,b} |K(ai,bk)|), k=1,nbn; b=1,dim; c=1,dim; (th de Gerschgorin) "
<< "\n # casMass_relax= 8 : S(ai) = sup (2.*|K(ai,ai)|, sum_{k,b} |K(ai,bk)|, k=1,nbn; b=1,dim) "
<< "\n # casMass_relax= 9 : S(cj) = sup_a (sup (2.*|K(ai,ai)|, sum_{k,b} |K(ai,bk)|, k=1,nbn; b=1,dim)), c=1,dim;"
<< "\n # "
<< "\n # Exemple de declaration: "
<< "\n ## parametre_calcul_de_la_masse_ casMass_relax= 3 "
<< "\n # "
<< "\n # tous les parametres sont facultatifs: les valeurs par defaut sont: "
<< "\n # casMass_relax= 3 "
<< "\n # "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # --- parametres de controle du re-calcul de la masse ----------- "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # ** PARAMETRES FACULTATIF ** "
<< "\n # on indique tout d'abord le mot cle : parametre_recalcul_de_la_masse_ "
<< "\n # suivi des parametres associes eventuelles "
<< "\n # "
<< "\n # "
<< "\n # "
<< "\n # Attention: "
<< "\n # si option_recalcul_mass= -1 : la matrice masse est recalculee a chaque iteration "
<< "\n # si option_recalcul_mass= 0 : la matrice masse est calculée au debut du calcul de l'increment et ensuite elle reste fixe "
<< "\n # il n'y a pas de parametres associes supplementaire "
<< "\n # si option_recalcul_mass= 1 : apres un calcul au debut, mise a jour a chaque maxi de l'energie cinetique "
<< "\n # il n'y a pas de parametres associes supplementaire "
<< "\n # si option_recalcul_mass= 2 : idem le cas 1, mais on garde la valeur maxi de la raideur entre la nouvelle et l'ancienne"
<< "\n # il n'y a pas de parametres associes supplementaire "
<< "\n # si option_recalcul_mass= 3 : recalcul apres ncycles, valeur qui doit etre indiquee apres le mot cle ncycle_calcul= "
<< "\n # de la maniere suivante par exemple pour un calcul tous les 100 "
<< "\n ## parametre_recalcul_de_la_masse_ ncycle_calcul= 100 "
<< "\n # par defaut ncycle_calcul= 100 "
<< "\n # si option_recalcul_mass= 4 : recalcule de la matrice masse lorsque l'indicateur suivant "
<< "\n # \epsilon = MAX_{i=1}^{nbddl}{ ( 0.5*lambda * (Delta t)^2) (|Delta ddot X_i|) / (|Delta X_i|) } "
<< "\n # est superieur a 1*fac_epsilon "
<< "\n # par defaut fac_epsilon = 1, on peut le modifier a l'aide du parametre fac_epsilon= une valeur "
<< "\n # ( ce type de controle provient de l'amortissement visqueux critique) "
<< "\n # Exemple de declaration: "
<< "\n ## parametre_recalcul_de_la_masse_ fac_epsilon= 1. "
<< "\n # "
<< "\n # tous les parametres sont facultatifs: les valeurs par defaut sont: "
<< "\n # fac_epsilon= 1. ; "
<< "\n # "
<< "\n # Remarque: a la fin de ligne relative au mot cle parametre_recalcul_de_la_masse_ "
<< "\n # il est possible d'indiquer le mot cle et_pas_recalcul_masse_a_la_transition_ "
<< "\n # la presence de ce mot cle indique alors a l'algorithme de ne pas recalculer la masse "
<< "\n # lors de la transition cinetique - visqueux pour le cas: typeCalRelaxation= 4 "
<< "\n # pour les autres typeCalRelaxation la presence du mot cle "
<< "\n # et_recalcul_masse_a_la_transition_ n'a aucun effet "
<< "\n # le fonctionnement pas defaut est de recalculer la masse a la transition "
<< "\n # "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # --- parametres de pilotage de lambda ----------- "
<< "\n # ---------------------------------------------------------------------- "
<< "\n # ** PARAMETRES FACULTATIF ** "
<< "\n # trois types de pilotage sont disponibles au choix: "
<< "\n # "
<< "\n # 1) pilotage en fonction du temps physique: syntaxe "
<< "\n # lambdaAvecPonderation_temps_ suivi du nom d'une fonction 1D "
<< "\n # exemple de declaration: "
<< "\n # lambdaAvecPonderation_temps_ evolLambda # evolLambda est le nom d'une courbe 1D"
<< "\n # fonctionnement: a chaque increment de temps physique: "
<< "\n # lambda(t)= lambda(0) * evolLambda(t) "
<< "\n # avec lambda(0) le lambda fournit apres le mot clef: lambda="
<< "\n # "
<< "\n # 2) pilotage a l'aide d'une fonction nD: syntaxe "
<< "\n # lambdaAvecPonderation_Globale_ suivi du nom d'une fonction nD "
<< "\n # exemple de declaration: "
<< "\n # lambdaAvecPonderation_Globale_ fonct_cinetique # fonct_cinetique est le nom d'une courbe nD"
<< "\n # fonctionnement: a chaque iteration: "
<< "\n # lambda = lambda(0) * fonct_cinetique(grandeur(s) globale(s)) "
<< "\n # avec lambda(0) le lambda fournit apres le mot clef: lambda= "
<< "\n # "
<< "\n # 3) pilotage automatique: syntaxe "
<< "\n # pilotageAutolambda_ lambda_min= <une valeur> lambda_max= <une valeur> delta_lambda= <une valeur> "
<< "\n # exemple de declaration: "
<< "\n # pilotageAutolambda_ lambda_min= 0.6 lambda_max= 7 delta_lambda= 0.1 "
<< "\n # fonctionnement: on observe l'evolution des relaxations et lorsqu'il y a trop de relaxations "
<< "\n # lambda est augmente de delta_lambda a condition que la nouvelle valeur soit <= a lambda_max "
<< "\n # lambda_min ne sert pas pour l'instant "
<< "\n # Plus precisemment (version 6.799) soit la liste des iterations de relaxation "
<< "\n # le premier de la liste etant relatif a la derniere relaxation "
<< "\n # -> si le (3e - le premier) < 7 => augmentation de lambda "
<< "\n # -> si le (4e - le premier) < 10 => augmentation de lambda "
<< "\n # -> si le (5e - le premier) < 21 => augmentation de lambda "
<< "\n # "
<< "\n # ........................................................."
<< "\n # / parametre specifique au contact /"
<< "\n #........................................................."
<< "\n # ** PARAMETRES FACULTATIF ** "
<< "\n # on indique tout d'abord le mot cle : parametre_activation_du_contact_ suivi des parametres associes eventuelles "
<< "\n # type_activation_contact= 0 : la recherche de nouveaux contact est activee "
<< "\n # a la fin de chaque increment"
<< "\n # type_activation_contact= 1 : (valeur par defaut) la recherche de nouveaux contact est effectuee "
<< "\n # a la fin de chaque iteration "
<< "\n # type_activation_contact= 2 : la recherche de nouveaux contact est effectuee apres chaque amortissement cinetique "
<< "\n # et a la fin de chaque increment (s'il n'y a pas d'amortissement cinetique c'est "
<< "\n # equivalent au cas = 0 ) "
<< "\n # Exemple de declaration: "
<< "\n ## parametre_activation_du_contact_ type_activation_contact= 1 "
<< "\n # "
<< "\n # "
<< "\n # Dans certain cas, la construction de la matrice masse conduit a des 0, il existe un parametre "
<< "\n # pour gerer le remplacement ou non de 0 sur la diagonale "
<< "\n # "
<< "\n # exemple de declaration: "
<< "\n # choix_mini_masse_nul= 0 "
<< "\n # "
<< "\n # Signification du parametre: "
<< "\n # =0 -> pas de modification du terme nul "
<< "\n # =1 -> remplacement du terme nul par le maxi des valeurs de la matrice masse deja calculee "
<< "\n # =2 -> remplacement du terme nul par la moyenne des valeurs de la matrice masse deja calculee "
<< "\n # "
<< "\n # ................................................................"
<< "\n # / mise en place des conditions limites a chaque iteration /"
<< "\n #.................................................................."
<< "\n # mot cle : cL_a_chaque_iteration_ suivi de 0 ou 1 "
<< "\n # par defaut == 0, c-a-d que les CL sont mises en place "
<< "\n # au debut de chaque increment et ensuite reste fixe pendant les iterations "
<< "\n # == 1 : a chaque iteration, les CL sont de nouveaux imposes "
<< "\n # "
<< "\n # "
<< "\n # "
;
// appel de la classe mère
Algori::Info_com_parametres(entreePrinc);
sort << "\n" << endl;
};
//------- décomposition en 3 du calcul d'équilibre -------------
// a priori : InitAlgorithme et FinCalcul ne s'appellent qu'une fois,
// par contre : CalEquilibre peut s'appeler plusieurs fois, le résultat sera différent si entre deux calculs
// certaines variables ont-été changés
// de même : MiseAJourAlgo a pour objectif de mettre à jour (vérif de la cohérence) des variables internes
// entre deux appels de CalEquilibre
// initialisation
void AlgoriRelaxDyna::InitAlgorithme(ParaGlob * paraGlob,LesMaillages * lesMail,
LesReferences* lesRef,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD
,VariablesExporter* varExpor
,LesLoisDeComp* lesLoisDeComp,DiversStockage* diversStockage,
Charge* charge,LesCondLim* lesCondLim,LesContacts* lescontacts
,Resultats* resultats)
{ // INITIALISATION globale
tempsInitialisation.Mise_en_route_du_comptage(); // temps cpu
Transfert_ParaGlob_ALGO_GLOBAL_ACTUEL(RELAX_DYNA); // transfert info
mise à jour 7.010
2023-05-03 17:23:49 +02:00
#ifdef UTILISATION_MPI
// calcul de l'équilibrage initiale par le cpu 0
if (distribution_CPU_algo.Tableau_element_CPU_en_cours()->Taille() == 0 )
{distribution_CPU_algo.Calcul_Equilibrage_initiale(lesMail);
V 7.017 - intro para limite_inf_bIIb_ et limite_inf_Qeps_ dans les potentiels de grenoble - cor bug diff finie pour Qeps petit sur loi grenoble - modif prise en compte de la regularisation loi grenoble (passage systématique en différence finie pour Qeps très faible)
2023-07-06 08:42:57 +02:00
temps_transfert_court.Mise_en_route_du_comptage(); // comptage cpu
mise à jour 7.010
2023-05-03 17:23:49 +02:00
distribution_CPU_algo.Passage_Equilibrage_aux_CPU();
paraGlob->Init_tableau (distribution_CPU_algo.Tableau_element_CPU_en_cours());
};
#endif
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// avant toute chose, au cas où l'algo interviendrait après un autre algo
// on inactive tous les ddl existants
lesMail->Inactive_ddl();
// on regarde s'il s'agit d'un pb purement non méca, dans ce cas il faut cependant initialiser les positions
// à t et tdt pour le calcul de la métrique associée
{ const list <EnumElemTypeProblem >& type_pb = lesMail->Types_de_problemes();
bool purement_non_meca = true;
list <EnumElemTypeProblem >::const_iterator il,ilfin=type_pb.end();
for (il=type_pb.begin();il != ilfin; il++)
{switch (*il)
{ case MECA_SOLIDE_DEFORMABLE: case MECA_SOLIDE_INDEFORMABLE: case MECA_FLUIDE:
purement_non_meca=false; break;
default: break; // sinon on ne fait rien
};
};
if (purement_non_meca) // si pas de méca, on initialise les coordonnées à t et tdt avec celles de 0
{lesMail->Init_Xi_t_et_tdt_de_0();}
else // sinon a minima on active X1
{ lesMail->Active_un_type_ddl_particulier(X1);
};
};
// cas du chargement, on verifie egalement la bonne adequation des references
charge->Initialise(lesMail,lesRef,pa,*lesCourbes1D,*lesFonctionsnD);
// on indique que l'on ne souhaite pas le temps fin stricte
// (sinon erreur non gérée après un changement de delta t), que l'on suppose négligeable
// après plusieurs incréments
// modif 21 juin 2017 : non on veut une gestion stricte (je ne sais pas pourquoi il y avait
// une telle gestion, sans doute une suite à un calcul explicite classique
charge->Change_temps_fin_non_stricte(0); //(1);
// dans le cas où l'on calcul des contraintes et/ou déformation et/ou un estimateur d'erreur
// à chaque incrément, initialisation
tenuXVG.Change_taille(3);tenuXVG(1)=X1;tenuXVG(2)=V1;tenuXVG(3)=GAMMA1;
prepa_avec_remont = Algori::InitRemont(lesMail,lesRef,diversStockage,charge,lesCondLim,lescontacts,resultats);
if ( prepa_avec_remont)// remise enservice des ddl du pab
{lesMail->Inactive_ddl(); lesMail->Active_un_type_ddl_particulier(X1);};
// dans le cas de l'algorithme de relaxation dynamique, le temps est normalement imposé à 1., donc n'intervient pas
// en fait dans l'algo d'avancement temporel !!
// Cependant comme le temps intervient également dans le chargement on le laisse à la valeur indiqué par l'utilisateur
this->Gestion_pas_de_temps(true,lesMail,1); // 1 signifie qu'il y a initialisation
// 00 --<ARD>-- on crée les ddl d'accélération et de vitesse non actif mais libres
// car les forces intérieures et extérieures sont les entitées duales
// des déplacements, qui sont donc les seules grandeurs actives à ce stade
lesMail->Plus_Les_ddl_Vitesse( HSLIBRE);
lesMail->Plus_Les_ddl_Acceleration( HSLIBRE);
// on défini globalement que l'on a une combinaison des ddl X V GAMMA en même temps
cas_combi_ddl=1;
// mise en place éventuelle du bulk viscosity
lesMail->Init_bulk_viscosity(pa.BulkViscosity(),pa.CoefsBulk());
// mise a zero de tous les ddl et creation des tableaux a t+dt
// les ddl de position ne sont pas mis a zero ! ils sont initialise
// a la position courante
lesMail->ZeroDdl(true);
// on vérifie que les noeuds sont bien attachés à un élément sinon on met un warning si niveau > 2
if (ParaGlob::NiveauImpression() > 2)
lesMail->AffichageNoeudNonReferencer();
// init des ddl avec les conditions initiales
// les conditions limites initiales de vitesse sont prise en compte
// de manière identiques à des ddl quelconques, ce ne sont pas des ddl fixé !!
// par contre l'accélération initiale est déterminée à l'aide de l'équilibre initiale
// donc les conditions d'accélérations initiale seront écrasées
lesCondLim->Initial(lesMail,lesRef,lesCourbes1D,lesFonctionsnD,true,cas_combi_ddl);
// mise à jour des différents pointeur d'assemblage et activation des ddl
// a) pour les déplacements qui sont à ce stade les seuls grandeurs actives
// on définit un nouveau cas d'assemblage pour les Xi
// à travers la définition d'une instance de la classe assemblage
if (Ass1_ == NULL) Ass1_ = new Assemblage(lesMail->InitNouveauCasAssemblage(1));
Assemblage& Ass1 = *Ass1_;
lesMail->MiseAJourPointeurAssemblage(Ass1.Nb_cas_assemb());// mise a jour des pointeurs d'assemblage
int nbddl_X = lesMail->NbTotalDdlActifs(X1); // nb total de ddl de déplacement
// qui est le même pour les accélérations et les vitesses
// b) maintenant le cas des vitesses qui doivent donc être activées
// on définit un nouveau cas d'assemblage pour les Vi
if (Ass2_ == NULL) Ass2_ = new Assemblage(lesMail->InitNouveauCasAssemblage(1));
Assemblage& Ass2 = *Ass2_;
lesMail->Inactive_un_type_ddl_particulier(X1); // on inactive les Xi
lesMail->Active_un_type_ddl_particulier(V1); // on active les Vi
lesMail->MiseAJourPointeurAssemblage(Ass2.Nb_cas_assemb()); // mise a jour des pointeurs d'assemblage
// c) idem pour les accélérations
// on définit le numéro de second membre en cours
// on définit un nouveau cas d'assemblage pour les pour GAMMAi
if (Ass3_ == NULL) Ass3_ = new Assemblage(lesMail->InitNouveauCasAssemblage(1));
Assemblage& Ass3 = *Ass3_;
lesMail->Inactive_un_type_ddl_particulier(V1); // on inactive les Vi
lesMail->Active_un_type_ddl_particulier(GAMMA1); // on active les GAMMAi
lesMail->MiseAJourPointeurAssemblage(Ass3.Nb_cas_assemb()); // mise a jour des pointeurs d'assemblage
// d) activation de tous les ddl, maintenant ils peuvent être les 3 actifs
lesMail->Active_un_type_ddl_particulier(X1);
lesMail->Active_un_type_ddl_particulier(V1);
// en fait ces trois pointeurs d'assemblage ne sont utils que pour la mise en place des conditions
// limites
// mise à jour du nombre de cas d'assemblage pour les conditions limites
// c-a-d le nombre maxi possible (intégrant les autres pb qui sont résolu en // éventuellement)
lesCondLim->InitNombreCasAssemblage(lesMail->Nb_total_en_cours_de_cas_Assemblage());
// définition d'un tableau globalisant les numéros d'assemblage de X V gamma
t_assemb.Change_taille(3);
t_assemb(1)=Ass1.Nb_cas_assemb();t_assemb(2)=Ass2.Nb_cas_assemb();t_assemb(3)=Ass3.Nb_cas_assemb();
// récupération des tableaux d'indices généraux des ddl bloqués, y compris les ddls associés
icas = 1; // pour indiquer au module Tableau_indice que l'on travaille avec l'association X V GAMMA
li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb,lesRef,charge->Temps_courant(),icas);
// on définit quatre tableaux qui serviront à stocker transitoirement les X V GAMMA correspondant au ddl imposés
int ttsi = li_gene_asso.size();
X_Bl.Change_taille(ttsi),V_Bl.Change_taille(ttsi),G_Bl.Change_taille(ttsi);
// def vecteurs globaux
vglobin.Change_taille(nbddl_X); // puissance interne : pour ddl accélération
vglobex.Change_taille(nbddl_X); // puissance externe
vglobaal.Change_taille(nbddl_X,0.); // puissance totale
// même si le contact n'est pas encore actif, il faut prévoir qu'il le deviendra peut-être !
if (lesMail->NbEsclave() != 0)
vcontact.Change_taille(nbddl_X); // puissance de contact
// 04 --<ARD>-- 6 vecteur pour une manipulation globale des positions vitesses et accélérations
// vecteur qui globalise toutes les positions de l'ensemble des noeuds
// dans le cas où on initialise les ddl_tdt(n) avec les ddl_t(n)+delta_ddl(n-1), def de grandeurs
X_t.Change_taille(nbddl_X);X_tdt.Change_taille(nbddl_X); delta_X.Change_taille(nbddl_X);
var_delta_X.Change_taille(nbddl_X);
// vecteur qui globalise toutes les vitesses de l'ensemble des noeuds
vitesse_t.Change_taille(nbddl_X);vitesse_tdt.Change_taille(nbddl_X);
// vecteur qui globalise toutes les accélérations
acceleration_t.Change_taille(nbddl_X);acceleration_tdt.Change_taille(nbddl_X) ;
// calcul des énergies
if (pa.Amort_visco_artificielle()) // dans le cas d'un amortissement artificiel
forces_vis_num.Change_taille(nbddl_X,0.);
E_cin_tdt = 0.; E_int_t = 0.; E_int_tdt = 0.; // init des différentes énergies
E_ext_t = 0.; E_ext_tdt = 0.; bilan_E = 0.; // " et du bilan
save_X_t.Change_taille(nbddl_X); // vecteur de sauvegarde
F_int_t.Change_taille(nbddl_X); F_ext_t.Change_taille(nbddl_X); // forces généralisées int et ext au pas précédent
F_int_tdt.Change_taille(nbddl_X); F_ext_tdt.Change_taille(nbddl_X); // forces généralisées int et ext au pas actuel
residu_final.Change_taille(nbddl_X); // pour la sauvegarde du résidu pour le post-traitement
// initialisation du compteur d'increments de charge
// au premier passage il y a un traitement spécial: pas d'incrémentation du temps, pas de sauvegarde
// in s'agit d'un incrément de mise en place, l'équilibre se faisant à la fin
// ici en relaxation dynamique on démarre à -1 car pour -1 et 0 on ne fait rien, on initialise (voir le calcul)
compteur_demarrage = -1; // pas vraiment utile, car sera redéfini dans CalEquilibre
icharge = 0; // par défaut
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// definition des elements de frontiere, ces elements sont utilises pour le contact
lesMail->CreeElemFront();
// calcul éventuel des normales aux noeuds -> init des normales pour t=0
lesMail->InitNormaleAuxNoeuds(); //utilisé pour la stabilisation des membranes par ex
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// --- init du contact ---
// doit-être avant la lecture d'un restart, car il y a une initialisation de conteneurs qui est faites
// qui ensuite est utilisée en restart
// par exemple il faut initialiser les frontières et la répartition esclave et maître
// pour préparer la lecture de restart éventuel
if (lesMail->NbEsclave() != 0)
{ // definition des elements de frontiere, ces elements sont utilises pour le contact
lesMail->Mise_a_jour_boite_encombrement_elem_front(TEMPS_t); //
// initialisation des zones de contacts éventuelles
lescontacts->Init_contact(*lesMail,*lesRef,lesFonctionsnD);
// verification qu'il n'y a pas de contact avant le premier increment de charge
lescontacts->Verification();
// definition des elements de contact eventuels
// et imposition éventuel de certaines des conditions de contact (dépend du modèle de contact)
lescontacts->DefElemCont(0.); // au début le déplacement des noeuds est nul
if (pa.ContactType() == 4) // cas particulier du type 4 de contact où on utilise les forces internes
{// def d'un type générique, utilisé pour le transfert des forces internes, vers les conteneurs noeuds
Coordonnee coor(ParaGlob::Dimension()); // un type coordonnee typique
Grandeur_coordonnee gt(coor); // une grandeur typique de type Grandeur_coordonnee
// def d'un type quelconque représentatif pour un vecteur force à chaque noeud
TypeQuelconque typQ_gene_int(FORCE_GENE_INT,X1,gt);
lesMail->AjoutConteneurAuNoeud(typQ_gene_int);
};
};
//--cas de restart et/ou de sauvegarde------------
// tout d'abord récup du restart si nécessaire
// dans le cas ou un incrément différent de 0 est demandé -> seconde lecture à l'incrément
bool brestart=false; // booleen qui indique si l'on est en restart ou pas
if (this->Num_restart() != 0)
{ int cas = 2;
// ouverture de base info
entreePrinc->Ouverture_base_info("lecture");
this->Lecture_base_info(cas ,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage
,charge,lesCondLim,lescontacts,resultats,(this->Num_restart()));
icharge = this->Num_restart();//+1;
// récup du pas de temps, proposé par l'utilisateur, initialisation et vérif / pas critique
this->Gestion_pas_de_temps(true,lesMail,2); // 2 signifie cas courant
brestart = true;
// on oblige les ddls Vi GAMMAi a avoir le même statut que celui des Xi
// comme les conditions limites cinématiques peuvent être différentes en restart
// par rapport à celles sauvegardées, on commence par libérer toutes les CL imposées éventuelles
lesMail->Libere_Ddl_representatifs_des_physiques(LIBRE);
lesMail->ChangeStatut(cas_combi_ddl,LIBRE);
// dans le cas d'un calcul axisymétrique on bloque le ddl 3
if (ParaGlob::AxiSymetrie())
mise à jour 7.010
2023-05-03 17:23:49 +02:00
lesMail->Inactive_un_ddl_particulier(X3);
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// on valide l'activité des conditions limites et condition linéaires, pour le temps initial
// en conformité avec les conditions lues (qui peuvent éventuellement changé / aux calcul qui a donné le .BI)
////--debug
//cout << "\n debug1: AlgoriRelaxDyna::InitAlgorithme(..";
//lesMail->Noeud_LesMaille(1,197).Affiche();cout << endl;
////--fin debug
li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb,lesRef,charge->Temps_courant(),icas);
int ttsi = li_gene_asso.size();
X_Bl.Change_taille(ttsi);V_Bl.Change_taille(ttsi);G_Bl.Change_taille(ttsi);
lesCondLim->Validation_blocage (lesRef,charge->Temps_courant());
// mise à jour pour le contact s'il y du contact présumé
if (pa.ContactType())
lesMail->Mise_a_jour_boite_encombrement_elem_front(TEMPS_t);
};
// vérif de cohérence pour le contact
if ((pa.ContactType()) && (lesMail->NbEsclave() == 0)) // là pb
{cout << "\n *** erreur: il n'y a pas de maillage disponible pour le contact "
<< " la definition d'un type contact possible est donc incoherente "
<< " revoir la mise en donnees !! "<< flush;
Sortie(1);
};
// on regarde s'il y a besoin de sauvegarde
if (this->Active_sauvegarde() && (ParaGlob::param->TypeCalcul_maitre() == this->typeCalcul) )
{ // si le fichier base_info n'est pas en service on l'ouvre
entreePrinc->Ouverture_base_info("ecriture");
// dans le cas ou se n'est pas un restart on sauvegarde l'incrément actuel
// c'est-à-dire le premier incrément
// après s'être positionné au début du fichier
if (this->Num_restart() == 0)
{ (entreePrinc->Sort_BI())->seekp(0);
int cas = 1;
paraGlob->Ecriture_base_info(*(entreePrinc->Sort_BI()),cas);
this->Ecriture_base_info
(cas,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage
,charge,lesCondLim,lescontacts,resultats,OrdreVisu::INCRE_0);
}
else
{ // sinon on se place dans le fichier à la position du restart
// debut_increment a été définit dans algori (classe mère)
(entreePrinc->Sort_BI())->seekp(debut_increment);
}
};
// ajout d'un conteneur pour les coordonnées à l'itération 0
{Coordonnee coor(ParaGlob::Dimension()); // un type coordonnee typique
Grandeur_coordonnee gt(coor); // une grandeur typique de type Grandeur_coordonnee
// def d'un type quelconque représentatif à chaque noeud
TypeQuelconque typQ_gene_int(XI_ITER_0,X1,gt);
lesMail->AjoutConteneurAuNoeud(typQ_gene_int);
};
// init pour var glob
Transfert_ParaGlob_COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL(icharge);
//--fin cas de restart et/ou de sauvegarde--------
// on signale que l'on utilise un comportement matériel normal
lesLoisDeComp->Loi_simplifie(false);
// choix de la matrice de masse, qui est en fait celle qui correspond au ddl Xi
// ici le numéro d'assemblage est celui de X car on projette bien sur des vitesses virtuelles c-a-d ddl X*.
if (ParaGlob::NiveauImpression() > 2)
cout << "\n matrice masse principale: ";
mat_masse = Choix_matrice_masse(nbddl_X,mat_masse,lesMail,lesRef
,Ass1.Nb_cas_assemb(),lescontacts,lesCondLim);
if (ParaGlob::NiveauImpression() > 2)
cout << "\n matrice masse secondaire: ";
mat_masse_sauve = Choix_matrice_masse(nbddl_X,mat_masse_sauve,lesMail,lesRef
,Ass1.Nb_cas_assemb(),lescontacts,lesCondLim);
// choix de la résolution
if (mat_masse->Type_matrice() == DIAGONALE)
// dans le cas d'une matrice diagonale on force la résolution directe quelque soit l'entrée
mat_masse->Change_Choix_resolution(DIRECT_DIAGONAL,pa.Type_preconditionnement());
else
mat_masse->Change_Choix_resolution(pa.Type_resolution(),pa.Type_preconditionnement());
// initialisation par exemple au niveau des maillages (éléments) pour les calculs futurs de la matrice masse
InitCalculMatriceMasse(lesMail,*mat_masse,Ass1,*mat_masse_sauve,lesFonctionsnD);
(*mat_masse_sauve) = (*mat_masse); // sauvegarde avant la mise en place des conditions limites
if ((casMass_relax > 5) && (casMass_relax < 10)) // cas où on utilise la matrice de raideur réelle
{ // choix des matrices de raideur de sustitution éventuelles : par défaut matglob = tab_mato(1)
Tableau < Mat_abstraite*> tab_mato(1);
// choix de la matrice de raideur du système linéaire (calcul éventuellement d'une largeur de bande ad oc)
Choix_matriciel(nbddl_X,tab_mato,lesMail,lesRef,Ass1.Nb_cas_assemb(),lesCondLim);
matglob = tab_mato(1);
// choix de la résolution pour la première matrice
matglob->Change_Choix_resolution(pa.Type_resolution(),pa.Type_preconditionnement());
};
// on vérifie que le ddl statique a bien été définit
// on définit le Ddl_enum_etendu correspondant à la masse au noeud pour la méthode
if (type_calcul_mass == 1)
{if (masseRelax_1.Nom_vide() )
masseRelax_1 = (Ddl_enum_etendu(string("masse_relax_dyn")));
// introduction des conteneurs de masse particulier à la méthode, sert de manière générale
// également pour le post-traitement
{ List_io < Ddl_enum_etendu > lienu; // def uniquement pour le passage de param
lienu.push_back(masseRelax_1);
lesMail->AjoutConteneurAuNoeud(lienu);
};
}
else if (type_calcul_mass == 2)
{ if ( masseRelax_2.EnuTypeQuelconque().EnumTQ() == RIEN_TYPEQUELCONQUE )
{ int dim = ParaGlob::Dimension();
// def des grandeurs courantes de type coordonnee
Coordonnee coor(dim); // un type coordonnee typique
// maintenant on définit une grandeur typique de type Grandeur_coordonnee
Grandeur_coordonnee gt(coor);
// def d'un type quelconque représentatif pour un vecteur pseudomasse
TypeQuelconque typQ(MASSE_RELAX_DYN,X1,gt);
// on ajoute le conteneur au noeud
lesMail->AjoutConteneurAuNoeud(typQ);
// on met à jout le conteneur statique
masseRelax_2.ChangeGrandeur(typQ);
};
}
else
{ cout << "\n *** erreur :type_calcul_mass est different de 1 ou 2 ";
Sortie(1);
};
// dans le cas où l'on utilise de l'amortissement numérique visqueux classique
// en fait ici cela correspond à : forces_vis_num = c [M] V(n)
// pour le second membre il est également nécessaire de tenir compte du terme -[C] V(n)
// il nous faut donc la matrice de viscosité
// forces_vis_num: forces visqueuses d'origines numériques
if ((typeCalRelaxation == 2) || pa.Amort_visco_artificielle() || (typeCalRelaxation == 3)
|| (typeCalRelaxation == 4))
{ bool initial = true; // def de la matrice (place et valeurs)
mat_C_pt = Cal_mat_visqueux_num_expli(*mat_masse,mat_C_pt,delta_X,initial,vitesse_tdt);
forces_vis_num.Change_taille(nbddl_X);
if (Arret_A_Equilibre_Statique()) // si on veut un équilibre statique, on sauvegarde les forces statiques
{ if (vglob_stat != NULL)
{vglob_stat->Change_taille(vglobaal.Taille());}
else
{vglob_stat = new Vecteur(vglobaal.Taille());}
};
};
// mise en place des conditions limites
// ---- initialisation des sauvegardes sur matrice et second membre
// ce qui ne correspond à rien ici normalement
lesCondLim->InitSauve(Ass3.Nb_cas_assemb());
//
lesCondLim->ImposeConLimtdt(lesMail,lesRef,*mat_masse,vglobaal,Ass3.Nb_cas_assemb()
,cas_combi_ddl,vglob_stat);
// puis on prépare (si besoin est en fonction du type de matrice) la résolution
if (!(lesCondLim->ExisteCondiLimite()))
mat_masse->Preparation_resol();
type_incre = OrdreVisu::PREMIER_INCRE; // pour la visualisation au fil du calcul
// on indique (si ça n'a pas été déjà fait) que l'on va utiliser la relaxation cinétique
if (typeCalRelaxation != 2)
{amortissement_cinetique=true; // paramètre définit dans Algori
Algori::InitialiseAmortissementCinetique(); // initialisation des compteurs pour l'amortissement
}
else {amortissement_cinetique=false; };
// dans le cas où les fonctions niveauF_grandeurGlobale ou/et niveauF_temps sont non nulles
// idem pour niveauLambda_grandeurGlobale ou/et niveauLambda_temps
// on finit l'affectation
if (niveauLambda_temps != NULL)
niveauLambda_temps->Affectation_fonctions(*lesCourbes1D);
if (niveauLambda_grandeurGlobale != NULL)
{niveauLambda_grandeurGlobale->Affectation_fonctions(*lesFonctionsnD);
// on vérifie que c'est bien une fonction scalaire
int nbcomp = niveauLambda_grandeurGlobale->C_proport()->NbComposante();
if (nbcomp != 1)
{cout << "\n erreur avec la fct " << niveauLambda_grandeurGlobale->Nom_fonction()
<< "definissant l'evolution de lambda, ce n'est pas une fct scalaire"
<< " elle retourne "<<nbcomp << "composantes ! verifier la mise en donnees ";
Sortie(1);
};
};
if (niveauF_temps != NULL)
niveauF_temps->Affectation_fonctions(*lesCourbes1D);
if (niveauF_grandeurGlobale != NULL)
{niveauF_grandeurGlobale->Affectation_fonctions(*lesFonctionsnD);
// on vérifie que c'est bien une fonction scalaire
int nbcomp = niveauF_grandeurGlobale->C_proport()->NbComposante();
if (nbcomp != 1)
{cout << "\n erreur avec la fct " << niveauF_grandeurGlobale->Nom_fonction()
<< "definissant la bascule cinetique-visqueux, ce n'est pas une fct scalaire"
<< " elle retourne "<<nbcomp << "composantes ! verifier la mise en donnees ";
Sortie(1);
};
};
// mise à jour variable globale: ici en initialisation
switch (typeCalRelaxation)
{ case 1: visqueux_activer = false; break;
case 2: visqueux_activer = true; break;
case 4: visqueux_activer = false; break; // au début en cinétique
default:
cout <<"\n typeCalRelaxation = " << typeCalRelaxation
<< " cas non pris en compte dans AlgoriRelaxDyna::Lecture_Base_info_Parametre( ";
Sortie(1);
break;
};
Transfert_ParaGlob_AMOR_CINET_VISQUEUX();
// ---- on renseigne les pointeurs de fonctions nD éventuels
// on récupère éventuellement les fonctions qui correspondent à fct_nD_option_recalcul_mass
int nb_opt_cal_masse = nom_fct_nD_option_recalcul_mass.Taille();
for (int i=1;i<= nb_opt_cal_masse;i++)
{if (nom_fct_nD_option_recalcul_mass(i).length())
{fct_nD_option_recalcul_mass(i) = lesFonctionsnD->Trouve(nom_fct_nD_option_recalcul_mass(i));
if (fct_nD_option_recalcul_mass(i)->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales
{// on vérifie qu'en retour on a un scalaire
fct_nD_option_recalcul_mass(i)->Valeur_pour_variables_globales(); // pour simplifier
if (fct_nD_option_recalcul_mass(i)->NbComposante() != 1)
{cout << "\n *** erreur parametre algorithme: la fonction nD "
<< nom_fct_nD_option_recalcul_mass(i)
<< " qui permet de calculer la valeur du parametre recalcul_mass("<<i<<") "
<< " ne retourne pas un scalaire unique !! a priori ce n'est pas correct ";
Sortie(1);
};
}
else
{cout << "\n *** erreur parametre algorithme: la fonction nD "
<< nom_fct_nD_option_recalcul_mass(i)
<< " qui permet de calculer la valeur du parametre recalcul_mass("<<i<<") "
<< " utilise des variables autres que globales, ce n'est pas possible ! ";
Sortie(1);
};
};
};
// mise à jour au cas où
Algori::MiseAJourAlgoMere(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp,diversStockage
,charge,lesCondLim,lescontacts,resultats);
tempsInitialisation.Arret_du_comptage(); // temps cpu
};
// mise à jour
void AlgoriRelaxDyna::MiseAJourAlgo
(ParaGlob * paraGlob,LesMaillages * lesMail
,LesReferences* lesRef,LesCourbes1D* lesCourbes1D
,LesFonctions_nD* lesFonctionsnD
,VariablesExporter* varExpor
,LesLoisDeComp* lesLoisDeComp,DiversStockage* diversStockage
,Charge* charge,LesCondLim* lesCondLim,LesContacts* lescontacts
,Resultats* resultats
)
{ // INITIALISATION globale
tempsMiseAjourAlgo.Mise_en_route_du_comptage(); // temps cpu
Transfert_ParaGlob_ALGO_GLOBAL_ACTUEL(RELAX_DYNA); // transfert info
// activation des ddl
lesMail->Inactive_ddl(); // on commence par inactiver tous les ddl
lesMail->Active_un_type_ddl_particulier(X1); // puis on active les ddl qu'ils faut ici
lesMail->Active_un_type_ddl_particulier(V1); // on active les Vi
lesMail->Active_un_type_ddl_particulier(GAMMA1); // on active les GAMMAi
// ici en relaxation dynamique on démarre à -1 car pour -1 et 0 on ne fait rien, on initialise (voir le calcul)
compteur_demarrage = -1; // mais en fait pas vraiment utile, car sera redéfini dans CalEquilibre
// mise à jour au cas où
Algori::MiseAJourAlgoMere(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,diversStockage,charge,lesCondLim,lescontacts,resultats);
// init indicateur
Transfert_ParaGlob_AMOR_CINET_VISQUEUX();
tempsMiseAjourAlgo.Arret_du_comptage(); // temps cpu
};
// calcul
// si tb_combiner est non null -> un tableau de 2 fonctions
// - la première fct dit si on doit valider ou non le calcul à convergence ok,
// - la seconde dit si on doit sortir de la boucle ou non à convergence ok
//
// si la validation est effectuée, la sauvegarde pour le post-traitement est également effectuée
// en fonction de la demande de sauvegard,
// sinon pas de sauvegarde pour le post-traitement à moins que l'on a demandé un mode debug
// qui lui fonctionne indépendamment
void AlgoriRelaxDyna::CalEquilibre(ParaGlob * paraGlob,LesMaillages * lesMail
,LesReferences* lesRef,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD
,VariablesExporter* varExpor
,LesLoisDeComp* lesLoisDeComp,DiversStockage* diversStockage
,Charge* charge,LesCondLim* lesCondLim,LesContacts* lescontacts
,Resultats* resultats,Tableau < Fonction_nD* > * tb_combiner)
{
tempsCalEquilibre.Mise_en_route_du_comptage(); // temps cpu
Transfert_ParaGlob_ALGO_GLOBAL_ACTUEL(RELAX_DYNA); // transfert info
// on vérifie que le pas de temps est correcte: on propose le deltat actuel qui sera modifié s'il n'est pas bon
double deltat_actuel = pa.Deltat();
pa.Modif_Detat_dans_borne(deltat_actuel);
// récup des entités
Assemblage& Ass1 = *Ass1_;
// Assemblage& Ass2 = *Ass2_;
Assemblage& Ass3 = *Ass3_;
// préparation pour les aspects validation du calcul et sortie contrôlée des incréments
int validation_calcul = 1; // init : par défaut la validation est effective si le calcul converge
int sortie_boucle_controler_par_fctnD = 0; // init : par défaut la sortie n'est pas contrôler
// boucle sur les increments de charge qui gère également les incréments de temps
// tant que la fin du chargement n'est pas atteinte
// dans le cas du premier chargement on calcul de toute manière, ce qui permet
// de calculer meme si l'utilisateur indique un increment de charge supérieur
// au temps final
// init de la valeur de lambda
lambda = Calcul_Lambda();
// on impose une accélération nulle compte tenu du fait qu'elle est ici purement numérique
acceleration_tdt.Zero();
vitesse_tdt.Zero(); // idem vitesses (pourrait poser des pb si l'on veut utiliser des vrais vitesses aux noeuds ?)
// ********* à améliorer
// un vecteur de travail qui devra être mieux dimensionné par la suite
Vecteur V_ext(F_int_tdt);
double max_delta_X=0.; // le maxi du delta X
double max_var_delta_X=0.; // idem d'une itération à l'autre
// def d'un type générique, utilisé pour le transfert des forces internes, vers les conteneurs noeuds
int dim = ParaGlob::Dimension();if (ParaGlob::AxiSymetrie()) dim--;
Coordonnee coor(dim); // un type coordonnee typique
Grandeur_coordonnee gt(coor); // une grandeur typique de type Grandeur_coordonnee
// def d'un type quelconque représentatif pour un vecteur force à chaque noeud
TypeQuelconque typQ_gene_int(FORCE_GENE_INT,X1,gt);
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// on définit un type générique qui sert pour passer aux noeuds les positions à l'itération 0
TypeQuelconque typQ_XI_ITER_0(XI_ITER_0,X1,gt);
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
bool arret=false; // booleen pour arrêter indépendamment de la charge
mise à jour 7.010
2023-05-03 17:23:49 +02:00
bool arret_pilotage=false; // pour arrêt du calcul au niveau du pilotage
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
bool premier_calcul = true; // utilisé pour l'initialisation de l'incrément avec le pas précédent
int indicCycleContact = 0; // def d'un indicateur donnant la situation dans le cycle de contact
// un booléen pour uniquement gérer le fait que dans la boucle globale on fait le test après le test du while
// d'où pour le temps final, on peut sortir de la boucle globale sans avoir convergé, ce qui n'est pas bon
bool pas_de_convergence_pour_l_instant=1;
// ici en relaxation dynamique on démarre à -1 car pour -1 et 0 on ne fait rien, on initialise (voir le calcul)
compteur_demarrage=-1; // donc compteur pour les premiers incréments
icharge++; // on incrémente le chargement -> donne le num d'incr du prochain incr chargé
while ( ((!charge->Fin(icharge,((!pas_de_convergence_pour_l_instant)&&(compteur_demarrage < 1)))) // on n'affiche la sin
// que si on a eu convergence
|| pas_de_convergence_pour_l_instant ||(compteur_demarrage <= 1)
)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
&& (charge->Fin(icharge,true)!=1) // si on a dépassé le temps fin on s'arrête
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// si on a dépassé le nombre d'incrément permis on s'arrête dans tous les cas
&& ((compteur_demarrage < 1)? (charge->Fin(icharge,false)!=2) : (charge->Fin(icharge,true)!=2))
mise à jour 7.010
2023-05-03 17:23:49 +02:00
&& (charge->Fin(icharge,true)!=3) // idem si on a dépassé le nombre d'essai d'incrément permis
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// 1er appel avec true: pour affichage et second avec false car c'est déjà affiché
&& (!pa.EtatSortieEquilibreGlobal())
)
{ double maxPuissExt; // maxi de la puissance des efforts externes
double maxPuissInt; // maxi de la puissance des efforts internes
double maxReaction; // maxi des reactions
int inReaction = 0; // pointeur d'assemblage pour le maxi de reaction
Transfert_ParaGlob_NORME_CONVERGENCE(ConstMath::grand);// on met une norme grande par défaut au début
lesLoisDeComp->MiseAJour_umat_nbincr(icharge); // init pour les lois Umat éventuelles
// --- en fait on doit faire 2 passages pour initialiser l'algo:
// * premier passage calcul de l'accélération à t=0
// --> on utilise un deltat = 0 et t = 0,
// * second passage --> on met un bon pas de temps, mais toujours pas
// d'accroissement du chargement ni de calcul de la relaxation dynamique
// (car il n'y a pas de référence), par contre on va
// calculer toutes les infos : accelération et vitesse à t = 0,
// * c'est seulement à partir de compteur_demarrage > 0 que
// l'on commence à calculer pour un t non nul et donc un accroissement
// du chargement non nul (mais il pouvait y avoir précédemment un chargement à t=0)
// bilan: on traite différamment selon que c'est le premier passage le second et les autres
bool aff_incr = true; // par défaut, est ensuite renseigné dans le switch qui suit
bool change_statut = false; // init des changements de statut
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// mise à jour du calcul éventuel des normales aux noeuds -> mise à jour des normales à t
// mais ici, on calcule les normales à tdt, et on transfert à t
// comme on est au début de l'incrément, la géométrie à tdt est identique à celle à t
// sauf "au premier incrément", si l'algo est un sous algo d'un algo combiné
// et que l'on suit un précédent algo sur un même pas de temps
// qui a aboutit à une géométrie à tdt différente de celle de t
// du coup cela permet d'utiliser la nouvelle géométrie pour ce premier incrément
lesMail->MiseAjourNormaleAuxNoeuds_de_tdt_vers_T();
// passage aux noeuds des vecteurs globaux: F_INT, F_EXT
Algori::Passage_aux_noeuds_F_int_t_et_F_ext_t(lesMail);
// qui, pour provoquer un arrêt de la boucle sur les incrément
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
switch (compteur_demarrage)
{case -1:
{// gestion du pas de temps, mis transitoirement à 0
this->Gestion_pas_de_temps(false,lesMail,0); // 0 signifie cas du passage à vide
break;
}
case 0:
{// on remet les choses dans l'ordre et on continue sur le cas courant
// mais c'est toujours une passe d'initialisation
this->Gestion_pas_de_temps(false,lesMail,1); // 1 signifie premier pas
break;
}
default: // cas normal ou compteur_demarrage > 0
{// gestion du pas de temps
this->Gestion_pas_de_temps(false,lesMail,2); // 2 signifie cas courant
if (!indicCycleContact) // modification de la charge et du pas de temps qu'au premier passage
// mais pas après un changement de statut
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{ bool modif_temps = Pilotage_du_temps(charge,arret_pilotage); // appel du Pilotage
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
//-- si le temps a changé il faut de nouveau appeler la gestion du pas de temps
// car il y a des grandeurs reliées au pas de temps qui y sont calculé
if (modif_temps)
this->Gestion_pas_de_temps(true,lesMail,2); // 2 signifie cas courant
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (arret_pilotage) break; // pilotage -> arret du calcul, on sort du switch
#ifdef UTILISATION_MPI
// seule le process 0 s'occupe de la sortie
if (ParaGlob::Monde()->rank() == 0)
{
#endif
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// affichage de l'increment de charge
aff_incr=pa.Vrai_commande_sortie(icharge,temps_derniere_sauvegarde); // pour simplifier
if (aff_incr)
{ cout << "\n======================================================================"
<< "\nINCREMENT DE CHARGE : " << icharge
<< " intensite " << charge->IntensiteCharge()
<< " t= " << charge->Temps_courant()
<< " dt= " << ParaGlob::Variables_de_temps().IncreTempsCourant()
<< "\n======================================================================";
};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
#ifdef UTILISATION_MPI
// seule le process 0 s'occupe de la sortie
}
else {aff_incr = false;};
#endif
// -- initialisation des coordonnees et des ddl a tdt en fonctions des
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// ddl imposes et de l'increment du chargement: change_statut sera recalculé ensuite
lesCondLim->MiseAJour_tdt
(pa.Multiplicateur(),lesMail,charge->Increment_de_Temps(),lesRef,charge->Temps_courant()
,lesCourbes1D,lesFonctionsnD,charge->MultCharge(),change_statut,cas_combi_ddl);
// dans le cas ou il y a changement de statut il faut remettre à jour
// le dimensionnement de variables intermédiaires
if (change_statut)
{ li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb
,lesRef,charge->Temps_courant(),icas);
int ttsi = li_gene_asso.size();
X_Bl.Change_taille(ttsi);V_Bl.Change_taille(ttsi);G_Bl.Change_taille(ttsi);
// ---- initialisation des sauvegardes sur matrice et second membre
// ce qui ne correspond à rien ici normalement
lesCondLim->InitSauve(Ass3.Nb_cas_assemb());
};
// --- calcul des éléments de contact: (correspond à la définition de la surface de contact)
// dans le cas où type_activation_contact ne s'effectue pas à chaque itération
if (type_activation_contact != 1)
{ // definition ou mise à jour, des elements de contact eventuels
// - imposition (en fonction de l'algorithme de contact) de conditions particulières de penetration (nulle par exemple)
if (pa.ContactType()) // traitement différent lors du premier incément calculé (premier passage) puis des passages courants
{ if (premier_calcul)
{ // ici delta_X est nulle au premier passage ... on laisse quand même par cohérence avec les autres algo
lescontacts->DefElemCont(delta_X.Max_val_abs());
} // au début il n'y a pas de déplacement à priori
else
{ lescontacts->SuppressionDefinitiveElemInactif(); // on supprime les éléments inactifs testés à l'incr prec dans Actualisation()
// lescontacts->Nouveau(lesMail->Max_var_dep_t_a_tdt()); // on rajoute éventuellement de nouveau élément de contact
lescontacts->Nouveau(delta_X.Max_val_abs()); // idem mais je pense plus rapide
};
};
};
}
else
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{ if (arret_pilotage) break; // pilotage -> arret du calcul
#ifdef UTILISATION_MPI
// seule le process 0 s'occupe de la sortie
if (ParaGlob::Monde()->rank() == 0)
#endif
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
cout << "\n============================================================================="
<< "\n ....... re-analyse du contact ........ "
<< "\nINCREMENT DE CHARGE : " << icharge << " intensite " << charge->IntensiteCharge()
<< " t= " << charge->Temps_courant()
<< " dt= " << ParaGlob::Variables_de_temps().IncreTempsCourant()
<< "\n============================================================================";
};
};
}; // fin du switch sur compteur_demarrage
mise à jour 7.010
2023-05-03 17:23:49 +02:00
#ifdef UTILISATION_MPI
temps_transfert_court.Mise_en_route_du_comptage(); // comptage cpu
broadcast(*ParaGlob::Monde(), arret_pilotage, 0);
temps_transfert_court.Arret_du_comptage(); // fin comptage cpu
#endif
if (arret_pilotage) break; // si dans le switch précédent on a un arret de pilotage qui est demandé
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
lesLoisDeComp->MiseAJour_umat_nbincr(icharge); // init pour les lois Umat éventuelles
// mise à jour éventuelle de la matrice de raideur en fonction de l'existence du contact et du type de modèle de contact
// if (pa.ContactType())
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// Mise_a_jour_Choix_matriciel_contact(mato,Ass.Nb_cas_assemb(),lescontacts);
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// --- récupération (initialisation) des ddl position, vitesse et accélération
// récupe X_t initiale
lesMail->Vect_loc_vers_glob(TEMPS_t,X1,X_t,X1);
save_X_t=X_t; // sauvegarde
// dans le cas particulier d'une succesion de sous-algo, avec non validation du précédent algo
// on utilise pour démarrer la situation finale précédemment calculée
if ((tb_combiner != NULL) && ((*tb_combiner)(3) != NULL))
{// on initialise les valeurs globales à t, par contre on ne change pas les grandeurs locales à t !!
lesMail->Vect_loc_vers_glob(TEMPS_tdt,X1,X_tdt,X1);
X_t = X_tdt;
// on récupère la situation
lesMail->Vect_loc_vers_glob(TEMPS_tdt,V1,vitesse_tdt,V1);
lesMail->Vect_loc_vers_glob(TEMPS_tdt,GAMMA1,acceleration_tdt,GAMMA1);
vitesse_t = vitesse_tdt;
acceleration_t = acceleration_tdt;
};
lesMail->Vect_loc_vers_glob(TEMPS_t,V1,vitesse_t,V1);
lesMail->Vect_loc_vers_glob(TEMPS_t,GAMMA1,acceleration_t,GAMMA1);
// on commence ici la boucle avec relaxation
int relax_vit_acce = 0; // variable de retour de la relaxation cinétique
// qui indique s'il y a relaxation ou pas, sert éventuellement pour le calcul de la masse
if (amortissement_cinetique)
{Algori::InitialiseAmortissementCinetique(); // initialisation des compteurs pour l'amortissement au cas ou
list_iter_relax.clear();
};
// on met à jour l'indicateur visqueux_activer
switch (typeCalRelaxation)
{ case 1: visqueux_activer = false; break;
case 2: visqueux_activer = true; break;
case 4: visqueux_activer = false; break; // au début en cinétique
default: break; // on ne doit pas arriver ici !! et a déjà été testé à la lecture
};
// force_recalcul_masse : un indicateur de retour pour la méthode Cinetique_ou_visqueux:
bool force_recalcul_masse= false; // par défaut on ne force pas
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// on démarre avec le compteur à 0 et on sauvegarde la position finale à l'itération 0
lesMail->Quelconque_glob_vers_local(X1,X_tdt,typQ_XI_ITER_0);
V 7.017 - intro para limite_inf_bIIb_ et limite_inf_Qeps_ dans les potentiels de grenoble - cor bug diff finie pour Qeps petit sur loi grenoble - modif prise en compte de la regularisation loi grenoble (passage systématique en différence finie pour Qeps très faible)
2023-07-06 08:42:57 +02:00
// boucle de convergence sur un increment
Vecteur * sol; // pointeur du vecteur solution
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
for (compteur = 0; (compteur<= pa.Iterations())&&(!pa.EtatSortieEquilibreGlobal()); compteur++)
//---//\\//\\// début de la boucle sur les itérations d'équilibres //\\//\\//
{
// calcul de l'increment
// initialisation des deux partie du second membre
vglobin.Zero();
vglobex.Zero();
if (pa.ContactType())
vcontact.Zero();
vglobaal.Zero(); // puissance totale
// init pour var glob
Transfert_ParaGlob_COMPTEUR_ITERATION_ALGO_GLOBAL(compteur);
// prise en compte du cas particulier ou l'utilisateur demande
// une mise à jour des conditions limites à chaque itération
if (cL_a_chaque_iteration)
{// -- initialisation des coordonnees et des ddl a tdt en fonctions des
lesCondLim->MiseAJour_tdt
(pa.Multiplicateur(),lesMail,charge->Increment_de_Temps(),lesRef,charge->Temps_courant()
,lesCourbes1D,lesFonctionsnD,charge->MultCharge(),change_statut,cas_combi_ddl);
// dans le cas ou il y a changement de statut il faut remettre à jour
// le dimensionnement de variables intermédiaires
if (change_statut)
{ li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb
,lesRef,charge->Temps_courant(),icas);
int ttsi = li_gene_asso.size();
X_Bl.Change_taille(ttsi);V_Bl.Change_taille(ttsi);G_Bl.Change_taille(ttsi);
// ---- initialisation des sauvegardes sur matrice et second membre
// ce qui ne correspond à rien ici normalement
lesCondLim->InitSauve(Ass3.Nb_cas_assemb());
};
// on récupère les nouvelles positions globalement de manière à pouvoir calculer le delta_X pour le contact
lesMail->Vect_loc_vers_glob(TEMPS_tdt,X1,X_tdt,X1);
save_X_t=X_t; // sauvegarde
};
// passage aux noeuds des vecteurs globaux: F_INT, F_EXT
Algori::Passage_aux_noeuds_F_int_t_et_F_ext_t(lesMail);
// renseigne les variables définies par l'utilisateur via les valeurs déjà calculées par Herezh
Algori::Passage_de_grandeurs_globales_vers_noeuds_pour_variables_globales(lesMail,varExpor,Ass1.Nb_cas_assemb(),*lesRef);
varExpor->RenseigneVarUtilisateur(*lesMail,*lesRef);
lesMail->CalStatistique(); // calcul éventuel de statistiques
// on met à jour la valeur de lambda
lambda = Calcul_Lambda();
lesMail->Force_Ddl_aux_noeuds_a_une_valeur(R_X1,0.0,TEMPS_tdt,true); // mise à 0 des ddl de réactions, qui sont uniquement des sorties
lesMail->Force_Ddl_etendu_aux_noeuds_a_zero(Ddl_enum_etendu::Tab_FN_FT()); // idem pour les composantes normales et tangentielles
// affichage ou non de l'itération
bool aff_iteration = (pa.freq_affich_iter() > 0) ?
mise à jour 7.010
2023-05-03 17:23:49 +02:00
(aff_incr && (compteur % pa.freq_affich_iter()==0) &&(compteur!=0)) : false ;
#ifdef UTILISATION_MPI
// seule le process 0 s'occupe de la sortie
if (ParaGlob::Monde()->rank() != 0) aff_iteration=false;
#endif
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
/* // --- imposition des ddls bloqués
// initialisation des coordonnees et des ddl a tdt en fonctions des
// ddl imposes et de l'increment du chargement et des conditions linéaires imposées
// bool change_statut = false; // init des changements de statut
// lesCondLim->MiseAJour_tdt
// (pa.Multiplicateur(),lesMail,charge->Increment_de_Temps(),lesRef,charge->Temps_courant()
// ,lesCourbes1D,charge->MultCharge(),change_statut,cas_combi_ddl);
// mise en place des conditions linéaires
// on ne met pas en place les conditions linéaires ici, car en fait on veut avoir la répercution du précédent calcul
// or, due au fait que ce sont des conditions linéaires, elles interviennent sur plusieurs ddl en même temps
// dans le repère globale. Mais quand on met les CL linéaires, on bloque une direction dans le nouveau repère !!
// or le couplage X V gamma n'est valide que dans le repère global !! donc on appplique les CL linéaires après ce couplage
// ***** pour l'instant, ce fonctionnement est ok pour les CL en gamma, il faudra regarder comment cela ce passe
// ***** dans le cas de CL linéaire en V et X */
// récupération au niveau global des ddl locaux à tdt avec conditions limite
// on ne s'occupe pas du vecteur accélération qui lui va être calculé !, car il dépend de
// l'équilibre dynamique (ici pseudo-dynamique)
// pour le vecteur accélération, seules les ddl avec CL sont différents de la précédente
// récupération
lesMail->Vect_loc_vers_glob(TEMPS_tdt,X1,X_tdt,X1);
lesMail->Vect_loc_vers_glob(TEMPS_tdt,V1,vitesse_tdt,V1);
// lesMail->Vect_loc_vers_glob(TEMPS_tdt,GAMMA1,acceleration_tdt,GAMMA1);
// maintenant on met les conditions limites sur les ddls bloqués secondaires c-a-d associés
// aux ddl bloqués par l'utilisateur, leur calcul dépend de l'algorithme d'où un calcul global
list <LesCondLim::Gene_asso>::iterator ie,iefin=li_gene_asso.end(); // def d'un iterator adoc
int ih=1; // indice
for(ie=li_gene_asso.begin(),ih=1;ie!=iefin;ie++,ih++)
// comme les valeurs des X V Gamma vont être écrasé par le calcul global, on utilise
// des conteneurs intermédiaires
{//trois cas
LesCondLim::Gene_asso & s = (*ie); // pour simplifier
int ix=s.pointe(1); // début des Xi
int iv=s.pointe(2); // début des Vi
int ig=s.pointe(3); // début des gammai
// on utilise le schéma des différences finis, mais avec un pas de temps = 1.
// centrés pour calculer les valeurs des ddl dans le cas de ddl bloqué
switch (PremierDdlFamille(s.ty_prin))
{ case X1: // cas ou les Xi sont imposés, calcul approché des V_tdt
{ X_Bl(ih) = X_tdt(ix); // on récupère la nouvelle position
// pour les vitesses en fait il faudrait la valeur à t+2*deltat, donc on fait l'approximation t+1/2 deltat
V_Bl(ih) = (X_tdt(ix)-X_t(ix));
// pour l'accélération, a priori elle n'intervient pas dans le calcul d'équilibre, on la met donc
G_Bl(ih) = acceleration_t(ig); // identique à celle à t
break;
}
case V1: // cas où les Vi sont imposées, calcul directe des Xi_tdt et Gammai_tdt
{ V_Bl(ih) = vitesse_tdt(iv);
G_Bl(ih) = (vitesse_tdt(iv)-vitesse_t(iv))*2. - acceleration_t(ig);
X_Bl(ih) = X_t(ix) + vitesse_t(iv) + 0.5 * acceleration_t(ig);
break;
}
case GAMMA1: // cas où les gammai sont imposées, calcul directe des Xi_tdt et Vi_tdt
{ G_Bl(ih) = acceleration_tdt(ig);
X_Bl(ih) = X_t(ix) + vitesse_t(iv) + 0.5 * acceleration_t(ig);
V_Bl(ih) = vitesse_t(iv)+0.5*(acceleration_tdt(ig)+acceleration_t(ig));
break;
}
default:
cout << "\n ERREUR bizarre de resolution :PremierDdlFamille(s.ty_prin)= " << Nom_ddl(PremierDdlFamille(s.ty_prin))
<< "\n AlgoriRelaxDyna::CalEquilibre(...) ";
Sortie(1);
};
};
// ---- calcul des champs de vitesse et de position
// position calculée de manière exacte par rapport au schéma
X_tdt = X_t + vitesse_t + 0.5 * acceleration_t;
// pour la vitesse, il s'agit d'une prédiction explicite, qui sera affinée après résolution
vitesse_tdt = vitesse_t + acceleration_t;
// pilotage éventuelle des positions et vitesses
Pilotage_maxi_X_V(X_t,X_tdt,vitesse_t,vitesse_tdt);
// -- maintenant on met réellement en place les CL a partir de la sauvegarde
for(ie=li_gene_asso.begin(),ih=1;ie!=iefin;ie++,ih++)
{LesCondLim::Gene_asso & s = (*ie); // pour simplifier
int ix=s.pointe(1); // début des Xi
int iv=s.pointe(2); // début des Vi
int ig=s.pointe(3); // début des gammai
X_tdt(ix) = X_Bl(ih);
vitesse_tdt(iv) = V_Bl(ih);
acceleration_tdt(ig) = G_Bl(ih);
};
// ---- passage des valeurs calculées aux niveaux des maillages
lesMail->Vect_glob_vers_local(TEMPS_tdt,X1,X_tdt,X1);
lesMail->Vect_glob_vers_local(TEMPS_tdt,V1,vitesse_tdt,V1);
// accélération à t : seules celles correspondantes au CL ont variées
// lesMail->Vect_glob_vers_local(TEMPS_t,GAMMA1,acceleration_t,GAMMA1);
lesLoisDeComp->MiseAJour_umat_nbiter(compteur); // init pour les lois Umat éventuelles
/* //debug ------
// { cout << "\n ***1 AlgoriRelaxDyna::CalEquilibre( ";
// lesMail->Noeud_LesMaille(1,483).Affiche();
// cout << endl ;
// };
////fin debug ------ */
// mise en place des conditions linéaires
lesCondLim->MiseAJour_condilineaire_tdt
(pa.Multiplicateur(),lesMail,charge->Increment_de_Temps(),lesRef,charge->Temps_courant()
,lesCourbes1D,lesFonctionsnD,charge->MultCharge(),change_statut,cas_combi_ddl);
// il n'y a pas de changement de largeur de bande pour les conditions linéaires, car on a mis la largeur maxi
// au moment de la création de la matrice masse
// recherche de nouveaux contacts si c'est demandé
if ((type_activation_contact == 1)&&(compteur_demarrage > 0))
{ if (pa.ContactType()) // réexamen du contact pour voir
// il faut mettre le contact ici, car il utilise le déplacement de t à tdt
{ lesMail->Mise_a_jour_boite_encombrement_elem_front(TEMPS_tdt); //s'il n'y a pas
if (premier_calcul)
{ lescontacts->DefElemCont(delta_X.Max_val_abs());
premier_calcul=false;} // au début il n'y a pas de déplacement à priori, on prend 2. * le delta noeud mini
else
{ lescontacts->SuppressionDefinitiveElemInactif(); // on supprime les éléments inactifs testés à l'incr prec dans Actualisation()
lescontacts->Nouveau(delta_X.Max_val_abs()); // idem mais je pense plus rapide
};
};
};
// -+-+ si la recherche de nouveaux contact n'est pas effectuée à chaque itération
// actualisation des conditions de contact qui peuvent changer la largeur de bande,
// quand un noeud glisse d'une facette sur une voisine, peut changer la position du noeud
// qui est projeté sur la facette dans le cas de l'algorithme cinématique
if (((compteur != 0)&&(pa.ContactType()))&&(compteur_demarrage > 0))
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{ lescontacts->Actualisation(0); // en particulier: pour le type 4 on a projection
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// des noeuds sur les facettes maîtres
// mise à jour éventuelle des répercussions du contact sur les noeuds en contact
AlgoriRelaxDyna::Repercussion_algo_sur_cinematique(lescontacts,X_tdt,vitesse_tdt);
// lesMail->Vect_loc_vers_glob(TEMPS_tdt,X1,X_tdt,X1);
// lesMail->Vect_loc_vers_glob(TEMPS_tdt,V1,vitesse_tdt,V1);
// mise à jour éventuelle de la matrice de raideur en fonction du contact
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// Mise_a_jour_Choix_matriciel_contact(mato,Ass.Nb_cas_assemb(),lescontacts);
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
};
// -+-+ sinon l'actualisation du contact s'effectue à la fin de l'itération (un peu plus loin)
// ---- calcul des puissances internes et externe
// appel du calcul de la puissance interne et des énergies
// dans le cas d'un calcul inexploitable arrêt de la boucle
if (!SecondMembreEnerg(lesMail,Ass1,vglobin)) break;
// calcul des maxi des puissances internes
maxPuissInt = vglobin.Max_val_abs();
F_int_tdt_prec = F_int_tdt; // sauvegarde des valeurs de l'itération précédente
F_int_tdt = vglobin; // sauvegarde des forces généralisées intérieures
if (pa.ContactType()==4) // dans le cas d'un contact de type 4
// on transfert les forces internes aux noeuds
lesMail->Quelconque_glob_vers_local(X1,F_int_tdt,typQ_gene_int);
// mise en place du chargement impose, c-a-d calcul de la puissance externe
// si pb on sort de la boucle
if (!(charge->ChargeSecondMembre_Ex_mecaSolid
(Ass1,lesMail,lesRef,vglobex,pa,lesCourbes1D,lesFonctionsnD)))
{ Change_PhaseDeConvergence(-10);break;};
//// autre essai !!
//if (pa.ContactType()==4) // dans le cas d'un contact de type 4
//{V_ext = F_int_tdt; V_ext += F_ext_tdt;
// lesMail->Quelconque_glob_vers_local(X1,V_ext,typQ_gene_int);
//};
// calcul des reactions de contact éventuelles (contact et frottement)
if ((pa.ContactType())&&(compteur_demarrage > 0)) // et des énergies développées pendant le contact
{// dans le cas où le calcul est inexploitable (pb dans le calcul) arrêt de la boucle
if (!SecondMembreEnergContact(lescontacts,Ass1,vcontact,aff_iteration)) break;
};
F_ext_tdt_prec = F_ext_tdt; // sauvegarde des valeurs de l'itération précédente
if (pa.ContactType())
vglobex += vcontact;
maxPuissExt = vglobex.Max_val_abs();
F_ext_tdt = vglobex; // sauvegarde des forces généralisées extérieures
// second membre total
vglobaal += vglobex ;vglobaal += vglobin ;
// on calcul la matrice de masse qui est très particulière dans le cas de la relaxation dynamique
// c'est la partie la plus importante de l'algo: l'adaptation de la masse en continue
//-- *** on utilise le volume de chaque élément, donc il doit avoir été calculé
// -- ** ce qui est le cas au moment du calcul des forces internes
CalculEnContinuMatriceMasse
(relax_vit_acce,lesMail,Ass1,compteur
,diversStockage,lesRef,X1,premier_calcul,lescontacts
,force_recalcul_masse,lesFonctionsnD);
// calcul des reactions de contact pour les noeuds esclaves
// dans le repere absolu ( pour la sortie des infos sur le contact)
// et test s'il y a decollement de noeud en contact (pour les contacts actifs)
// mais ici, il n'y a pas de modification des éléments de contact (elles ont été faites dans lescontacts->Actualisation())
bool decol=false; // création et init de decol
if ((pa.ContactType())&&(compteur_demarrage > 0))
lescontacts->CalculReaction(vglobin,decol,Ass1.Nb_cas_assemb(),aff_iteration);
// - definition éventuelle de conditions limites linéaires de contact, en fonction des éléments du contact existant à ce niveau
// (certains contacts (par pénalisation par exemple) ne produise pas de conditions linéaires)
const Tableau <Condilineaire>* tabCondLine=NULL;
if (((pa.ContactType()==1) || (pa.ContactType()==3))&&(compteur_demarrage > 0))
tabCondLine = &(lescontacts->ConditionLin(Ass1.Nb_cas_assemb()));
// -- dans le cas d'un amortissement visqueux critique ou
// dans le cas où l'on utilise de l'amortissement numérique le second membre est modifiée
if ((typeCalRelaxation == 2)
|| (typeCalRelaxation == 3)
|| (typeCalRelaxation == 4)
|| (pa.Amort_visco_artificielle())
)
{ if (Arret_A_Equilibre_Statique()) // si on veut un équilibre statique, on sauvegarde les forces statiques, dans (*vglob_stat) qui est différent de (*vglobaal)
// en amortissement cinétique, ce n'est pas le cas
(*vglob_stat) = vglobaal;
};
if ((typeCalRelaxation == 2)
|| (pa.Amort_visco_artificielle()&& !(Cinetique_ou_visqueux(force_recalcul_masse)))
|| (typeCalRelaxation == 3)
|| ((typeCalRelaxation == 4) && !(Cinetique_ou_visqueux(force_recalcul_masse)))
)
{
/* if (Arret_A_Equilibre_Statique()) // si on veut un équilibre statique, on sauvegarde les forces statiques
// (*vglob_stat) = (*vglobaal); */
Cal_mat_visqueux_num_expli(*mat_masse_sauve,mat_C_pt,delta_X,false,vitesse_tdt);
/* else { mat_C_pt->Initialise(0.);};
// if (typeCalRelaxation == 2)
// CalculEnContinuMatriceViscositeCritique(*mat_masse_sauve,*mat_C_pt,delta_X,vitesse_tdt); */
vglobaal -= mat_C_pt->Prod_mat_vec(vitesse_t,forces_vis_num);
};
// initialisation des sauvegardes sur second membre (uniquement pour les gammai)
lesCondLim->InitSauve(Ass3.Nb_cas_assemb());
////---- debug
//cout << "\n debug algo relax ";
//vglobaal.Affiche();
////--- fin debug
// on récupère les réactions avant changement de repère et calcul des torseurs de réaction
lesCondLim->ReacAvantCHrepere(vglobaal,lesMail,lesRef,Ass3.Nb_cas_assemb(),cas_combi_ddl);
// -->> expression de la matrice masse et du second membre dans un nouveau repere
// mais ici on n'impose pas les conditons, on fait simplement le changement de repère
//modif_repere = indicateur si réellement il y a un changement
bool modif_repere = lesCondLim->CoLinCHrepere_int(*mat_masse,vglobaal,Ass3.Nb_cas_assemb(),vglob_stat);
if ((pa.ContactType()==1)&&(compteur_demarrage > 0)) // idem pour le contact conduisant à des conditions linéaires
modif_repere = modif_repere ||
lesCondLim->CoLinCHrepere_ext(*mat_masse,vglobaal,*tabCondLine,Ass3.Nb_cas_assemb(),vglob_stat);
// sauvegarde des reactions pour les ddl bloques (simplement)
// ***dans le cas statique il semble (cf. commentaire dans l'algo) que ce soit inutile donc a voir
// ***donc pour l'instant du a un bug je commente
// sauvegarde des reactions aux ddl bloque (uniquement pour les Xi)
// lesCondLim->ReacApresCHrepere(vglobin,lesMail,lesRef,Ass3.Nb_cas_assemb(),cas_combi_ddl);
// mise en place des conditions limites sur la matrice masse et le second membre
lesCondLim->ImposeConLimtdt(lesMail,lesRef,*mat_masse,vglobaal
,Ass3.Nb_cas_assemb(),cas_combi_ddl,vglob_stat);
// puis on prépare (si besoin est en fonction du type de matrice) la résolution
// s'il n'y a pas de condition linéaire et pas de chg de repère,
// on peut déjà inverser par exemple la matrice masse
if ((!modif_repere) && (!(lesCondLim->ExisteCondiLimite())))
mat_masse->Preparation_resol();
// (blocage de toutes les conditions lineaires, quelque soit leur origines ext ou int donc contact éventuel)
lesCondLim->CoLinBlocage(*mat_masse,vglobaal,Ass3.Nb_cas_assemb(),vglob_stat);
// calcul du maxi des reactions (pour les xi)
maxReaction = lesCondLim->MaxEffort(inReaction,Ass3.Nb_cas_assemb());
// sortie d'info sur l'increment concernant les réactions
if ( aff_iteration)
{ if (ParaGlob::NiveauImpression() > 2) cout << "\n IT:" << compteur;
InfoIncrementReac(lesMail,inReaction,maxReaction,Ass3.Nb_cas_assemb());
};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
bool arretResidu = false; // pour gérer le cas particulier ou on veut un arrêt et sur le résidu et sur le déplacement
bool demande_de_break=false; // pour gestion du break en tenant compte ou non du MPI
#ifdef UTILISATION_MPI
// seule le process 0 s'occupe de la convergence
if (ParaGlob::Monde()->rank() == 0)
{
#endif
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// examen de la convergence si nécessaire, utilisant le résidu
if (ArretEquilibreStatique() && (compteur>1) && (compteur_demarrage != 0) )// cas d'une convergence en utilisant le résidu
{ double toto=0.; int itera = 0; // valeur par defaut pour ne pas se mettre dans un cas itératif de type algo de Newton
bool arret_demande = false; // normalement n'intervient pas ici, car il n'y a pas de prise en compte d'iteration
// bool affiche = aff_incr && (ParaGlob::NiveauImpression() > 2);
arret = Convergence(aff_iteration ,toto,vglobaal,maxPuissExt,maxPuissInt,maxReaction,itera,arret_demande);
if (arret)
{ // sortie des itérations sauf si l'on est en loi simplifiée
if (lesLoisDeComp->Test_loi_simplife() )
{lesLoisDeComp->Loi_simplifie(false); // cas d'une loi simplifié on remet normal
arret = false;
}
else
{if ((ArretEquilibreStatique() == 2)
&& ((typeCalRelaxation == 1)||(typeCalRelaxation == 3)
|| ((typeCalRelaxation == 4) && (Cinetique_ou_visqueux(force_recalcul_masse)))
)
)
{ arretResidu=true;} // cas relaxation avec amortissement cinétique
mise à jour 7.010
2023-05-03 17:23:49 +02:00
else {demande_de_break=true;}//break;}; // cas normal,
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
};
};
};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
#ifdef UTILISATION_MPI
};
temps_transfert_court.Mise_en_route_du_comptage(); // comptage cpu
broadcast(*ParaGlob::Monde(), arretResidu, 0);
broadcast(*ParaGlob::Monde(), arret, 0);
temps_transfert_court.Arret_du_comptage(); // fin comptage cpu
#endif
if (demande_de_break)
break;
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// 4 --<ARD>-- calcul des nouvelles accélérations
mise à jour 7.010
2023-05-03 17:23:49 +02:00
#ifdef UTILISATION_MPI
// seule le process 0 fait la résolution globale
if (ParaGlob::Monde()->rank() == 0)
{
#endif
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// resolution simple (fonction du type de matrice)
// ou non suivant modif_repere
tempsResolSystemLineaire.Mise_en_route_du_comptage(); // temps cpu
residu_final = vglobaal; // sauvegarde pour le post-traitement
// ici il y a deux conditions équivalentes, mais on les laisse car on pourrait peut-être
// avoir un changement de repère sans
if ((!modif_repere) && (!(lesCondLim->ExisteCondiLimite())))
{mat_masse->Simple_Resol_systID_2 (vglobaal,acceleration_tdt,pa.Tolerance()
,pa.Nb_iter_nondirecte(),pa.Nb_vect_restart());}
else // cas où on a eu des changements de repère, il faut une résolution complète
{acceleration_tdt = (mat_masse->Resol_systID(vglobaal,pa.Tolerance()
,pa.Nb_iter_nondirecte(),pa.Nb_vect_restart()));};
tempsResolSystemLineaire.Arret_du_comptage(); // temps cpu
// affichage éventuelle du vecteur solution : accélération
if (ParaGlob::NiveauImpression() >= 10)
{ string entete = " affichage du vecteur solution acceleration ";
acceleration_tdt.Affichage_ecran(entete); };
// retour des accélération dans les reperes generaux, dans le cas où ils ont ete modifie
// par des conditions linéaires
lesCondLim->RepInitiaux( acceleration_tdt,Ass3.Nb_cas_assemb());
mise à jour 7.010
2023-05-03 17:23:49 +02:00
#ifdef UTILISATION_MPI
}
else // s'il s'agit d'un process de calcul élémentaire
{sol = &vglobaal; // il faut affecter sol pour récupérer ensuite la solution
};
// le process 0 transmet aux autres process le vecteur résultat
temps_transfert_long.Mise_en_route_du_comptage(); // comptage cpu
sol->Broadcast(0);
temps_transfert_long.Arret_du_comptage(); // fin comptage cpu
#endif
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// effacement du marquage de ddl bloque du au conditions lineaire imposée par l'entrée
lesCondLim->EffMarque();
if (pa.ContactType()) lescontacts->EffMarque();
// -- calcul de la vitesse exacte à t+dt selon le schéma DFC avec un delta t =1.
vitesse_tdt = vitesse_t + 0.5 * (acceleration_t+acceleration_tdt);
// -- on remet réellement en place les CL a partir de la sauvegarde sur les vitesses (seulement)
for(ie=li_gene_asso.begin(),ih=1;ie!=iefin;ie++,ih++)
{LesCondLim::Gene_asso & s = (*ie); // pour simplifier
vitesse_tdt(s.pointe(2)) = V_Bl(ih);
};
// passage des accélérations et des vitesses calculées aux niveaux des maillages
lesMail->Vect_glob_vers_local(TEMPS_tdt,GAMMA1,acceleration_tdt,GAMMA1);
lesMail->Vect_glob_vers_local(TEMPS_tdt,V1,vitesse_tdt,V1);
// dans le cas ou la recherche de nouveaux contacts est effectuée à chaque itération
if (((type_activation_contact== 1) && (pa.ContactType()))&&(compteur_demarrage > 0))
{ // actualisation des éléments de contact et éventuellement inactivation d'éléments
mise à jour 7.010
2023-05-03 17:23:49 +02:00
lescontacts->Actualisation(0); // si on n'a plus de projection
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// on inactive les éléments de contact qui se relache: testé soit via la réaction
lescontacts->RelachementNoeudcolle(); // ou via la sortie d'une zone d'accostage (dépend de l'algo)
};
if (compteur_demarrage != -1)
{// mise à jour de delta_X, var_delta_X et passage en global des maxi
// delta_X = X_tdt; delta_X -= X_t; // X_tdt - X_t
Algori::Cal_Transfert_delta_et_var_X(max_delta_X,max_var_delta_X);
CalEnergieAffichage(1.,vitesse_tdt,*mat_masse_sauve,delta_X
,icharge,brestart,acceleration_tdt,forces_vis_num);
if (compteur_demarrage==0)// dans le cas du premier incrément on considère que la balance vaut l'énergie
// cinétique initiale, car vu que l'on ne met pas de CL à t=0, E_cin_0 est difficile à calculer
{E_cin_0 = E_cin_tdt - bilan_E + E_int_tdt - E_ext_tdt; };
};
// calcul de l'amortissement cinétique dans le cas ou ce n'est pas de l'amortissement visqueux critique
if ((typeCalRelaxation == 1)||(typeCalRelaxation == 3)
|| ((typeCalRelaxation == 4) && (Cinetique_ou_visqueux(force_recalcul_masse)))
)
{ relax_vit_acce = AmortissementCinetique(delta_X,1.,X_tdt,*mat_masse_sauve,compteur,vitesse_tdt);
// il faut re-updater les vitesses
// if (Abs(relax_vit_acce) == 1)
if (relax_vit_acce == 1) // il y a eu relaxation
{lesMail->Vect_glob_vers_local(TEMPS_tdt,V1,vitesse_tdt,V1);
list_iter_relax.push_front(compteur); // on sauvegarde l'iter relaxé
};
// examen de la convergence éventuelle, utilisant le déplacement et/ou le résidu
if (Pilotage_fin_relaxation_et_ou_residu(relax_vit_acce,1,compteur,arretResidu,arret))
break;
};
// on incrémente les valeurs globales à t, par contre on ne change pas les grandeurs locales à t !!
X_t = X_tdt;
vitesse_t = vitesse_tdt;
acceleration_t = acceleration_tdt;
// dans le cas du mode debug on sort éventuellement les infos au fil du calcul (un peu bricolé)
if ((mode_debug > 0)||(pa.EtatSortieEtatActuelDansCVisu()))
{ bool a_sortir_debug=false; // pour éviter une division par 0 du test (compteur % mode_debug == 0)
if (pa.EtatSortieEtatActuelDansCVisu()) {a_sortir_debug=true;}
else if (compteur % mode_debug == 0) {a_sortir_debug=true;};
if (a_sortir_debug)
{ // on passe les déplacements de tdt à t
lesMail->Vect_glob_vers_local(TEMPS_t,X1,X_tdt,X1);
// visualisation éventuelle au fil du calcul: essentiellement la déformée
tempsCalEquilibre.Arret_du_comptage(); // on arrête le compteur pour la sortie
VisuAuFilDuCalcul(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage,charge
,lesCondLim,lescontacts,resultats,type_incre,(icharge*1000000+compteur));
// on remet les choses dans l'ordre initial
tempsCalEquilibre.Mise_en_route_du_comptage(); // on remet en route le compteur
lesMail->Vect_glob_vers_local(TEMPS_t,X1,save_X_t,X1);
};
};
// cas du calcul des énergies, passage des grandeurs de tdt à t
if (compteur_demarrage != -1) Algori::TdtversT();
// si compteur_demarrage = -1 ou 0 on a un traitement spécial, on sort de la boucle
// car ce sont en fait des passages d'initialisation
if (compteur_demarrage < 1)
break;
//---//\\//\\// fin de la boucle sur les itérations d'équilibres //\\//\\//
};
// effacement du marquage de ddl bloque du au conditions lineaire imposée par l'entrée
// car ici on n'est pas passé par cette endroit si l'on a convergé
lesCondLim->EffMarque();
if (pa.ContactType()) lescontacts->EffMarque();
// mise à jour des indicateurs contrôlés par le tableau: *tb_combiner
if (tb_combiner != NULL) // cas d'un contrôle via des fonctions nD
{if ((*tb_combiner)(1) != NULL)
validation_calcul = (*tb_combiner)(1)->Valeur_pour_variables_globales()(1);
if ((*tb_combiner)(2) != NULL)
sortie_boucle_controler_par_fctnD = (*tb_combiner)(2)->Valeur_pour_variables_globales()(1);
};
// gestion de la fin des itérations
mise à jour 7.010
2023-05-03 17:23:49 +02:00
#ifdef UTILISATION_MPI
// seule le process 0 a fait la résolution globale
// il gère seul également la convergence, mais il doit tenir au courant les autres process
Algori::Passage_indicConvergenceAuxProcCalcul();
// ce qui permet le déroulement correct de la suite pour tous les process
#endif
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
if ( ((compteur_demarrage > 0) && (!Pilotage_fin_iteration_implicite(compteur))))
{ // cas d'une non convergence
pas_de_convergence_pour_l_instant = 1;
// comme on incrémente pas les ddl on remet cependant les ddl
// et les grandeurs actives de tdt aux valeurs de ceux de t
lesMail->TversTdt();
// ------ cas particulier où on a une divergence qui demande de remonter sur plus d'un incrément
if (Algori::PhaseDeConvergence() == -8)
{ int nb_incr_en_arriere = 3; // !!!!! nombre actuellement arbitraire -> par la suite mettre dans les para
if (Controle_retour_sur_un_increment_enregistre(nb_incr_en_arriere,icharge))
{ // cas ou on a réussi à trouver un incrément sauvegardé adoc = maintenant icharge
int cas=2;
this->Lecture_base_info(cas ,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage
,charge,lesCondLim,lescontacts,resultats,icharge);
// comme les conditions limites cinématiques peuvent être différentes en restart
// on libére toutes les CL imposées éventuelles
lesMail->Libere_Ddl_representatifs_des_physiques(LIBRE);
lesMail->ChangeStatut(cas_combi_ddl,LIBRE);
// dans le cas d'un calcul axisymétrique on bloque le ddl 3
if (ParaGlob::AxiSymetrie())
mise à jour 7.010
2023-05-03 17:23:49 +02:00
lesMail->Inactive_un_ddl_particulier(X3);
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// on valide l'activité des conditions limites et condition linéaires, pour le temps initial
// en conformité avec les conditions lues (qui peuvent éventuellement changé / aux calcul qui a donné le .BI)
lesCondLim->Validation_blocage (lesRef,charge->Temps_courant());
li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb,lesRef,charge->Temps_courant(),icas);
int ttsi = li_gene_asso.size();
X_Bl.Change_taille(ttsi);V_Bl.Change_taille(ttsi);G_Bl.Change_taille(ttsi);
// on remet à jour les éléments pour le contact s'il y du contact présumé
if (pa.ContactType())
lesMail->Mise_a_jour_boite_encombrement_elem_front(TEMPS_t);
brestart = true; // on signale que l'on repars avec un restart
};
// sinon on ne fait rien, on se contente du pilotage de divergence existant
};
// ------ fin cas particulier où on a une divergence qui demande de remonter sur plus d'un incrément
}
else
{// --- sinon calcul correcte
pas_de_convergence_pour_l_instant = 0;
if ((ParaGlob::NiveauImpression() > 0)&&(compteur_demarrage > 0))
cout << " \n ... convergence en " << compteur << " iterations "<< endl ;
// traitement du contact dans le cas où son activation n'a pas été faite à chaque itération
if ((pa.ContactType()) && (type_activation_contact != 1))
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{ lescontacts->Actualisation(0); // actualisation du contact en fonction du dernier incrément
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// mise à jour éventuelle de la matrice de raideur en fonction du contact
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// Mise_a_jour_Choix_matriciel_contact(mato,Ass.Nb_cas_assemb(),lescontacts);
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// réexamen du contact pour voir s'il n'y a pas de nouveau element de contact
// en fait on fera au plus deux passages supplémentaire, sinon la boucle peut être infini,
// à la fin du second passage, on regarde s'il y a décollement, si oui on relâche et on refait un passage
// sinon on valide
//I)--------------
//cout << "\n avant test: indicCycleContact= " << indicCycleContact << endl;
if (indicCycleContact == 0 )
{ if (lescontacts->Nouveau(lesMail->Max_var_dep_t_a_tdt()))
{indicCycleContact=1;} // on a de nouveau contact on refait le deuxième cycle
else
// sinon, on n'a pas de nouveau contact, on regarde s'il y a du relachement
{ indicCycleContact=2;};
}
else if (indicCycleContact == 1)
{indicCycleContact=2;} // pour regarder le relachement
else
{indicCycleContact=0;}; // pas de newcontact, ni de relachement donc c'est ok
//II)---------------
if (indicCycleContact == 2)
{ if (lescontacts->RelachementNoeudcolle())
{indicCycleContact=3;} // pour une dernière boucle d'équilibre
else
{indicCycleContact=0;};// pas de relachement donc ok
};
//cout << "\n après test: indicCycleContact= " << indicCycleContact << endl;
}
else
// concerne ici: soit le cas où c'est un problème sans contact (=contact non activé)
// soit le cas où l'activation du contact s'effectue à chaque itération
{ indicCycleContact = 0; }; //pour le prochain increment
if (!(pa.ContactType()) || (indicCycleContact == 0))
{// on impose une accélération nulle compte tenu du fait qu'elle est ici purement numérique
acceleration_tdt.Zero();
vitesse_tdt.Zero(); // idem vitesses (pourrait poser des pb si l'on veut utiliser des vrais vitesses aux noeuds ?)
// passage des accélérations et des vitesses calculées aux niveaux des maillages
lesMail->Vect_glob_vers_local(TEMPS_tdt,GAMMA1,acceleration_tdt,GAMMA1);
lesMail->Vect_glob_vers_local(TEMPS_tdt,V1,vitesse_tdt,V1);
// si on est sans validation, on ne fait rien, sinon on valide l'incrément
// avec sauvegarde éventuelle
// lesMail->TdtversT(); // ******* essai
// // on valide l'activité des conditions limites et condition linéaires
// lesCondLim->Validation_blocage (lesRef,charge->Temps_courant());
if (validation_calcul)
{// actualisation des ddl et des grandeurs actives de t+dt vers t
lesMail->TdtversT();
lescontacts->TdtversT();
// cas du calcul des énergies, passage des grandeurs de tdt à t
if (compteur_demarrage != -1) Algori::TdtversT();
// actualisation des éléments de contact et éventuellement suppression
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (pa.ContactType()) lescontacts->Actualisation(0); // des éléments qui ne sont plus en contact
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
if (compteur_demarrage != -1)
{// on valide l'activité des conditions limites et condition linéaires
lesCondLim->Validation_blocage (lesRef,charge->Temps_courant());
//s'il y a remonté des sigma et/ou def aux noeuds et/ou calcul d'erreur
bool change =false; // calcul que s'il y a eu initialisation
if(prepa_avec_remont) {change = Algori::CalculRemont(lesMail,type_incre,icharge);};
if (change) // dans le cas d'une remonté il faut réactiver les bon ddls
{lesMail->Inactive_ddl(); lesMail->Active_un_type_ddl_particulier(tenuXVG);};
// 3 mai 2019: je commente la suite, car 1) ne sert pas ici car on initialise jamais le prochain
// incr avec le précédent, 2) cela change le delta_X qui ne représente plus le vrai delta X du coup
// tous les calculs qui utilise ce delta X sont faux (matrice viscosité par exemple !)
// // on regarde si l'on doit initialiser le prochain incrément avec le pas précédent
// if (pa.IniIncreAvecDeltaDdlPrec()) // si oui on sauvegarde le delta_X actuel / delta_t
// {// récupération (initialisation) des ddl position à t et t+dt
// lesMail->Vect_loc_vers_glob(TEMPS_t,X1,X_t,X1);
// lesMail->Vect_loc_vers_glob(TEMPS_tdt,X1,X_tdt,X1);
// // calcul de la variation de ddl / delta t
// delta_X = X_tdt; delta_X -= X_t; delta_X /= pa.Deltat();
// };
// cas du calcul des énergies, passage des grandeurs de tdt à t
Algori::TdtversT();
if (compteur_demarrage > 0)
{
tempsCalEquilibre.Arret_du_comptage(); // on arrête le compteur pour la sortie
// sauvegarde de l'incrément si nécessaire
Ecriture_base_info(2,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage,charge
,lesCondLim,lescontacts,resultats,type_incre,icharge);
// enregistrement du num d'incrément et du temps correspondant
list_incre_temps_calculer.push_front(Entier_et_Double(icharge,pa.Variables_de_temps().TempsCourant()));
// visualisation éventuelle au fil du calcul
VisuAuFilDuCalcul(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage,charge
,lesCondLim,lescontacts,resultats,type_incre,icharge);
tempsCalEquilibre.Mise_en_route_du_comptage(); // on remet en route le compteur
};
}; //- fin de la condition if (compteur_demarrage != -1)
}; // -- fin du test: if (validation_calcul)
if (compteur_demarrage != -1) // de nouveau le test qui est indépendant de la valeur de "validation_calcul" (??)
premier_calcul = false;// -- on enregistre le fait que l'on n'est plus au premier passage
brestart = false; // dans le cas où l'on était en restart, on passe l'indicateur en cas courant
if (validation_calcul)
{ if (compteur_demarrage > 0) // si > 0 : calcul normal donc on incrémente icharge
icharge++;
// init de var glob
Transfert_ParaGlob_COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL(icharge);
};
}; // -- fin du cas ou le calcul converge
}; // fin du calcul correcte : if // gestion de la fin des itérations
// cas particulier où la sortie de la boucle est pilotée
if (sortie_boucle_controler_par_fctnD
&& (!pas_de_convergence_pour_l_instant)
&& (compteur_demarrage > 0)
)
break;
compteur_demarrage++; // dans tous les cas on incrémente le compteur de démarrage
////------- debug
//cout << "\n debug calequilibre relaxdyna";
//cout << "\n (!charge->Fin(icharge,!pas_de_convergence_pour_l_instant))= "<< (!charge->Fin(icharge,!pas_de_convergence_pour_l_instant))
// << "\n pas_de_convergence_pour_l_instant= "<< pas_de_convergence_pour_l_instant
// << " compteur_demarrage= "<< compteur_demarrage
// << " charge->Fin(icharge,true= "<< charge->Fin(icharge,true)
// << " charge->Fin(icharge,false)= "<< charge->Fin(icharge,false)
// << " (!pa.EtatSortieEquilibreGlobal())= "<< (!pa.EtatSortieEquilibreGlobal())
// << "\n le resultat = " << ( ((!charge->Fin(icharge,!pas_de_convergence_pour_l_instant))
// || pas_de_convergence_pour_l_instant ||(compteur_demarrage < 1)
// )
// && (charge->Fin(icharge,true)!=2) // si on a dépassé le nombre d'incrément permis on s'arrête dans tous les cas
// && (charge->Fin(icharge,false)!=3) // idem si on a dépassé le nombre d'essai d'incrément permis
// // 1er appel avec true: pour affichage et second avec false car c'est déjà affiché
// && (!pa.EtatSortieEquilibreGlobal())
// )
// << "\n\n\n "
// << endl;
//--- fin debug
// while ( ((!charge->Fin(icharge,!pas_de_convergence_pour_l_instant)) // on n'affiche la sin
// // que si on a eu convergence
// || pas_de_convergence_pour_l_instant ||(compteur_demarrage < 1)
// )
// && (charge->Fin(icharge,true)!=2) // si on a dépassé le nombre d'incrément permis on s'arrête dans tous les cas
// && (charge->Fin(icharge,false)!=3) // idem si on a dépassé le nombre d'essai d'incrément permis
// // 1er appel avec true: pour affichage et second avec false car c'est déjà affiché
// && (!pa.EtatSortieEquilibreGlobal())
// )
}; // -- fin du while sur les incréments de charge
// test de fin de calcul effectue dans charge via : charge->Fin()
// on remet à jour le nombre d'incréments qui ont été effectués qui a été incrémenté en fin de while
// donc qui ne donne pas le vrai incrément futur
if (validation_calcul)
{if (compteur_demarrage > 0) // si <= 1 cela signifie qu'il n'a pas eu d'incrémentation auparavant
icharge--;
Transfert_ParaGlob_COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL(icharge);
// dans le cas d'une suite de sous algo on signale qu'il y a validation
if ((tb_combiner != NULL) && ((*tb_combiner)(1) != NULL))
(*tb_combiner)(3) = NULL;
}
else // si on ne valide pas le calcul, on reinitialise la charge
// c-a-d l'avancement en temps, incrément et facteur multiplicatif
// de manière à avoir les mêmes conditions de départ pour le prochain calcul
{ charge->Precedant(true);
// dans le cas d'une suite de sous algo on signale qu'il n'y a pas validation
if ((tb_combiner != NULL) && ((*tb_combiner)(1) != NULL))
(*tb_combiner)(3) = (*tb_combiner)(1);
} ;
////------ debug
//cout << "\n ********* debug calequilibre RD: "
// << " double deltat_actuel = pa.Deltat()= " << pa.Deltat() << flush;
////--- fin debug
tempsCalEquilibre.Arret_du_comptage(); // temps cpu
};
// dernière passe
void AlgoriRelaxDyna::FinCalcul(ParaGlob * paraGlob,LesMaillages * lesMail,
LesReferences* lesRef,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD
,VariablesExporter* varExpor
,LesLoisDeComp* lesLoisDeComp,DiversStockage* diversStockage,
Charge* charge,LesCondLim* lesCondLim,LesContacts* lescontacts
,Resultats* resultats)
{ // passage finale dans le cas d'une visualisation au fil du calcul
Transfert_ParaGlob_ALGO_GLOBAL_ACTUEL(RELAX_DYNA); // transfert info
type_incre = OrdreVisu::DERNIER_INCRE;
VisuAuFilDuCalcul(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage,charge
,lesCondLim,lescontacts,resultats,type_incre,icharge);
// sauvegarde de l'incrément si nécessaire
Ecriture_base_info(2,lesMail,lesRef,lesCourbes1D,lesFonctionsnD
,lesLoisDeComp,diversStockage,charge,lesCondLim,lescontacts
,resultats,type_incre,icharge);
// enregistrement du num d'incrément et du temps correspondant
list_incre_temps_calculer.push_front(Entier_et_Double(icharge,pa.Variables_de_temps().TempsCourant()));
};
// initialisation pour le calcul de la matrice masse dans le cas de l'algorithme de relaxation dynamique
void AlgoriRelaxDyna::InitCalculMatriceMasse(LesMaillages * lesMail
,Mat_abstraite& mat_mass,Assemblage& Ass
,Mat_abstraite& mat_mass_sauve
,LesFonctions_nD* lesFonctionsnD)
{
switch (type_calcul_mass)
{ case 1: // cas historique
{// on définit un vecteur de dimension le nombre de noeud = celui de mat_mass / dim
v_mass.Change_taille(mat_mass.Nb_ligne() / ParaGlob::Dimension());
// boucle sur les elements et on ne fait qu'une initialisation
int nbMailMax = lesMail->NbMaillage();
for (int nbMail =1; nbMail<= nbMailMax; nbMail++)
{ int nemax = lesMail->Nombre_element(nbMail);
for (int ne=1; ne<= nemax;ne++)
{ Element & el = lesMail->Element_LesMaille(nbMail,ne); // l'element
((ElemMeca&) el).InitCalculMatriceMassePourRelaxationDynamique (casMass_relax);
};
};
break;
}
case 2: // ---cas ou on utilise la matrice de raideur réelle
// premier dimentionnement
{mat_mass.Initialise(0.);
// ici le vecteur v_mass à la taille du nb de ddl
v_mass.Change_taille(mat_mass.Nb_ligne());
v_mass.Zero(); // idem: un vecteur intermédiaire de calcul
if ( ((casMass_relax > 3) && (casMass_relax < 6)) //4-5 cas où on a besoin de deux assemblages
||((casMass_relax > 7) && (casMass_relax < 10))
)
v_mass1.Change_taille(mat_mass.Nb_ligne() / ParaGlob::Dimension());
break;
};
default:
{ cout << "\n **** erreur le cas type_calcul_mass= " << type_calcul_mass
<< "\n pour l'instant n'est pas implante "
<< "\n AlgoriRelaxDyna::InitCalculMatriceMasse( ... ";
Sortie(1);
};
};
};
// calcul de la matrice masse dans le cas de l'algorithme de relaxation dynamique
void AlgoriRelaxDyna::CalculMatriceMasse
(const int& relax_vit_acce,LesMaillages * lesMail,Assemblage& Ass
, int compteur, const DiversStockage* divStoc,LesReferences* lesRef
,const Enum_ddl & N_ddl ,LesContacts* lescontacts
,LesFonctions_nD* lesFonctionsnD)
{
int dima = ParaGlob::Dimension();
// dans le cas où on est en axisymétrie il faut diminuer de 1
if (ParaGlob::AxiSymetrie())
dima--;
int nbddl_X = mat_masse->Nb_ligne();
int nbNoe = nbddl_X / dima;
switch (type_calcul_mass)
{ case 1: // cas historique, on effectue le calcul effectif de la matrice masse
{// on initialise la matrice masse à 0
mat_masse->Initialise(0.);
// vue que le v_mass a une taille qui n'est pas la même pour tous les paramétres, on vérifie
// on définit un vecteur de dimension le nombre de noeud = celui de mat_mass / dim
v_mass.Change_taille(nbNoe);
v_mass.Zero(); // idem: un vecteur intermédiaire qui sert pour la construction de la matrice masse
// de dimension le nombre de noeud
// initialisation des updates du ddl_étendu de masse sur les noeuds: on met à 0 le ddl_etendu correspondant,
// le compteur, comptant le nombre de fois où le noeuds est modifiés, est mis à 0
lesMail->InitUpdateAuNoeud(masseRelax_1);
// --- boucle sur les elements, pour définir la grandeur "masse_relax_dyn" aux noeuds de chaque éléments
for (int nbMail =1; nbMail<= lesMail->NbMaillage(); nbMail++)
for (int ne=1; ne <= lesMail->Nombre_element(nbMail); ne++)
{ Element & el = lesMail->Element_LesMaille(nbMail,ne); // l'element
// Tableau<Noeud *>& taN = el.Tab_noeud(); // tableau de noeuds de l'el
// on calcul aux noeuds, la valeur des masses nécessaires
((ElemMeca&) el).CalculMatriceMassePourRelaxationDynamique (alph,beta,lambda,gamma,theta,casMass_relax);
// cout << "numéro =" << ne << endl;
// mat_masse_loc->Affiche();
// mat_masse_loc->Affiche1(1,8,1,1,8,1);
};
// si on est dans le cas 3 ou 4, on fait la moyenne aux noeuds
if ((casMass_relax == 3) || (casMass_relax == 4))
lesMail->MoyenneCompteurAuNoeud(masseRelax_1);
// maintenant on balaie les noeuds pour définir les valeurs de la matrice masse
// on boucle sur les noeuds du maillage, on regarde si les ddl de déplacement sont actifs
// si oui, on renseigne la matrice masse
lesMail->Scalaire_loc_vers_glob(GAMMA1,v_mass,masseRelax_1);
// construction des termes de la matrice
for (int inn = 1; inn<= nbNoe; inn++)
switch (dima)
// { case 3: mat_d.set_element((inn-1)*dima+3,(inn-1)*dima+3,v_mass(inn));
// case 2: mat_d.set_element((inn-1)*dima+2,(inn-1)*dima+2,v_mass(inn));
// case 1: mat_d.set_element((inn-1)*dima+1,(inn-1)*dima+1,v_mass(inn));
// };
{ case 3: (*mat_masse)((inn-1)*dima+3,(inn-1)*dima+3) = v_mass(inn);
case 2: (*mat_masse)((inn-1)*dima+2,(inn-1)*dima+2) = v_mass(inn);
case 1: (*mat_masse)((inn-1)*dima+1,(inn-1)*dima+1) = v_mass(inn);
};
// transfert dans la matrice de masse
// MatDiag & mat_d = *((MatDiag*) & mat_masse); // récup de la matrice diagonale
break;
}
case 2: // ---cas ou on utilise la matrice de raideur réelle
// **** le calcul effectif de la matrice masse se fait ici
{mat_masse->Initialise(0.);
v_mass.Change_taille(nbddl_X);
v_mass.Zero(); // idem: v_mass1 un vecteur intermédiaire de calcul
if ( ((casMass_relax > 3) && (casMass_relax < 6)) //4-5 cas où on a besoin de deux assemblages
||((casMass_relax > 7) && (casMass_relax < 10))
)
{v_mass1.Change_taille(nbddl_X);v_mass1.Zero();};
if ((casMass_relax > 5) && (casMass_relax < 10))
matglob->Initialise(0.);
// boucle sur les elements
int nbMailMax = lesMail->NbMaillage();
for (int nbMail =1; nbMail<= nbMailMax; nbMail++)
{ int nemax = lesMail->Nombre_element(nbMail);
for (int ne=1; ne<= nemax;ne++)
{ //calcul de la raideur local et du residu local
ElemMeca & el = *((ElemMeca*) &lesMail->Element_LesMaille(nbMail,ne)); // l'element
Tableau<Noeud *>& taN = el.Tab_noeud(); // tableau de noeuds de l'el
////------- debug ---
// if (ne == 15)
// cout << "\n debug AlgoriRelaxDyna::CalculMatriceMasse( ";
////resu.raid->Affiche();cout << endl;
////------ fin debug -----
Element::ResRaid resu = el.Calcul_implicit(pa);
////------- debug ---
//cout << "\n debug AlgoriRelaxDyna::CalculMatriceMasse( ";
//resu.raid->Affiche();cout << endl;
////------ fin debug -----
// assemblage
if (casMass_relax < 2) // 0 1 : on ne retient que la diagonale
{ Ass.AssembDiagonale (v_mass,*(resu.raid),el.TableauDdl(),taN);}
else if (casMass_relax < 6) // sinon c'est la somme des valeurs absolues de la ligne
{ Ass.AssembDiagoMajorValPropre (v_mass,*(resu.raid),el.TableauDdl(),taN);}
// maintenant on s'occupe du cas où on utilise une matrice de raideur complète assemblée
else if ((casMass_relax > 5) && (casMass_relax < 10))
{if (pa.Symetrie_matrice())
Ass.AssembMatSym (*matglob,*(resu.raid),el.TableauDdl(),taN); // de la raideur
else
Ass.AssembMatnonSym (*matglob,*(resu.raid),el.TableauDdl(),taN); // de la raideur
}
else
{ cout << "\n *** attention le cas casMass_relax = " << casMass_relax
<< " n'est pas encore pris en compte !! se plaindre !! ";
Sortie(1);
};
// cas où on fait un deuxième assemblages
if ((casMass_relax > 3) && (casMass_relax < 6)) //4-5
{ Ass.AssembDiagonale (v_mass1,*(resu.raid),el.TableauDdl(),taN);};
////------- debug ---
//{cout << "\n debug AlgoriRelaxDyna::CalculMatriceMasse( ";
// cout << "\n -- info vecteur masse diagonale sur la raideur uniquement -- ";
// Coordonnee3 moy=v_mass.MinMaxMoy(true);
//};
////------ fin debug -----
};
};
if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 5 )) || (permet_affichage > 2))
{// on affiche les min max du vecteur
cout << "\n -- info vecteur masse diagonale sur la raideur uniquement -- ";
Coordonnee3 moy=v_mass.MinMaxMoy(true);
};
// prise en compte éventuelle de la raideur des éléments de contact
if (pa.ContactType()) // et des énergies développées pendant le contact
{// on récupère la liste des éléments de contact
LaLIST <ElContact>& listElContact = lescontacts->LesElementsDeContact();
LaLIST <ElContact>::iterator il,ilfin=listElContact.end();
// on gère les exceptions éventuelles en mettant le bloc sous surveillance
try
{// boucle sur les elements de contact
for (il=listElContact.begin();il != ilfin; il++)
{
if ((*il).Actif()) // on n'intervient que si le contact est actif
{//calcul de la raideur locale et du residu local
ElContact& elcontact = (*il); // pour simplifier
Tableau<Noeud *>& taN = elcontact.TabNoeud_pour_assemblage(); // tableau de noeuds
// calcul effectif de la raideur et résidu locals
Element::ResRaid resu = elcontact.SM_K_charge_contact();
// cout << "numéro =" << ne << endl;
// (resu.res)->Affiche();
// (resu.raid)->Affiche(); //1,8,1,1,8,1);
// assemblage
// dans le cas où le retour est un pointeur nul, cela signifie qu'il n'y a pas de second membre et raideur calculés
if (resu.res != NULL)
{ if (casMass_relax < 2) // 0 1 : on ne retient que la diagonale
{ Ass.AssembDiagonale (v_mass,*(resu.raid),elcontact.TableauDdlCont(),taN);}
else if ((casMass_relax > 1) && (casMass_relax < 6)) // sinon c'est la somme
// des valeurs absolues de la ligne
{ Ass.AssembDiagoMajorValPropre (v_mass,*(resu.raid),elcontact.TableauDdlCont(),taN);}
else
{ cout << "\n *** attention le cas du contact avec le casMass_relax = " << casMass_relax
<< " (c-a-d entre 6 a 9) "
<< " n'est pas encore pris en compte !! se plaindre !! "
<< " \n --> cas de la prise en compte de la raideur des elements de contact ";
Sortie(1);
};
// cas où on fait un deuxième assemblages
if ((casMass_relax > 3) && (casMass_relax < 6)) //4-5
{ Ass.AssembDiagonale (v_mass1,*(resu.raid),elcontact.TableauDdlCont(),taN);};
};
}; //-- fin du test sur l'activité
}; //-- fin de la boucle sur les éléments de contact
} //-- fin du try
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... )
{ if (ParaGlob::NiveauImpression() >= 1)
{cout << "\n warning: exception generee par un element de contact mais dont la prise en compte "
<< " n'est pas prevu !, on ne fait rien et on continue le calcul"
<< " \n --> cas de la prise en compte de la raideur des elements de contact ";
if (ParaGlob::NiveauImpression() >= 4) cout << "\n AlgoriRelaxDyna::CalculMatriceMasse(..";
};
};
}; // fin du test: if (pa.ContactType())
if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 5 )) || (permet_affichage > 2))
{// on affiche les min max du vecteur
cout << "\n -- info vecteur masse diagonale sur la raideur + contact -- ";
Coordonnee3 moy=v_mass.MinMaxMoy(true);
};
// --- dans le cas où on utilise la matrice de raideur réelle on calcul la diagonale
// condensée
if ((casMass_relax > 5) && (casMass_relax < 10))
{
if (casMass_relax > 7) // ici on a besoin de v_mass1
{for (int i=1;i<= nbddl_X; i++)
{v_mass(i)=0.;v_mass1(i)=0.;
Vecteur li = matglob->Ligne(i); // récup de la ligne
int taille_li = li.Taille();
for (int j=1;j<=taille_li;j++)
v_mass(i) += Dabs(li(j));
v_mass1(i) += 2. * Dabs((*matglob)(i,i));
};
}
else // on a besoin uniquement de v_mass
{for (int i=1;i<= nbddl_X; i++)
{v_mass(i)=0.;
Vecteur li = matglob->Ligne(i); // récup de la ligne
int taille_li = li.Taille();
for (int j=1;j<=taille_li;j++)
v_mass(i) += Dabs(li(j));
};
};
};
// maintenant pour certain cas, on va balayer sur la dimension (dima) les termes diagonaux de la masse
// et retenir le maxi selon les dima directions pour chaque noeud
// et ainsi construire la matrice masse correspondant à un pas de pseudo-temps critique de 1.
for (int inoeu = 1; inoeu <= nbNoe; inoeu++)
{ // on passe en revue les dima ddl, pour récupérer la raideur la plus élevée
int deb_indice = (inoeu-1)*dima; // deb_indice+1 = l'indice de x1
double raidmax = Dabs(v_mass(deb_indice+1)); // on récup la raideur (x1,x1)
for (int ix=2;ix<=dima;ix++) // on parcours les deux autres indices (si 3D), pour ne garder que la plus grande raideur
if (raidmax < Dabs(v_mass(deb_indice+ix)))
{raidmax = Dabs(v_mass(deb_indice+ix));};
// arrivée ici, raidmax contient le maxi selon les dima directions
// on calcule la matrice masse
if ((casMass_relax == 1) || (casMass_relax == 3) || (casMass_relax == 5)
|| (casMass_relax == 7) || (casMass_relax == 9)
)
{for (int jx=1;jx<=dima;jx++)
{(*mat_masse)(deb_indice+jx,deb_indice+jx) = 0.5 * lambda * raidmax;}
}
else // cas == 0 ou 2 ou 4 : dans ce cas on a une raideur différentes pour chaque direction
// == 6 ou 7
{ for (int jx=1;jx<=dima;jx++)
{(*mat_masse)(deb_indice+jx,deb_indice+jx) = 0.5 * lambda * v_mass(deb_indice+jx);}
};
};
// dans le cas casMass_relax == 4 ou 5, on utilise le maxi entre deux assemblages
// et casMass_relax == 8 ou 9
if ( (casMass_relax == 5) // pour le second second assemblage on récup le maxi sur les dima direction
|| (casMass_relax == 9)
)
{ for (int inoeu = 1; inoeu <= nbNoe; inoeu++)
{ // on passe en revue les dima ddl, pour récupérer la raideur la plus élevée
int deb_indice = (inoeu-1)*dima; // deb_indice+1 = l'indice de x1
double raidmax = v_mass1(deb_indice+1); // on récup la raideur (x1,x1)
for (int ix=2;ix<=dima;ix++) // on parcours les deux autres indices
// (si 3D), pour ne garder que la plus grande raideur
if (raidmax < v_mass1(deb_indice+ix))
{raidmax = v_mass1(deb_indice+ix);};
// arrivée ici, raidmax contient le maxi du second assemblage selon les dima directions
// maintenant pour la masse, on retiend la maxi entre les deux assemblages
for (int jx=1;jx<=dima;jx++)
(*mat_masse)(deb_indice+jx,deb_indice+jx) = DabsMaX((*mat_masse)(deb_indice+jx,deb_indice+jx)
, lambda * raidmax); // 1.*lambda et non 0.5, car on doit * 2
};
}
else if ((casMass_relax == 4) || (casMass_relax == 8) )
{ for (int i=1;i<=nbddl_X;i++) // ici on retiend simplement le maximum entre les deux assemblages
(*mat_masse)(i,i) = DabsMaX((*mat_masse)(i,i), lambda *v_mass1(i));
};
//v_mass.Affiche();
if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 5 )) || (permet_affichage > 2))
{// on affiche les min max du vecteur
cout << "\n -- info matrice masse finale -- ";
Coordonnee3 moy=mat_masse->MinMaxMoy(true);
};
break;
};
default:
{ cout << "\n **** erreur le cas type_calcul_mass= " << type_calcul_mass
<< "\n pour l'instant n'est pas implante "
<< "\n AlgoriRelaxDyna::CalculMatriceMasse( ... ";
Sortie(1);
};
};
// prise en compte de masses ponctuelles
// ---- cas des masses concentrées -----
Algori::Ajout_masses_ponctuelles(lesMail,Ass,(*mat_masse),divStoc,lesRef,N_ddl,lesFonctionsnD);
// ---- fin cas des masses concentrées -----
// maintenant on va transférer les valeurs aux noeuds, ce qui pourra être utile pour post-traiter
v_mass.Change_taille(nbddl_X); // au cas où car il n'a pas la même dim pour les différents cas de masse
for (int i=1;i<=nbddl_X;i++)
v_mass(i) = (*mat_masse)(i,i); // on utilise v_mass comme stockage intermédiaire
lesMail->Quelconque_glob_vers_local(X1,v_mass,masseRelax_2);
// si jamais on a une masse nulle on l'affiche éventuellement
if (ParaGlob::NiveauImpression() > 0)
for (int i=1;i<=nbddl_X;i++)
{if (Abs(v_mass(i)) < ConstMath::trespetit)
{// on récupère les informations
// retrouver le ddl correspondant a un pointeur de position
// d'assemblage, le nb du noeud et du maillage
// insol = le pointeur d'assemblage;
// ddl = le ddl en sortie; a t+dt si elle il existe
// sinon la valeur a t
// casAssemb : donne le cas d'assemblage qui est a considérer
int nbNoeud=0; int nbMaillage=0; // init
Ddl toto = lesMail->NoeudIndice(i,nbNoeud,nbMaillage
,Ass.Nb_cas_assemb());
int ne = (i-1)/dima+1;int coor = i - (ne-1)*dima;
cout << "\n *** attention : le noeud "<<nbNoeud
<< " du maillage "<< nbMaillage
<<" a une masse nulle en coordonnee " << coor ;
}
};
};
// calcul de la matrice masse dans le cas de l'algorithme de relaxation dynamique avec optimisation en continu
// de la matrice masse
// force_recalcul_masse : un indicateur pour forcer le calcul éventuellement
// en retour est mis à false s'il était true en entrée
void AlgoriRelaxDyna::CalculEnContinuMatriceMasse
(const int& relax_vit_acce,LesMaillages * lesMail,Assemblage& Ass
, int compteur, const DiversStockage* divStoc,LesReferences* lesRef
,const Enum_ddl & N_ddl
,bool premier_calcul,LesContacts* lescontacts
,bool & force_recalcul_masse,LesFonctions_nD* lesFonctionsnD)
{ // on regarde en fonction de option_recalcul_mass et des parametres associés éventuels s'il faut calculer la matrice masse
bool calcul_a_effectuer = false;
double epsilon = 0.;
// --- étude des différents cas de contrôle du recalcul de masse ---
int option_recalcul_mass_en_cours = option_recalcul_mass(1);
// init principal (convient s'il n'y a pas d'option spécifique visqueuse
if (fct_nD_option_recalcul_mass(1) != NULL)
{option_recalcul_mass_en_cours = (fct_nD_option_recalcul_mass(1)->Valeur_pour_variables_globales())(1);
if (option_recalcul_mass_en_cours != option_recalcul_mass(1))
{if (ParaGlob::NiveauImpression() > 2)
cout << "\n >>> changement : option_recalcul_mass "<< option_recalcul_mass(1) <<" en "
<< option_recalcul_mass_en_cours << std::flush;
option_recalcul_mass(1) = option_recalcul_mass_en_cours; // on sauvegarde par cohérence
};
};
// >> on regarde s'il y a plusieurs options
if (option_recalcul_mass.Taille() == 2)
{ // cas où on peut avoir du visqueux et/ou du cinétique
if (visqueux_activer)
{if (fct_nD_option_recalcul_mass(2) != NULL)
{option_recalcul_mass_en_cours = (fct_nD_option_recalcul_mass(2)->Valeur_pour_variables_globales())(1);
if (option_recalcul_mass_en_cours != option_recalcul_mass(2))
{if (ParaGlob::NiveauImpression() > 2)
cout << "\n >>> changement : option_recalcul_mass en visqueux "<< option_recalcul_mass(2) <<" en "
<< option_recalcul_mass_en_cours << std::flush;
option_recalcul_mass(2) = option_recalcul_mass_en_cours; // on sauvegarde par cohérence
}
}
else // sinon c'est la valeur fixe qu'il faut prendre
{option_recalcul_mass_en_cours = option_recalcul_mass(2);};
};
};
// --- fin étude des différents cas de contrôle du recalcul de masse ---
//cout << "\n option_recalcul_mass " << option_recalcul_mass_en_cours;
// on fait un premier passage pour voir si le calcul est a effectuer
switch (option_recalcul_mass_en_cours)
{case -1: // cas où la matrice est recalculée à chaque itération
{ calcul_a_effectuer = true;
break;
}
case 0: // cas où la matrice est calculée au début et ensuite elle est gardée constante
{ if (compteur == 0) // cas du premier calcul
calcul_a_effectuer = true;
break;
}
case 1: // apres un calcul au debut, mise a jour a chaque maxi de l'énergie cinétique
{ if ( (compteur == 0) // cas du premier calcul
|| (relax_vit_acce == 1)) // cas où on est à un maxi de l'énergie cinétique
calcul_a_effectuer = true;
break;
}
case 2: // idem le cas 1, mais on garde la valeur maxi de la raideur entre la nouvelle et l'ancienne
{ if ( (compteur == 0) // cas du premier calcul
|| (relax_vit_acce == 1)) // cas où on est à un maxi de l'énergie cinétique
calcul_a_effectuer = true;
break;
}
case 3: // recalcul apres ncycles, valeur qui est indiquée par ncycle_calcul
{ if ( (compteur == 0) // cas du premier calcul
|| ((compteur % ncycle_calcul)==0) ) // cas où on est au début de ncycle_calcul
calcul_a_effectuer = true;
break;
}
case 4: // recalcule de la matrice masse lorsque l'indicateur suivant "
// \epsilon = MAX_{i=1}^{nbddl}{ ( 0.5*lambda * (Delta t)^2) (|Delta ddot X_i|) / (|Delta X_i|) } "
// est superieur a 1*fac_epsilon "
{ int nbddl_X = delta_X.Taille();
epsilon = 0.;
for (int i=1; i<= nbddl_X;i++)
{ double denominateur = Dabs(delta_X(i));
if(denominateur > ConstMath::petit)
epsilon = MaX(epsilon,Dabs(acceleration_t(i))/denominateur);
};
epsilon *= 0.5 * lambda;
// on effectue le calcul conditionnel
if ( (compteur == 0) // cas du premier calcul
|| (epsilon > fac_epsilon) ) // cas ou la condition de recalcul est satisfaite
calcul_a_effectuer = true;
break;
}
case 5: // cas où on calcule la matrice au tout début du calcul et ensuite jamais
{ if ((compteur == 0) && premier_calcul) // cas du tout premier calcul
calcul_a_effectuer = true;
break;
}
default:
{ cout << "\n **** erreur le cas option_recalcul_mass= " << option_recalcul_mass_en_cours
<< "\n pour l'instant n'est pas implante "
<< "\n AlgoriRelaxDyna::CalculEnContinuMatriceMasse( ... ";
Sortie(1);
};
break;
};
// dans le cas ou on est à la première itération, le contact est updaté et
// on peut statuer s'il y a un pb avec des conditions linéaires de contact
if (((compteur <= 1 )&&(pa.ContactType()==1)))
{// mise à jour éventuelle du type et/ou de la taille de la matrice de masse
// retourne un pointeur sur la nouvelle matrice, s'il y a eu une modification
// l'ancienne est supprimée
// sinon retour d'un pointeur NULL
Mat_abstraite* mat_inter = Mise_a_jour_type_et_taille_matrice_masse_en_explicite
(delta_X.Taille(),mat_masse,lesMail,lesRef,(Ass.Nb_cas_assemb()).n,lescontacts);
if (mat_inter != NULL)
{ mat_masse = mat_inter; // c'est le nouveau stockage
// il faut dans ce cas également changer le stockage de la copie
delete (mat_masse_sauve);
mat_masse_sauve = mat_masse->NouvelElement(); // création d'un nouvel élément identique
calcul_a_effectuer=true;
}
else if (mat_masse->Place() != mat_masse_sauve->Place())
{// c'est le cas où la taille à diminuée sans changement de type de matrice
delete (mat_masse_sauve);
mat_masse_sauve = mat_masse->NouvelElement(); // création d'un nouvel élément identique
calcul_a_effectuer=true;
}
};
// on prend en compte un cas de forçage éventuel
if (force_recalcul_masse)
{calcul_a_effectuer = true;
force_recalcul_masse = false;
};
// on calcul la matrice masse si nécessaire
if (calcul_a_effectuer)
{ CalculMatriceMasse(relax_vit_acce,lesMail,Ass,compteur,divStoc
,lesRef,N_ddl,lescontacts,lesFonctionsnD);
// ici on va quand même éviter d'avoir des masses nulles sur la diagonale
// on commence par calculer la valeur moyenne
Coordonnee3 moy=mat_masse->MinMaxMoy((ParaGlob::NiveauImpression() > 3));
if ((ParaGlob::NiveauImpression() > 2) && (Dabs(moy(1))< ConstMath::petit))
{ cout << "\n warning la masse diagonale minimale est "<<moy(1)
<< " ??";
};
// on corrige la matrice masse avec le maxi ?
bool modif = false;
switch (choix_mini_masse_nul)
{case 1: modif = mat_masse->Limitation_min_diag(ConstMath::petit, moy(2)); break;
case 2: modif = mat_masse->Limitation_min_diag(ConstMath::petit, moy(3)); break;
default: // par défaut on ne fait rien
break;
};
// bool modif = mat_masse->Limitation_min_diag(ConstMath::petit, moy(1));
if (modif && (ParaGlob::NiveauImpression() > 2))
{ cout << "\n warning modification d'une masse diagonal initialement calculee a 0 "
<< " et mise a la valeur";
switch (choix_mini_masse_nul)
{case 1: cout << " maxi "; break;
case 2: cout << " moyenne "; break;
default: // par défaut on ne fait rien
break;
};
cout << " de la diag matrice = (min,max,moy= "; moy.Affiche_1(cout);
cout << ") "<< endl;
if (ParaGlob::NiveauImpression() >= 8)
{ string entete = " affichage de la matrice masse utilisee en relaxation dynamique";
cout << "\n " << entete;
mat_masse->Affiche();
cout << endl;
};
};
*mat_masse_sauve = *mat_masse;
}
else { *mat_masse = *mat_masse_sauve ;}; // on reprend l'ancienne sans les CL, si on n'a rien calculé
// info si voulu, dans le cas où option_recalcul_mass_en_cours ==-1 c'est un calcul à chaque itération donc on n'affiche rien
if ((calcul_a_effectuer && (ParaGlob::NiveauImpression() > 2)) && (option_recalcul_mass_en_cours != -1))
{ if (option_recalcul_mass_en_cours != 4)
{ cout << "\n recalcul de la pseudo-masse (type: "<< option_recalcul_mass_en_cours <<") "; }
else { cout << "\n recalcul de la pseudo-masse; epsilon=" << epsilon << endl ; };
};
// affichage éventuelle de la matrice de masse
if (ParaGlob::NiveauImpression() >= 10)
{ string entete = " affichage de la matrice masse utilisee en relaxation dynamique";
cout << "\n " << entete;
mat_masse->Affiche();
cout << endl;
};
};
//---- gestion des commndes interactives --------------
// écoute et prise en compte d'une commande interactive
// ramène true tant qu'il y a des commandes en cours
bool AlgoriRelaxDyna::ActionInteractiveAlgo()
{ cout << "\n commande? ";
return false;
};
// sortie du schemaXML: en fonction de enu
void AlgoriRelaxDyna::SchemaXML_Algori(ofstream& sort,const Enum_IO_XML enu) const
{
switch (enu)
{ case XML_TYPE_GLOBAUX :
{
break;
}
case XML_IO_POINT_INFO :
{
break;
}
case XML_IO_POINT_BI :
{
break;
}
case XML_IO_ELEMENT_FINI :
{
break;
}
case XML_ACTION_INTERACTIVE :
{sort << "\n <!-- ********** algorithme dynamique relaxation dynamique ******************** -->"
<< "\n<xs:complexType name=\"INIT\" >"
<< "\n <xs:annotation>"
<< "\n <xs:documentation> initialisation de l'algo "
<< "\n </xs:documentation>"
<< "\n </xs:annotation>"
<< "\n</xs:complexType>";
sort << "\n<xs:complexType name=\"EXECUTION\" >"
<< "\n <xs:annotation>"
<< "\n <xs:documentation> execution de l'ensemble de l'algo, sans l'initialisation et la derniere passe "
<< "\n </xs:documentation>"
<< "\n </xs:annotation>"
<< "\n</xs:complexType>";
sort << "\n<xs:complexType name=\"FIN_ALGO_EXPLI\" >"
<< "\n <xs:annotation>"
<< "\n <xs:documentation> fin de l'algo "
<< "\n </xs:documentation>"
<< "\n </xs:annotation>"
<< "\n</xs:complexType>";
break;
}
case XML_STRUCTURE_DONNEE :
{
break;
}
};
};
// dans le cas d'une relaxation de type 4 (mixte cinétique et visqueuse, on test pour savoir
// s'il faut du visqueux ou du cinétique
// ramène true si on continue en cinétique
// false si on bascule
// force_recalcul_masse : un indicateur de retour qui s'applique uniquement si
// et_recalcul_masse_a_la_transition = 1
// et que l'on est juste à la transition, sinon il est toujours false
bool AlgoriRelaxDyna::Cinetique_ou_visqueux(bool & force_recalcul_masse)
{ bool retour = true; // init en cinétique par défaut
force_recalcul_masse = false; // init par défaut
if (visqueux_activer)
// dans le cas ou la viscosité a été activé, on ne revient pas en arrière
{Transfert_ParaGlob_AMOR_CINET_VISQUEUX();
return false;
};
if (typeCalRelaxation==4)
{// dans le cas où il n'y a pas de fonction d'évolution, par défaut on utilise la précision
if ((niveauF_temps == NULL) && (niveauF_grandeurGlobale==NULL))
{ // on récupère le pointeur correspondant à NORME_CONVERGENCE
const void* pointe = (ParaGlob::param->GrandeurGlobal(NORME_CONVERGENCE));
TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
// le type est scalaire double
Grandeur_scalaire_double& gr
= *((Grandeur_scalaire_double*) gr_quelc->Grandeur_pointee()); // pour simplifier
double norme_actuelle = *(gr.ConteneurDouble());
double precision_seuil = Algori::Precision_equilibre(); // initialisation de la précision demandé
if (norme_actuelle * proportion_cinetique > precision_seuil)
retour = true; // on continue en cinétique
else
retour = false; // on demande la bascule
}
else // sinon cela signifie qu'au moins une fonction est définie
{// cas d'une dépendance au temps
if (niveauF_temps != NULL)
{ // on récupère le temps
const VariablesTemps& tempo = pa.Variables_de_temps();
double temps = tempo.TempsCourant();
// on filtre avec retour et la version int de la fonction
retour = retour && (int (niveauF_temps->C_proport()->Valeur(temps)) );
};
// cas d'une dépendance à des variables globales
if (niveauF_grandeurGlobale != NULL)
{ // on appel la fonction
Tableau <double> & tab_val = niveauF_grandeurGlobale->C_proport()->Valeur_pour_variables_globales();
// on n'utilise que la première valeur
// on filtre avec retour et la version int de la fonction
retour = retour && (int (tab_val(1)));
};
};
};
// if (ParaGlob::NiveauImpression()>2 && !retour)
// cout << "\n ==> relaxation dynamique: bascule de cinetique a visqueux " << endl;
if (retour == false)
{visqueux_activer = true; // init pour le prochain passage
if (et_recalcul_masse_a_la_transition)
force_recalcul_masse = true; // on indique que l'on vient juste de changer
if (ParaGlob::NiveauImpression() >= 1)
{ cout << "\n\n $$$$$$$$$ bascule vers amortissement visqueux $$$$$$$$ "
<< "\n $$$$$$$$$ (iteration: "<<compteur << " ) $$$$$$$$ \n" << endl ;
};
}
Transfert_ParaGlob_AMOR_CINET_VISQUEUX();
return retour;
};
// calcul du facteur lambda en fonction de ce qui est demandé (automatique ou via une fonction etc.)
double AlgoriRelaxDyna::Calcul_Lambda()
{ double retour = lambda; // init via la valeur par défaut
double retour_initial = lambda; // pour voir si on a réellement changé
bool a_changer=false; // indicateur pour savoir si on doit augmenter
// --- on regarde s'il y a une dépendance à une fonction d'évolution
if ( (niveauLambda_temps != NULL) || (niveauLambda_grandeurGlobale != NULL))
{// on a une fonction d'évolution donc on n'utilise pas le lambda_initiale
// retour = 1.; // init pour les multiplications
// dans ce cas on initialise avec la valeur initiale de lambda
retour = lambda_initial;
// on regarde s'il y a une fonction d'évolution
if (niveauLambda_temps != NULL)
{ // on récupère le temps
const VariablesTemps& tempo = pa.Variables_de_temps();
double temps = tempo.TempsCourant();
retour *= (niveauLambda_temps->C_proport()->Valeur(temps));
a_changer=true;
};
// cas d'une dépendance à des variables globales
if (niveauLambda_grandeurGlobale != NULL)
{ // on appel la fonction
Tableau <double> & tab_val = niveauLambda_grandeurGlobale->C_proport()->Valeur_pour_variables_globales();
double facteur = tab_val(1); // pour simplifier
if (facteur != 1.) // cela veut qu'il y a du changement voulu
{ a_changer=true;
retour *= facteur;
};
};
};
// --- on regarde s'il y a une gestion automatique de lambda
if (pilotage_auto_lambda)
{int iter_maxi = 25; // une limite arbitraire pour l'instant figé
bool a_augmenter=false; // indicateur pour savoir si on doit augmenter
// on va parcourir les itérations de relaxation
List_io<int>::iterator ili,ilifin=list_iter_relax.end();
// l'idée est d'explorer sur 5*5 iter = 25,
// si on a 1 relax tous les 2 iter, 3 fois à suivre, cela veut dire que l'on a 3 relax sur 6 iter
// et donc 3 enregistrements tel que le (3ième - le premier) < 5 par exemple
// ex de suite: 14, 12, 10 --> 14-10 = 4 < 7
// si on a moins de 3 relax on sort
// idem pour 1 relax tous les 3 iter, 3 fois à suivre
// ex: 19, 16, 13 --> 19-13 = 6 < 10
// idem pour 1 relax tous les 3 iter, 4 fois à suivre
// ex: 17, 14, 11, 8 --> 17-8 = 9 < 21
if (list_iter_relax.size() > 3)
{ili=list_iter_relax.begin();
int incre_deb = *ili; // c'est le plus grand car on empile par l'avant
int compt_1R_2iter = 1;
for (;ili != ilifin; ili++,compt_1R_2iter++)
{ if (compt_1R_2iter == 3)
{if ((incre_deb - (*ili)) < 7)
{// on doit augmenter lambda
retour += delta_lambda;
// il faut que l'on supprime tous les éléments qui sont sup et = à iter
// peut-être pas suffisant list_iter_relax.erase(ili,ilifin);
a_augmenter = a_changer = true;
//--- debug
//cout << "\n debug AlgoriRelaxDyna::Calcul_Lambda() " << lambda << " --> " << retour ;
//--- fin debug
};
};
if (compt_1R_2iter == 4)
{if ((incre_deb - (*ili)) < 10)
{// on doit augmenter lambda
retour += delta_lambda;
// il faut que l'on supprime tous les éléments qui sont sup et = à iter
// peut-être pas suffisant list_iter_relax.erase(ili,ilifin);
a_augmenter = a_changer = true;
//--- debug
//cout << "\n debug AlgoriRelaxDyna::Calcul_Lambda() " << lambda << " --> " << retour ;
//--- fin debug
};
};
if (compt_1R_2iter == 5)
{if ((incre_deb - (*ili)) < 21)
{// on doit augmenter lambda
retour += delta_lambda;
// il faut que l'on supprime tous les éléments qui sont sup et = à iter
// peut-être pas suffisant list_iter_relax.erase(ili,ilifin);
a_augmenter = a_changer = true;
//--- debug
//cout << "\n debug AlgoriRelaxDyna::Calcul_Lambda() " << lambda << " --> " << retour ;
//--- fin debug
};
};
// si on a augmenté,
if (a_augmenter)
{ list_iter_relax.clear(); // en test pour l'instant car c'est radical
break;
}
};
// si on a plus d'incrément que le nombre limite, on supprime les 3 plus vieux
if (list_iter_relax.size() > iter_maxi)
for (int j= 1; j<=3;j++)
{ili=list_iter_relax.end(); ili--;
list_iter_relax.erase(ili);
};
}
else if (list_iter_relax.size() == 2)
// sinon cas de 2 éléments, on supprime éventuellement les éléments s'ils sont trop éloignés
{ ili=list_iter_relax.begin();
int incre_deb = *ili;
ili++;
if ((incre_deb - (*ili)) > iter_maxi)
list_iter_relax.erase(ili);
};
// on regarde les mini et maxi
if (retour > lambda_max) retour = lambda_max;
if (retour < lambda_min) retour = lambda_min;
};
// si on a augmenté,
if ((a_changer)&&(retour != retour_initial))
{ if (ParaGlob::NiveauImpression() >= 1)
{ cout << "\n modification de lambda " << lambda << " --> " << retour
<<" (delta= "<< (retour-retour_initial)<<") " ;
};
};
// retour
return retour;
};
// fonction virtuelle permettant de mettre à jour les infos aux noeuds
// à cause de certains contact (ex: cas_contact = 4)
// Les vecteurs Xtdt et Vtdt sont mis à jour par la méthode
// la vitesse locale du noeud est mis à jour en fonction de la position locale et de l'algo
void AlgoriRelaxDyna::Repercussion_algo_sur_cinematique(LesContacts* lesContacts,Vecteur& Xtdt,Vecteur& Vtdt)
{ // on récupère la liste des noeuds éventuellement modifiés par le contact
list <Noeud * > li_noe; // init
lesContacts->Liste_noeuds_position_changer(li_noe);
// maintenant on parcours la liste (qui peut-être vide )
list <Noeud * >::iterator il,ilfin = li_noe.end();
for (il = li_noe.begin();il != ilfin; il++)
{Noeud & noe = *(*il); // pour simplifier
const Coordonnee& Mtdt = noe.Coord2(); // récup des positions
const Coordonnee& Mt = noe.Coord1();
//cout << "\n Mtdt: "<< Mtdt << ", Mt: "<< Mt;
int dim = Mtdt.Dimension();
int nb_assemb = (Ass1_->Nb_cas_assemb()).n; // cas d'assemblage pour les X mais
// qui correspond à la même position
// dans le vecteur global que les vitesses
// on met à jour localement les vitesse, et globalement positions et vitesses
switch (dim)
{ case 3:if (!(noe.Ddl_fixe(X3)))
{ double v3 = Mtdt(3) - Mt(3);
noe.Change_val_tdt(V3,v3);
int posi = noe.Pointeur_assemblage(X3,nb_assemb);
Xtdt(posi) = Mtdt(3); Vtdt(posi) = v3;
}
case 2:if (!(noe.Ddl_fixe(X2)))
{ double v2 = Mtdt(2) - Mt(2);
noe.Change_val_tdt(V2,v2);
int posi = noe.Pointeur_assemblage(X2,nb_assemb);
//cout << "\n posi= "<<posi<< " Mtdt(2)="<<Mtdt(2);
Xtdt(posi) = Mtdt(2); Vtdt(posi) = v2;
}
case 1:if (!(noe.Ddl_fixe(X1)))
{ double v1 = Mtdt(1) - Mt(1);
noe.Change_val_tdt(V1,v1);
int posi = noe.Pointeur_assemblage(X1,nb_assemb);
Xtdt(posi) = Mtdt(1); Vtdt(posi) = v1;
}
default:
break;
};
};
};
Reference in a new issue Copy permalink