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

1707 lines
99 KiB
C++
Raw Permalink 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_tchamwa.h"
// CONSTRUCTEURS :
AlgoriTchamwa::AlgoriTchamwa () : // par defaut
Algori()
,alphaa(0.5),gammaa(0.5),lambdaa(1.),betaa(0.)
,delta_t(0.),deltat2(0.),unsurdeltat(0.),lambdaDeltat(0.),alphaDelta_t(0.),betaDeltat2(0.),gammaDelta_t(0.)
,delta_t_critique(0.)
,type_cal_equilibre(1),CGamma_pourVarPhi(NULL),nom_courbe_CGamma_pourVarPhi()
,degre_num(1),degre_deno(1) // utilisation des approximants de pade
,type4_inc_deb(),type4_inc_deplace() // exploratoire
,npas_moyacc(0),npas_effectue(0),valmin(0.),valmax(ConstMath::tresgrand),moy_acc(),moy_acc_en_calcul() // exploratoire
// --------------------------------------------------------------------------------------
// -- variables de transferts internes entre: InitAlgorithme, CalEquilibre, FinCalcul --
// --------------------------------------------------------------------------------------
,Ass1_(NULL),Ass2_(NULL),Ass3_(NULL)
,cas_combi_ddl(),icas(),prepa_avec_remont(false)
,brestart(false),type_incre(OrdreVisu::PREMIER_INCRE)
,vglobin(),vglobex(),vglobaal(),vcontact()
,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)
// ------------------------------------------------------------------------------------------
// -- fin variables de transferts internes entre: InitAlgorithme, CalEquilibre, FinCalcul --
// ------------------------------------------------------------------------------------------
{ //vglobaal = &vglobin;
// message d'erreur
cout << "\n constructeur par defaut de AlgoriTchamwa, 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
AlgoriTchamwa::AlgoriTchamwa (const bool avec_typeDeCal
,const list <EnumSousTypeCalcul>& soustype
,const list <bool>& avec_soustypeDeCal
,UtilLecture& entreePrinc) :
Algori(DYNA_EXP_TCHAMWA,avec_typeDeCal,soustype,avec_soustypeDeCal,entreePrinc)
,alphaa(0.5),gammaa(0.5),lambdaa(1.),betaa(0.)
,delta_t(0.),deltat2(0.),unsurdeltat(0.),lambdaDeltat(0.),alphaDelta_t(0.),betaDeltat2(0.),gammaDelta_t(0.)
,delta_t_critique(0.)
,type_cal_equilibre(1),CGamma_pourVarPhi(NULL),nom_courbe_CGamma_pourVarPhi()
,degre_num(1),degre_deno(1) // utilisation des approximants de pade
,type4_inc_deb(),type4_inc_deplace() // exploratoire
,npas_moyacc(0),npas_effectue(0),valmin(0.),valmax(ConstMath::tresgrand),moy_acc(),moy_acc_en_calcul() // exploratoire
// --------------------------------------------------------------------------------------
// -- variables de transferts internes entre: InitAlgorithme, CalEquilibre, FinCalcul --
// --------------------------------------------------------------------------------------
,Ass1_(NULL),Ass2_(NULL),Ass3_(NULL)
,cas_combi_ddl(),icas(),prepa_avec_remont(false)
,brestart(false),type_incre(OrdreVisu::PREMIER_INCRE)
,vglobin(),vglobex(),vglobaal(),vcontact()
,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)
// ------------------------------------------------------------------------------------------
// -- fin variables de transferts internes entre: InitAlgorithme, CalEquilibre, FinCalcul --
// ------------------------------------------------------------------------------------------
{ //vglobaal = &vglobin;
// lecture des paramètres attachés au type de calcul
switch (entreePrinc.Lec_ent_info())
{ case 0 :
{ lecture_Parametres(entreePrinc);
betaa = (*phi_1) - gammaa * lambdaa;
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
AlgoriTchamwa::AlgoriTchamwa (const AlgoriTchamwa& algo):
Algori(algo),alphaa(0.5),gammaa(0.5),lambdaa(1.),betaa(algo.betaa)
,delta_t(0.),deltat2(0.),unsurdeltat(0.),lambdaDeltat(0.),alphaDelta_t(0.),betaDeltat2(0.),gammaDelta_t(0.)
,delta_t_critique(0.)
,type_cal_equilibre(algo.type_cal_equilibre),CGamma_pourVarPhi(algo.CGamma_pourVarPhi)
,nom_courbe_CGamma_pourVarPhi(algo.nom_courbe_CGamma_pourVarPhi)
,degre_num(algo.degre_num),degre_deno(algo.degre_deno) // utilisation des approximants de pade
,type4_inc_deb(),type4_inc_deplace() // exploratoire
,npas_moyacc(0),npas_effectue(0),valmin(0.),valmax(ConstMath::tresgrand),moy_acc(),moy_acc_en_calcul() // exploratoire
// --------------------------------------------------------------------------------------
// -- variables de transferts internes entre: InitAlgorithme, CalEquilibre, FinCalcul --
// --------------------------------------------------------------------------------------
,Ass1_(NULL),Ass2_(NULL),Ass3_(NULL)
,cas_combi_ddl(),icas()
,prepa_avec_remont(false)
,brestart(false),type_incre(OrdreVisu::PREMIER_INCRE)
,vglobin(),vglobex(),vglobaal(),vcontact()
,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)
// ------------------------------------------------------------------------------------------
// -- fin variables de transferts internes entre: InitAlgorithme, CalEquilibre, FinCalcul --
// ------------------------------------------------------------------------------------------
{ //vglobaal = &vglobin;
phi_1 = & paraTypeCalcul(1);
// on regarde s'il s'agit d'une courbe locale ou d'une courbe globale pour CGamma_pourVarPhi
if (CGamma_pourVarPhi->NomCourbe() == "_")
{// comme il s'agit d'une courbe locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
// CGamma_pourVarPhi = Courbe1D::New_Courbe1D(non_courbe,algo.CGamma_pourVarPhi->Type_courbe());
CGamma_pourVarPhi = Courbe1D::New_Courbe1D(*algo.CGamma_pourVarPhi);
};
};
// destructeur
AlgoriTchamwa::~AlgoriTchamwa ()
{ 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 (CGamma_pourVarPhi != NULL)
if (CGamma_pourVarPhi->NomCourbe() == "_") delete CGamma_pourVarPhi;
};
// execution de l'algorithme dans le cas dynamique explicite,
void AlgoriTchamwa::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
{ switch (type_cal_equilibre)
{ case 4:
// 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_Equilibre4(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
break;
case 2:
// 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_Equilibre2(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
break;
case 1: case 3: // dernier né: permet de faire plusieurs appels,
{// 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 );
};
break;
}
};
}
// { if ((type_cal_equilibre==1) || (type_cal_equilibre==3))
// {Calcul_Equilibre(paraGlob,lesMail,lesRef,lesCourbes1D,lesLoisDeComp
// ,toto,charge,lesCondLim,titi,tutu );
// }
// else if (type_cal_equilibre==2)
// {Calcul_Equilibre2(paraGlob,lesMail,lesRef,lesCourbes1D,lesLoisDeComp
// ,toto,charge,lesCondLim,titi,tutu );
// }
// else if (type_cal_equilibre==4)
// {Calcul_Equilibre4(paraGlob,lesMail,lesRef,lesCourbes1D,lesLoisDeComp
// ,toto,charge,lesCondLim,titi,tutu );
// };
// }
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
{// 1- on fait le calcul d'équilibre
switch (type_cal_equilibre)
{ case 4:
// 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_Equilibre4(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
break;
case 2:
// 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_Equilibre2(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,divStock,charge,lesCondLim,lesContacts,resultats );
break;
case 1: case 3: // dernier né: permet de faire plusieurs appels,
{// 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 );
};
break;
}
};
// if ((type_cal_equilibre==1) || (type_cal_equilibre==3))
// {Calcul_Equilibre(paraGlob,lesMail,lesRef,lesCourbes1D,lesLoisDeComp
// ,toto,charge,lesCondLim,titi,tutu );
// }
// else if (type_cal_equilibre==2)
// {Calcul_Equilibre2(paraGlob,lesMail,lesRef,lesCourbes1D,lesLoisDeComp
// ,toto,charge,lesCondLim,titi,tutu );
// }
// else if (type_cal_equilibre==4)
// {Calcul_Equilibre4(paraGlob,lesMail,lesRef,lesCourbes1D,lesLoisDeComp
// ,toto,charge,lesCondLim,titi,tutu );
// };
// 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);
};
// lecture des paramètres du calcul
void AlgoriTchamwa::lecture_Parametres(UtilLecture& entreePrinc)
{ // dimensionnement
paraTypeCalcul.Change_taille(1);
MotCle motCle; // ref aux mots cle
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;
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 du premier paramètres de l'algorithme
entreePrinc.NouvelleDonnee(); // ligne suivante
// lecture du nom du paramètre et vérification
string st1;
if (strstr(entreePrinc.tablcar,"phi=")!=NULL)
{ *(entreePrinc.entree) >> st1 ;
if (st1 != "phi=")
{ cout << "\n erreur en lecture du parametre de l'algorithme de tchamwa "
<< "\n on attendait le mot : phi= , au lieu de " << st1
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
}
// lecture du parametre
*(entreePrinc.entree) >> paraTypeCalcul(1);
lecture_effective = true;
}
else // sinon paramètre par défaut
{paraTypeCalcul(1) = 1.01;};
// cas où le paramètre type_cal_equilibre existe
if (strstr(entreePrinc.tablcar,"typeCalEqui=")!=NULL)
{ // on lit le type de calcul d'équilibre
string st1;
*(entreePrinc.entree) >> st1 >> type_cal_equilibre;
if (st1 != "typeCalEqui=")
{ cout << "\n erreur en lecture du type de calcul d'equilibre "
<< "\n on attendait le mot : typeCalEqui= , au lieu de " << st1
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
// on vérifie la valeur lue
if ((type_cal_equilibre<1)||(type_cal_equilibre>4))
// if ((type_cal_equilibre!=1)&&(type_cal_equilibre!=2)&&(type_cal_equilibre!=3))
{ cout << "\n erreur en lecture du type de calcul d'equilibre "
<< "\n on attendait une valeur comprise entre 1 et 4 au lieu de: " << type_cal_equilibre
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
if (type_cal_equilibre == 2)
{// on lit la courbe de modération de phi en fonction de Gamma
*(entreePrinc.entree) >> st1 >> nom_courbe_CGamma_pourVarPhi;
if (st1 != "CGamma_pourVarPhi=")
{ cout << "\n erreur en lecture du mot cle pour la courbe de moderation "
<< "\n on attendait le mot : CGamma_pourVarPhi= , au lieu de " << st1
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
// on regarde s'il faut utiliser une moyenne d'accélération
if (strstr(entreePrinc.tablcar,"n_=")!=NULL)
{// on lit le nombre de pas
*(entreePrinc.entree) >> st1 >> npas_moyacc;
if (st1 != "n_=")
{ cout << "\n erreur en lecture du mot cle pour le nombre de pas "
<< "\n on attendait le mot : n_= , au lieu de " << st1
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
// on lit les valeurs mini et maxi
string st2;
*(entreePrinc.entree) >> st1 >> valmin >> st2 >> valmax;
if ((st1 != "valmin_=") || (st2 !="valmax_=" ))
{ cout << "\n erreur en lecture des mots cle pour le mini et maxi de la moyenne "
<< "\n on attendait les mots : valmin_= et valmax_=, au lieu de " << st1 << " " << st2
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
// on vérifie que les grandeurs lues sont positives
if ((valmax < 0.) || (valmin < 0.))
{ cout << "\n erreur en lecture: les grandeurs valmax et valmin doivent etre positives "
<< "\n alors que l'on a lue valmax= " << valmax << " valmin= " << valmin
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
};
}
else if (type_cal_equilibre == 3) // cas de l'utilisation des approximants de pade
{// on lie les degres du numerateur et du denominateur et on verifie
*(entreePrinc.entree) >> st1 >> degre_num;
if (st1 != "degre_num_")
{ cout << "\n erreur en lecture du mot cle degre_num_ on a lu " << st1
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
*(entreePrinc.entree) >> st1 >> degre_deno;
if (st1 != "degre_deno_")
{ cout << "\n erreur en lecture du mot cle degre_deno_ on a lu " << st1
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
// maintenant on vérifie que les types sont exploitables
if ((degre_num != 1) && ((degre_deno<1) || (degre_deno>2)))
{ cout << "\n erreur seul les degres de pade 1 - 1 ou 1 - 2 sont actuellement implante " << st1
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
}
else if (type_cal_equilibre == 4)
{// on lie la courbe de modération de phi en fonction de Gamma
*(entreePrinc.entree) >> st1 >> type4_inc_deb;
if (st1 != "type4_inc_deb_")
{ cout << "\n erreur en lecture du mot cle type4_inc_deb_ on a lu " << st1
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
*(entreePrinc.entree) >> st1 >> type4_inc_deplace;
if (st1 != "type4_inc_deplace_")
{ cout << "\n erreur en lecture du mot cle type4_inc_deplace_ on a lu " << st1
<< "\n AlgoriTchamwa::lecture_Parametres( ... ";
Sortie(1);
};
};
lecture_effective=true;
}
else // sinon paramètre par défaut
{type_cal_equilibre = 1;};
}
else // sinon on met une valeur par défaut qui ici correspond au cas
// des différences centrées pour l'accélération, mais décalé à gauche
// pour la vitesse (ce qui est aussi le cas quand phi est différent de 1.)
{ paraTypeCalcul(1) = 1.01;
type_cal_equilibre = 1;
}
// on relie le paramètre à phie
phi_1 = & paraTypeCalcul(1);
// 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);
};
// écriture des paramètres dans la base info
// = 1 : on écrit tout
// = 2 : on écrit uniquement les données variables (supposées comme telles)
void AlgoriTchamwa::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());
if (cas == 1)
{
// ecriture du parametre
*(sort) << " phi= " << paraTypeCalcul(1) << " " << "typeCalEqui= " << type_cal_equilibre << " ";
if (type_cal_equilibre == 2)
{*(sort) << " CGamma_pourVarPhi " << nom_courbe_CGamma_pourVarPhi;
};
};
// (*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 : fonctionnememt 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 AlgoriTchamwa::Lecture_Base_info_Parametre(UtilLecture& entreePrinc,const int& cas,bool choix)
{if (cas == 1)
{// dimensionnement
paraTypeCalcul.Change_taille(1);
if (choix)
{// cas d'une lecture normale
// récup du flot
ifstream * ent = entreePrinc.Ent_BI();
string toto;
// lecture du parametre
*(ent) >> toto ;
if (toto != "phi=")
{ cout << "\n erreur en lecture du parametre de l'algorithme explicite: tchamwa"
<< "\n on attendait le mot : phi= , au lieu de " << toto
<< "\n AlgoriTchamwa::Lecture_Base_info_Parametre( ... ";
Sortie(1);
}
*(ent) >> paraTypeCalcul(1) ;
// lecture du type d'équilibre
*(ent) >> toto ;
if (toto != "typeCalEqui=")
{ cout << "\n erreur en lecture du parametre typeCalEqui: tchamwa"
<< "\n on attendait le mot : typeCalEqui= , au lieu de " << toto
<< "\n AlgoriTchamwa::Lecture_Base_info_Parametre( ... ";
Sortie(1);
}
*(ent) >> type_cal_equilibre ;
if (type_cal_equilibre == 2)
{ *(ent) >> toto;
if (toto != " CGamma_pourVarPhi ")
{ cout << "\n erreur en lecture de la fonction CGamma_pourVarPhi,"
<< " on attendait CGamma_pourVarPhi et on a lue " << toto
<< "\n AlgoriTchamwa::Lecture_base_info_loi(...";
Sortie(1);
};
*(ent) >> nom_courbe_CGamma_pourVarPhi;
};
}
else
{// cas où la lecture n'est pas possible, attribution des valeurs par défaut
paraTypeCalcul(1) = 1.01;
type_cal_equilibre = 1;
}
// on relie le paramètre à phie
phi_1 = & paraTypeCalcul(1);
}
};
// création d'un fichier de commande: cas des paramètres spécifiques
void AlgoriTchamwa::Info_commande_parametres(UtilLecture& entreePrinc)
{ // écriture dans le fichier de commande
ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
sort << "\n#-----------------------------------------------------------------------------"
<< "\n| parametres (falcultatifs ) associes au calcul de dynamique explicite |"
<< "\n| avec le modele de Tchamwa - Wielgoz: (amortissement des freq numeriques) |"
<< "\n| phi (facteur d'attenuation) un peu sup à 1 (ex 1.01) |"
<< "\n#-----------------------------------------------------------------------------"
<< "\n"
<< "\n PARA_TYPE_DE_CALCUL"
<< "\n # ........................................................."
<< "\n # / facteur d'attenuation : phi = 1.01 /"
<< "\n #........................................................"
<< "\n phi= 1.01"
<< "\n # <<<<<< la suite concerne des versions experimentales >>>>>>>>> "
<< "\n # on peut egalement mettre sur la meme ligne le type d'algorithme "
<< "\n # par defaut =1, dans le cas du type 2 (version recherche), il faut egalement indiquer le nom de la "
<< "\n # courbe de modulation de phi en fonction de la norme de l'acceleration "
<< "\n # exemple: phi= 1.01 typeCalEqui= 2 CGamma_pourVarPhi= fonc "
<< "\n # un rafinement est possible en se referent a une moyenne des accelerations sur "
<< "\n # n pas. (moy(i) = (somme (pour r variant de 0 a n) delta gamma_temps_r(i)/ * delta t)/n)."
<< "\n # Le coefficient phi est pondere suivant la valeur de cette moyenne, entre 1 et la "
<< "\n # valeur demandee. (|moy(i)| > valmax ) ==> 1 , (|moy(i)| < valmin) ==> phi, entre valmin et valmax "
<< "\n # la progression est lineaire. En appelant phi_modif(i), la valeur de ce phi modifie, le coefficient phi_utilise "
<< "\n # retenue a pour valeur au noeud 'i' : "
<< "\n # phi_utilise(i) = (1. + (phi_modif(i)-1) * fonc(|gamma(i)|) "
<< "\n # Remarque: le phi_modif est mis a jour uniquement tous les n pas"
<< "\n # n par defaut = 0 et phi_modi=phi "
<< "\n # exemple d'utilisation de la valeur de la moyenne : "
<< "\n # "
<< "\n # phi= 1.03 typeCalEqui= 2 CGamma_pourVarPhi= fonc n_= 4 valmin_= 1.e1 valmax_= 1.e3 "
<< "\n # "
<< "\n # dans le cas ou le type = 3, on utilise une technique de pade (experimental) "
<< "\n # dans ce cas il faut indiquer deux parametres degre_num_ et degre_deno_ "
<< "\n # degre_num_ indique le degre du polynome numerateur, degre_deno_ celui du denominateur "
<< "\n # actuellement sont implantes: degre 1 - 1 et 1 - 2 "
<< "\n # dans le cas ou le type = 4, cas d'une structure 1D elastique, et une fenetre qui se deplace "
<< "\n # sur le maillage, a une vitesse definit par l'utilisateur) dans ce cas 2 parametres"
<< "\n # type4_inc_deb_ : nb increment a partir duquel la fenetre avance "
<< "\n # type4_inc_deplace_ : nb increment tous les combiens on deplace la fenetre d'un noeud ";
// appel de la classe mère
Algori::Info_com_parametres(entreePrinc);
sort << "\n" << endl;
};
// gestion et vérification du pas de temps et modif en conséquence si nécessaire
// 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
// cas = 3: il y a une demande de changement de pas de temps: nouveau_dt
// si nouveau_dt = 0, on ne fait rien
// 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 AlgoriTchamwa::Gestion_pas_de_temps(bool modif_pas_de_temps,LesMaillages * lesMail,int cas, double nouveau_dt)
{ 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
switch (cas)
{ case 1 :
{ // --<T W>-- récup du pas de temps, proposé par l'utilisateur
delta_t = pa.Deltat(); double delta_t_old = delta_t;
// --<T W>-- 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
double l_sur_c = lesMail->Longueur_arrete_mini_sur_c(TEMPS_0);
double delta_t_essai = l_sur_c / (*phi_1);
if (permet_affichage > 2)
cout << "\n pas de temps critique DFC: "<<l_sur_c
<< " pas de temps critique tchamwa: " << delta_t_essai << flush;
delta_t_critique = delta_t_essai;
// mise à jour éventuel du pas de temps et du pas de temps maxi 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
// maintenant on vérifie le pas de temps choisit
if (pa.Limit_temps_stable() && ( delta_t > delta_t_essai))
{ if ((ParaGlob::NiveauImpression() > 0)|| (permet_affichage > 0))
cout << "\n ** ATTENTION ** , l'increment de temps propose : " << delta_t
<< ", ne satisfait pas la condition de Courant (C-F-L), "
<< "\n on le regle donc pour que cette condition soit satisfaite: nouvel increment: "
<< delta_t_essai << endl;
delta_t = delta_t_essai;
if (delta_t_old != delta_t)
modif_deltat=true;
// modification de l'increment de temps dans les paramètres
switch (pa.Modif_Deltat(delta_t))
{case -1: { if ((ParaGlob::NiveauImpression() > 0)|| (permet_affichage > 0))
cout << "\n initialisation du pas de temps avec la condition de Courant impossible, le nouveau pas"
<< "\n de temps calcule est plus petit que le pas de temps mini limite: " << pa.Deltatmini();
Sortie(1); break;
}
case 0: { break;} // cas normal
case 1: { pa.Modif_Detat_dans_borne(delta_t);
if ((ParaGlob::NiveauImpression() > 0)|| (permet_affichage > 0))
cout << "\n >>>> rectification du pas de temps "<<delta_t_old<<" dans les bornes : nouvel increment:" << delta_t;
modif_deltat=true;break;
} // cas ou on dépasse la borne maxi
};
}
else
// sinon on regarde simplement si la modification de temps est possible
{ switch (pa.Modif_Deltat(delta_t))
{case -1: { cout << "\n initialisation du pas de temps impossible, le nouveau pas"
<< "\n de temps calcule est plus petit que le pas de temps mini limite: " << pa.Deltatmini();
Sortie(1); break;
}
case 0: { break;} // cas normal
case 1: { pa.Modif_Detat_dans_borne(delta_t);
if ((ParaGlob::NiveauImpression() > 0)|| (permet_affichage > 0))
cout << "\n >>>> rectification du pas de temps "<<delta_t_old<<" dans les bornes : nouvel increment:" << delta_t;
modif_deltat=true;break;
} // cas ou on dépasse la borne maxi
};
};
break;
}
case 2:
// on regarde si le pas de temps n'est pas devenu supérieur au pas critique
// ou au contraire que le pas de temps critique a augmenté
{ double l_sur_c = lesMail->Longueur_arrete_mini_sur_c(TEMPS_t);
double delta_t_essai = l_sur_c / (*phi_1) ;double ancien_pas=delta_t;
if (permet_affichage > 2)
cout << "\n pas de temps critique DFC: "<<l_sur_c
<< " pas de temps critique tchamwa: " << delta_t_essai << flush;
//on intervient: si le nouveau pas de temps critique est inférieur au pas de temps en cours
// ou si le nouveau pas de temps critique de l'algorithme est supérieur à l'ancien
if (( delta_t_essai < delta_t) || ( delta_t_essai > delta_t_critique))
{// mise à jour éventuel du pas de temps et du pas de temps maxi s'ils sont définit à partir du temps critique
bool modif = pa.Modif_Deltat_DeltatMaxi(delta_t_essai,l_sur_c); // mais pas de test vis-a-vis des bornes
// s'il y a eu modif du pas de temps on met à jour le pas courant
if (modif) {
delta_t = pa.Deltat(); modif_deltat=true;}
// maintenant on vérifie néanmoins le pas de temps choisit au cas où
// en particulier cette vérification n'est en fait pas utile si l'utilisateur
// à utilisé un facteur du pas critique < 1
if (pa.Limit_temps_stable() && ( delta_t > delta_t_essai))
{ if (pa.Coefficient_pas_critique_deltat() <= 1.)
{ delta_t = pa.Coefficient_pas_critique_deltat() * delta_t_essai;}
else { delta_t = delta_t_essai;};
// modification de l'increment de temps dans les paramètres
switch (pa.Modif_Deltat(delta_t))
{case -1: { if ((ParaGlob::NiveauImpression() > 0)|| (permet_affichage > 0))
cout << "\n **** modification du pas de temps impossible,"
<< " le nouveau pas de temps calcule est plus petit que le pas de temps mini limite: "
<< pa.Deltatmini();
Sortie(1); break;
}
case 0: { break;} // cas normal
case 1: { pa.Modif_Detat_dans_borne(delta_t); break;} // cas ou on dépasse la borne maxi
};
modif_deltat=true;
}
else // sinon on regarde simplement si la modification de temps est possible
{ switch (pa.Modif_Deltat(delta_t))
{case -1: { cout << "\n **** modification du pas de temps impossible,"
<< " le nouveau pas de temps calcule est plus petit que le pas de temps mini limite: "
<< pa.Deltatmini();
Sortie(1); break;
}
case 0: { break;} // cas normal
case 1: { pa.Modif_Detat_dans_borne(delta_t); break;} // cas ou on dépasse la borne maxi
};
};
if (modif_deltat && ((ParaGlob::NiveauImpression() > 4) || (permet_affichage > 2)))
cout << "\n --->>>> modif increment de temps de " << ancien_pas << " a " << delta_t;
////----debug
//if (delta_t > 0.5)
//cout <<"\n debug AlgoriTchamwa::Gestion_pas_de_temps "
// << endl;
////---fin debug
};
// sauvegarde du nouveau pas de temps critique
delta_t_critique = delta_t_essai;
break;
} //-- fin du cas == 2
case 3: // comme le cas deux, mais en plus on regarde par rapport à un temps proposé qui peut être limitatif
// on regarde si le pas de temps n'est pas devenu supérieur au pas critique
// ou au contraire que le pas de temps critique a augmenté
{ double l_sur_c = lesMail->Longueur_arrete_mini_sur_c(TEMPS_t);
double delta_t_essai = l_sur_c / (*phi_1) ;double ancien_pas=delta_t;
if (permet_affichage > 2)
cout << "\n pas de temps critique DFC: "<<l_sur_c
<< " pas de temps critique tchamwa: " << delta_t_essai << flush;
//on intervient: si le nouveau pas de temps critique est inférieur au pas de temps en cours
// ou si le nouveau pas de temps critique de l'algorithme est supérieur à l'ancien
if (( delta_t_essai < delta_t) || ( delta_t_essai > delta_t_critique))
{// mise à jour éventuel du pas de temps et du pas de temps maxi s'ils sont définit à partir du temps critique
bool modif = pa.Modif_Deltat_DeltatMaxi(delta_t_essai,l_sur_c); // mais pas de test vis-a-vis des bornes
// s'il y a eu modif du pas de temps on met à jour le pas courant
if (modif) {
delta_t = pa.Deltat(); modif_deltat=true;}
// maintenant on vérifie néanmoins le pas de temps choisit au cas où
// en particulier cette vérification n'est en fait pas utile si l'utilisateur
// à utilisé un facteur du pas critique < 1
if (pa.Limit_temps_stable() && ( delta_t > delta_t_essai))
{ if (pa.Coefficient_pas_critique_deltat() <= 1.)
{ delta_t = pa.Coefficient_pas_critique_deltat() * delta_t_essai;}
else { delta_t = delta_t_essai;};
// modification de l'increment de temps dans les paramètres
switch (pa.Modif_Deltat(delta_t))
{case -1: { if ((ParaGlob::NiveauImpression() > 0)|| (permet_affichage > 0))
cout << "\n **** modification du pas de temps impossible,"
<< " le nouveau pas de temps calcule est plus petit que le pas de temps mini limite: "
<< pa.Deltatmini();
Sortie(1); break;
}
case 0: { break;} // cas normal
case 1: { pa.Modif_Detat_dans_borne(delta_t); break;} // cas ou on dépasse la borne maxi
};
modif_deltat=true;
}
else // sinon on regarde simplement si la modification de temps est possible
{ switch (pa.Modif_Deltat(delta_t))
{case -1: { cout << "\n **** modification du pas de temps impossible,"
<< " le nouveau pas de temps calcule est plus petit que le pas de temps mini limite: "
<< pa.Deltatmini();
Sortie(1); break;
}
case 0: { break;} // cas normal
case 1: { pa.Modif_Detat_dans_borne(delta_t); break;} // cas ou on dépasse la borne maxi
};
};
};
// maintenant on regarde le temps proposé qui peut être également limitatif
if ((nouveau_dt != 0.) && (nouveau_dt < delta_t)) // il faut de nouveau réduire
{// on regarde si la modification de temps est possible
switch (pa.Modif_Deltat(nouveau_dt))
{case -1: { cout << "\n **** modification du pas de temps impossible,"
<< " le nouveau pas de temps calcule est plus petit que le pas de temps mini limite: "
<< pa.Deltatmini();
Sortie(1); break;
}
case 0: { modif_deltat = true; break;} // cas normal
case 1: { cout << "\n **** modification du pas de temps"<< pa.Deltatmini() << " anormale ,"
<< " on ne devrait pas passer ici !!! "
<< "\n AlgoriTchamwa::Gestion_pas_de_temps(...";
Sortie(1); break;} // cas ou on dépasse la borne maxi
};
delta_t = pa.Deltat();
};
if (modif_deltat && ((ParaGlob::NiveauImpression() > 4)||(permet_affichage > 2)))
cout << "\n --->>>> modif increment de temps de " << ancien_pas << " a " << delta_t;
////----debug
//if (delta_t > 0.5)
//cout <<"\n debug AlgoriTchamwa::Gestion_pas_de_temps "
// << endl;
////---fin debug
// sauvegarde du nouveau pas de temps critique
delta_t_critique = delta_t_essai;
break;
} //-- fin du cas == 3
default :
{cout << "\nErreur : valeur incorrecte du cas = " << cas << "\n";
cout << "AlgoriTchamwa::Gestion_pas_de_temps(... \n";
Sortie(1);
}
}; // fin du switch
// --<T W>-- mise à jour éventuelle des variables simplificatrices pour le temps
if ((cas==1) || ( modif_deltat))
{ deltat2 = delta_t * delta_t;
unsurdeltat = 1./delta_t; lambdaDeltat = lambdaa * delta_t;
alphaDelta_t = alphaa * delta_t; betaDeltat2 = betaa * (delta_t*delta_t);
gammaDelta_t = gammaa * delta_t;
};
// retour
return modif_deltat;
};
//------- 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 calcul
// certaines variables ont-été changés
// initialisation
void AlgoriTchamwa::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(DYNA_EXP_TCHAMWA); // 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);
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
charge->Change_temps_fin_non_stricte(1);
// --<T W>-- on se place dans le cadre de l'algorithme proposé par Tchamwa et Wielgoz
// 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 pb
{lesMail->Inactive_ddl(); lesMail->Active_un_type_ddl_particulier(X1);};
// 00 --<T W>-- 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);
// 01 --<T W>-- récup du pas de temps, proposé par l'utilisateur, initialisation et vérif / pas critique
this->Gestion_pas_de_temps(true,lesMail,1,0.); // 1 signifie qu'il y a initialisation
// 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 initials
// les conditions limites initiales de vitesse et d'accélération sont prise en compte
// de manière identiques à des ddl quelconques, ce ne sont pas des ddl fixé !!
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 trois tableau 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
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 --<DFC>-- 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
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
icharge = 0;
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
lesContacts->DefElemCont(0.); // au début le déplacement des noeuds est nul
};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
//--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
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,0.); // 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)
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);
// 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);
}
}
//--fin cas de restart et/ou de sauvegarde--------
// 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);
};
// 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*.
// s'il y a des conditions linéaires, elles sont prises en compte ici une fois pour toutes concernant la largeur de bande
mat_masse = Choix_matrice_masse(nbddl_X,mat_masse,lesMail,lesRef
,Ass1.Nb_cas_assemb(),lesContacts,lesCondLim);
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());
// on signale que l'on utilise un comportement matériel normal
lesLoisDeComp->Loi_simplifie(false);
// on calcul la matrice de masse qui est supposée identique dans le temps
// c'est-à-dire que l'on considère que la masse volumique est constante
Cal_matrice_masse(lesMail,Ass1,*mat_masse,diversStockage,lesRef,X1,lesFonctionsnD);
// on sauvegarde la matrice masse
(*mat_masse_sauve) = (*mat_masse);
// dans le cas où l'on utilise de l'amortissement numérique la matrice masse est modifiée
// en fait ici cela correspond à : ([M] + [C] delta t) avec [C] = nu [M]
// 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 (pa.Amort_visco_artificielle())
{ 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);
// dans le cas d'un test de convergence vers une solution statique pour laquelle il faut réserver de la place
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
if (amortissement_cinetique)
Algori::InitialiseAmortissementCinetique(); // initialisation des compteurs pour l'amortissement au cas ou
// 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 AlgoriTchamwa::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(DYNA_EXP_TCHAMWA); // 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
// mise à jour au cas où
Algori::MiseAJourAlgoMere(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
,diversStockage,charge,lesCondLim,lescontacts,resultats);
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 AlgoriTchamwa::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(DYNA_EXP_TCHAMWA); // transfert info
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// init var glob: du num d'itération. De manière arbitraire, en dynamique explicite
Transfert_ParaGlob_COMPTEUR_ITERATION_ALGO_GLOBAL(1); // on a toujours une seule itération
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// récup des entités
mise à jour 7.010
2023-05-03 17:23:49 +02:00
Assemblage& Ass1 = *Ass1_;//Assemblage& Ass2 = *Ass2_;
Assemblage& Ass3 = *Ass3_;
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +02:00
// 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
double max_delta_X=0.; // le maxi du delta X
double max_var_delta_X=0.; // idem d'une itération à l'autre
// boucle sur les increments de charge qui sont é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
bool arret=false; // booleen pour arrêter indépendamment de la charge
bool premier_calcul = true; // utilisé pour le contact
icharge++; // on incrémente le chargement -> donne le num d'incr du prochain incr chargé
while (((!charge->Fin(icharge))||(icharge == 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()))
{ 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
int inSol =0 ; // pointeur d'assemblage du maxi de variation de ddl
double maxDeltaDdl=0; // // maxi de variation de ddl
// initialisation de la variable puissance_précédente d'une itération à l'autre
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// double puis_precedente = 0.;
ajout de tous les algorithmes globaux
2021-09-26 14:31:23 +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);
// 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
// gestion du pas de temps, vérif / pas critique
this->Gestion_pas_de_temps(false,lesMail,2,0.); // 2 signifie cas courant
//--- ici on commence une procédure d'essai d'un temps: ce temps pourra être modifié en fonction du contact
// on considère uniquement 2 passes
bool seconde_passe_demandee = false; // tout d'abord la première passe
do
{ bool modif_temps = charge->Avance(); // avancement de la charge et donc du temps courant
//-- 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,0.); // 2 signifie cas courant
// calcul de l'increment
// initialisation des deux partie du second membre
vglobin.Zero();
vglobex.Zero();
if (pa.ContactType())
vcontact.Zero();
vglobaal.Zero(); // puissance totale
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
// 2_3 --<T W>-- 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,lesFonctionsnD,charge->MultCharge(),change_statut,cas_combi_ddl);
// 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
// dans le cas de l'existance de condition linéaire, la résolution détruit la matrice, il faut donc la récupéré à chaque temps
if (lesCondLim->ExisteCondiLimite())
(*mat_masse) = (*mat_masse_sauve);
// dans le cas ou il y a changement de statut il faut remettre à jour
// les conditions limites sur la matrice masse
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);
// récupération de la matrice masse sans conditions limites
(*mat_masse) = (*mat_masse_sauve);
// normalement sur la matrice visqueuse il n'y a rien n'a faire
// if (pa.Amort_visco_artificielle()) // initialisation de la matrice visqueuse
// { bool initial = false;
// Cal_mat_visqueux_num_expli(matrice_mas,mat_C_pt,delta_X,initial,vitesse_tdt);
// }
// 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();
};
// 1_0 --<T W>-- récupération (initialisation) des ddl position, vitesse et accélération
lesMail->Vect_loc_vers_glob(TEMPS_t,X1,X_t,X1);
lesMail->Vect_loc_vers_glob(TEMPS_t,V1,vitesse_t,V1);
lesMail->Vect_loc_vers_glob(TEMPS_t,GAMMA1,acceleration_t,GAMMA1);
// récupération au niveau global des ddl locaux à tdt avec conditions limite
// 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
// quelquesoit la valeur de phi du schéma de tchamwa, on utilise le schéma des différences finis
// centrés pour calculer les valeurs des ddl dans le cas de ddl bloqué
if (PremierDdlFamille(s.ty_prin) == X1)
// cas ou les Xi sont imposés, on calcul Vi et Gammai
{ X_Bl(ih) = X_tdt(ix);
V_Bl(ih) = (X_tdt(ix) - X_t(ix))*unsurdeltat ;
G_Bl(ih)= (V_Bl(ih)-vitesse_t(iv))*unsurdeltat ;
}
else if (PremierDdlFamille(s.ty_prin) == V1)
// cas ou les Vi sont imposés, calcul des Xi et Gammai
{ V_Bl(ih) = vitesse_tdt(iv);
G_Bl(ih) = (vitesse_tdt(iv)-vitesse_t(iv))*unsurdeltat ;
X_Bl(ih) = X_t(ix) + delta_t * vitesse_tdt(iv);
}
else if (PremierDdlFamille(s.ty_prin) == GAMMA1)
// cas ou les gammai sont imposés, calcul des Vi et Xi
{ G_Bl(ih) = acceleration_tdt(ig);
V_Bl(ih) = vitesse_t(iv) + delta_t * acceleration_t(iv);
X_Bl(ih) = X_t(ix) + delta_t * vitesse_t(iv)
+ deltat2 * acceleration_t(ig);
}
acceleration_t(ig) = G_Bl(ih); // pour le cas ou il y a relachement des conditions limites
// au prochain pas de temps
};
// 1_1 --<T W>-- calcul du champ de vitesse
vitesse_tdt = vitesse_t + lambdaDeltat * acceleration_t;
// 2_0 --<T W>-- calcul du champ de position à t+dt
X_tdt = X_t + alphaDelta_t * vitesse_t
+ betaDeltat2 * acceleration_t
+ gammaDelta_t * vitesse_tdt;
// -- 2-01 --<T W>-- -- traitement particulier à l'aide de pade du premier ordre 1-1
int nbddl_X = vitesse_t.Taille();
if ((type_cal_equilibre == 3) && (icharge>10))
{ if ((degre_num == 1) && (degre_deno == 1))
{for (int i=1;i<=nbddl_X;i++)
if (Dabs(vitesse_t(i)) > ConstMath::trespetit)
{double b=(vitesse_t(i)-acceleration_t(i)*X_t(i)/(2.*vitesse_t(i)));
double cc=-acceleration_t(i)/(2.*vitesse_t(i));
double denominateur=(1.+cc*delta_t);
double numerateur= (X_t(i)+b*delta_t);
X_tdt(i) = numerateur / denominateur;
// (X_t(i)+(vitesse_t(i)-acceleration_t(i)*X_t(i)/(2.*vitesse_t(i)))*delta_t)
// /(1.-acceleration_t(i)*delta_t/(2.*vitesse_t(i)));
vitesse_tdt(i)=(b*denominateur-cc*numerateur)/(denominateur*denominateur);
};
}
else if ((degre_num == 1) && (degre_deno == 2))
{for (int i=1;i<=nbddl_X;i++)
{ double r1=acceleration_t(i)*X_t(i)-2.*vitesse_t(i)*vitesse_t(i);
if ((Dabs(vitesse_t(i)) > ConstMath::trespetit)
&& (Dabs(r1) > ConstMath::trespetit))
{ double cc=(acceleration_t(i)*vitesse_t(i))/r1;
double b = cc*X_t(i)+ vitesse_t(i);
double d=-cc* acceleration_t(i)*0.5/vitesse_t(i);
double denominateur=(1.+cc*delta_t+d*deltat2);
double numerateur= (X_t(i)+b*delta_t);
X_tdt(i) = numerateur / denominateur;
vitesse_tdt(i)=(b*denominateur-(cc+2.*delta_t)*numerateur) /(denominateur*denominateur);
// vitesse_tdt(i)= vitesse_t(i)/(1.-acceleration_t(i)/vitesse_t(i) * delta_t);
};
};
};
};
// -- 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);
};
// on fait un premier calcul du delta_X, qui n'intègre pas le contact, sinon au travers de l'accélération
// et vitesse précédemment calculée
// delta_X.Zero(); 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);
// 2_1 --<T W>-- 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
// a priori la ligne qui suit ne sert à rien: cf l'algo en dfc, où ça a été supprimé
lesMail->Vect_glob_vers_local(TEMPS_t,GAMMA1,acceleration_t,GAMMA1);
// 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
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)
{//double diam_mini = lesMail->Min_dist2Noeud_des_elements(TEMPS_t);
//lesContacts->DefElemCont(2. * diam_mini); //0.1);
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
};
};
// calcul des reactions de contact éventuelles (contact et frottement) et les puissances associées
if (pa.ContactType()) // 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
bool aff_incr=pa.Vrai_commande_sortie(icharge,temps_derniere_sauvegarde); // pour simplifier
if (!SecondMembreEnergContact(lesContacts,Ass1,vcontact,aff_incr)) arret = true;
};
// maintenant on regard si éventuellement le contact demanderait un pas de temps localement plus petit
if (!seconde_passe_demandee) // cas de la première passe
{double noueau_temps_ideal = lesContacts->Pas_de_temps_ideal();
//if (noueau_temps_ideal != 0.)
// cout << "\n nnoueau_temps_ideal= " << noueau_temps_ideal << endl;
if (this->Gestion_pas_de_temps(false,lesMail,3,noueau_temps_ideal)) // demande éventuelle de modif du pas de temps
seconde_passe_demandee = false;
};
} while (seconde_passe_demandee);
if (arret) break; // pour récupérer un pb éventuel dans SecondMembreEnergContact(
// affichage de l'increment de charge
bool 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======================================================================";
};
lesLoisDeComp->MiseAJour_umat_nbincr(icharge); // init pour les lois Umat éventuelles
// 3 --<T W>-- calcul des puissances internes et externe, (hors contact)
// 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;};
if (pa.ContactType())
vglobex += vcontact;
maxPuissExt = vglobex.Max_val_abs();
F_ext_tdt = vglobex; // sauvegarde des forces généralisées extérieures
// 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 = vglobin; // sauvegarde des forces généralisées intérieures
// second membre total
vglobaal += vglobex ;vglobaal += vglobin ;
// 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())
lesContacts->CalculReaction(F_ext_tdt,decol,Ass1.Nb_cas_assemb(),aff_incr);
// - 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())
tabCondLine = &(lesContacts->ConditionLin(Ass1.Nb_cas_assemb()));
// dans le cas où l'on utilise de l'amortissement numérique le second membre est modifiée
if (pa.Amort_visco_artificielle())
{ 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); // init de C éventuelle
vglobaal -= mat_C_pt->Prod_mat_vec(vitesse_t,forces_vis_num);
};
// initialisation des sauvegardes sur second membre (uniquement pour les gammai)
// non en fait pour les gammai
lesCondLim->InitSauve(Ass3.Nb_cas_assemb());
// sauvegarde des reactions aux ddl bloque (uniquement pour les Xi)
// non en fait pour les gammai
// 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
bool modif_repere = lesCondLim->CoLinCHrepere_int(*mat_masse,(vglobaal),Ass3.Nb_cas_assemb(),vglob_stat);
if (pa.ContactType()==1) // 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
// lesCondLim->ReacApresCHrepere(vglobin,lesMail,lesRef,Ass3.Nb_cas_assemb(),cas_combi_ddl);
// mise en place des conditions limites sur les Xi
// lesCondLim->ImposeConLimtdt(lesMail,lesRef,vglobal,Ass1.Nb_cas_assemb(),cas_combi_ddl);
// sur les Vi
// lesCondLim->ImposeConLimtdt(lesMail,lesRef,vglobal,Ass2.Nb_cas_assemb(),cas_combi_ddl);
// sur les Gammai
lesCondLim->ImposeConLimtdt(lesMail,lesRef,vglobaal,Ass3.Nb_cas_assemb(),cas_combi_ddl,vglob_stat);
// (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_incr)
InfoIncrementReac(lesMail,inReaction,maxReaction,Ass3.Nb_cas_assemb());
// examen de la convergence si nécessaire, utilisant le résidu
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
if (ArretEquilibreStatique() && (icharge>1))// 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
arret = Convergence(aff_incr,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) { arretResidu=1;}
else { cout << "\n critere equilibre statique satisfait pour l'increment : "
<< icharge << " " <<endl;
break;
};
}; // cas normal,
};
};
// 4 --<T W>-- calcul des accélérations
residu_final = vglobaal; // sauvegarde pour le post-traitement
// resolution simple (fonction du type de matrice)
// ou non suivant modif_repere
tempsResolSystemLineaire.Mise_en_route_du_comptage(); // temps cpu
if (!modif_repere)
{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() > 9)
{ 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());
// effacement du marquage de ddl bloque du au conditions lineaire imposée par l'entrée
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);
};
// calcul du maxi de variation de ddl
maxDeltaDdl = acceleration_tdt.Max_val_abs(inSol);
// sortie d'info sur l'increment concernant les variations de ddl
if ( aff_incr &&(ParaGlob::NiveauImpression() > 1))
InfoIncrementDdl(lesMail,inSol,maxDeltaDdl,Ass3.Nb_cas_assemb());
// calcul des énergies et affichage des balances
// delta_X.Zero(); 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 (icharge==1)// 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 éventuelle de l'amortissement cinétique
int relax_vit_acce = AmortissementCinetique(delta_X,1.,X_tdt,*mat_masse_sauve,icharge,vitesse_tdt);
// s'il y a amortissement cinétique il faut re-updater les vitesses
if (Abs(relax_vit_acce) == 1) {lesMail->Vect_glob_vers_local(TEMPS_tdt,V1,vitesse_tdt,V1);}
// examen de la convergence éventuelle, utilisant le déplacement et/ou le résidu
if (Pilotage_fin_relaxation_et_ou_residu(relax_vit_acce,0,icharge,arretResidu,arret))
break;
// si on est sans validation, on ne fait rien, sinon on valide l'incrément
// avec sauvegarde éventuelle
if (validation_calcul)
{// passage des accélérations calculées aux niveaux des maillages
lesMail->Vect_glob_vers_local(TEMPS_tdt,GAMMA1,acceleration_tdt,GAMMA1);
// 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
Algori::TdtversT();
if (pa.ContactType())
{ // 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)
};
// // actualisation des éléments de contact et éventuellement suppression
// if (pa.ContactType()) lesContacts->Actualisation(); // des éléments qui ne sont plus en contact
// 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);};
// sauvegarde de l'incrément si nécessaire
tempsCalEquilibre.Arret_du_comptage(); // on arrête le compteur pour la sortie
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
icharge++;
Transfert_ParaGlob_COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL(icharge);
}; // -- fin du test: if (validation_calcul)
brestart = false; // dans le cas où l'on était en restart, on passe l'indicateur en cas courant
// test de fin de calcul effectue dans charge via : charge->Fin()
// --cas particulier où la sortie de la boucle est pilotée
// ici, cas d'un calcul explicite au sens classique, il n'y a pas de notion de convergence
// celle-ci est toujours considérée comme effective
if (sortie_boucle_controler_par_fctnD)
break;
};
// on remet à jour le nombre d'incréments qui ont été effectués:
if (validation_calcul)
{icharge--;
Transfert_ParaGlob_COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL(icharge);
}
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);} ;
// fin des calculs
tempsCalEquilibre.Arret_du_comptage(); // temps cpu
};
// dernière passe
void AlgoriTchamwa::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(DYNA_EXP_TCHAMWA); // transfert info
type_incre = OrdreVisu::DERNIER_INCRE;
VisuAuFilDuCalcul(paraGlob,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()));
// 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);
};
Reference in a new issue Copy permalink