// This file is part of the Herezh++ application. // // The finite element software Herezh++ is dedicated to the field // of mechanics for large transformations of solid structures. // It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600) // INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) . // // Herezh++ is distributed under GPL 3 license ou ultérieure. // // Copyright (C) 1997-2022 Université Bretagne Sud (France) // AUTHOR : Gérard Rio // E-MAIL : gerardrio56@free.fr // // This program is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, // or (at your option) any later version. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty // of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. // See the GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program. If not, see . // // For more information, please consult: . #include "AlgoriDynaExpli_zhai.h" // CONSTRUCTEURS : AlgoriDynaExpli_zhai::AlgoriDynaExpli_zhai () : // par defaut Algori() // liste de variables de travail déclarées ici pour éviter le passage de paramètre entre les // méthodes internes à la classe ,grand_phi(NULL),phi_minus(NULL),gamma_pt(NULL),beta(NULL) ,delta_t(0.),deltat2(0.),unsurdeltat(0.) ,undemiplusgrandphi_deltat2(0.),grandphi_deltat2(0.),unpluphi_delta_t(0.),phi_delta_t(0.) ,un_moins_gamma_delta_t(0.),gamma_delta_t(0.),undemi_moins_beta_deltat2(0.),beta_deltat2(0.) // -------------------------------------------------------------------------------------- // -- 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(),acceleration_prec() ,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 -- // ------------------------------------------------------------------------------------------ { }; // constructeur en fonction du type de calcul et du sous type // il y a ici lecture des parametres attaches au type AlgoriDynaExpli_zhai::AlgoriDynaExpli_zhai (const bool avec_typeDeCal ,const list & soustype ,const list & avec_soustypeDeCal ,UtilLecture& entreePrinc) : Algori(DYNA_EXP_ZHAI,avec_typeDeCal,soustype,avec_soustypeDeCal,entreePrinc) ,grand_phi(NULL),phi_minus(NULL),gamma_pt(NULL),beta(NULL) // liste de variables de travail déclarées ici pour éviter le passage de paramètre entre les // méthodes internes à la classe ,delta_t(0.),deltat2(0.),unsurdeltat(0.) ,undemiplusgrandphi_deltat2(0.),grandphi_deltat2(0.),unpluphi_delta_t(0.),phi_delta_t(0.) ,un_moins_gamma_delta_t(0.),gamma_delta_t(0.),undemi_moins_beta_deltat2(0.),beta_deltat2(0.) // -------------------------------------------------------------------------------------- // -- 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(),acceleration_prec() ,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 -- // ------------------------------------------------------------------------------------------ { // 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 AlgoriDynaExpli_zhai::AlgoriDynaExpli_zhai (const AlgoriDynaExpli_zhai& algo): Algori(algo) ,grand_phi(NULL),phi_minus(NULL),gamma_pt(NULL),beta(NULL) // liste de variables de travail déclarées ici pour éviter le passage de paramètre entre les // méthodes internes à la classe ,delta_t(0.),deltat2(0.),unsurdeltat(0.) ,undemiplusgrandphi_deltat2(0.),grandphi_deltat2(0.),unpluphi_delta_t(0.),phi_delta_t(0.) ,un_moins_gamma_delta_t(0.),gamma_delta_t(0.),undemi_moins_beta_deltat2(0.),beta_deltat2(0.) // -------------------------------------------------------------------------------------- // -- 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(),acceleration_prec() ,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 -- // ------------------------------------------------------------------------------------------ { // on relie les paramètre à phi_minus et à grandphi phi_minus = & paraTypeCalcul(1); grand_phi = & paraTypeCalcul(2); gamma_pt = & paraTypeCalcul(3); beta = & paraTypeCalcul(4); }; // destructeur AlgoriDynaExpli_zhai::~AlgoriDynaExpli_zhai () { 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_; }; // execution de l'algorithme dans le cas dynamique explicite, sans contact void AlgoriDynaExpli_zhai::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 ::const_iterator ili,ili_fin = soustypeDeCalcul->end(); list ::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 ::const_iterator ili,ili_fin = soustypeDeCalcul->end(); list ::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); }; //------- 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 AlgoriDynaExpli_zhai::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_ZHAI); // transfert info #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() ,distribution_CPU_algo.Tab_indique_CPU_en_cours() ,distribution_CPU_algo.Tableau_noeud_CPU_en_cours() ,distribution_CPU_algo.Tab_indique_noeud_CPU_en_cours()); }; #endif // 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 & type_pb = lesMail->Types_de_problemes(); bool purement_non_meca = true; list ::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); // ---- on se place dans le cadre de l'algorithme de différences finis centrées // 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);}; // 01 ---- récup du pas de temps, proposé par l'utilisateur, initialisation et vérif / pas critique this->Gestion_pas_de_temps(true,lesMail,1); // 1 signifie qu'il y a initialisation // ---- on se place dans le cadre de l'algorithme proposé par Chung-Lee // 00 ---- 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 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 ---- 6 vecteur pour une manipulation globale des positions vitesses et accélérations // vecteur qui globalise toutes les positions de l'ensemble des noeuds 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) ; acceleration_prec.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; // 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 // --- 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 }; //--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); // 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()) lesMail->Inactive_un_ddl_particulier(X3); // 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*. 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); 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 // on considère que le nombre de noeud n'a pas changé, ni le nombre de ddl void AlgoriDynaExpli_zhai::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_ZHAI); // 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 de l'équilibre // 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 AlgoriDynaExpli_zhai::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_ZHAI); // transfert info // 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 // 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 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 double max_delta_X=0.; // le maxi du delta X double max_var_delta_X=0.; // idem d'une itération à l'autre 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é // int icharge_precedent = icharge; // pour se souvenir du précédent icharge, ne sert que pour l'initialisation 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é && !arret) // on fait deux passages, une prédiction puis une correction for (int prediction = 1;prediction > -1; prediction--) { 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 // double puis_precedente = 0.; // 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 lesLoisDeComp->MiseAJour_umat_nbincr(icharge); // init pour les lois Umat éventuelles if(prediction) { // gestion du pas de temps, vérif / pas critique this->Gestion_pas_de_temps(false,lesMail,2); // 2 signifie cas courant 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); // 2 signifie cas courant }; // affichage de l'increment de charge bool aff_incr=pa.Vrai_commande_sortie(icharge,temps_derniere_sauvegarde); // pour simplifier if ((aff_incr)&& prediction) {cout << "\n======================================================================" << "\nINCREMENT DE CHARGE : " << icharge << " intensite " << charge->IntensiteCharge() << " t= " << charge->Temps_courant() << " dt= " << ParaGlob::Variables_de_temps().IncreTempsCourant() << "\n======================================================================"; }; // 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 ---- imposition des ddls bloqués // initialisation des coordonnees et des ddl a tdt en fonctions des // ddl imposes et de l'increment du chargement bool change_statut = false; // init des changements de statut if(prediction) { 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 // 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); // 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 ---- récupération (initialisation) des ddl position, vitesse et accélération if (prediction) // récup uniquement la première fois {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 ::iterator ie,iefin=li_gene_asso.end();; // def d'un iterator adoc int ih=1; // indice if (prediction) 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 des paramètre, 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 temsp }; if (prediction) {// 1_1 ---- calcul du champ de vitesse vitesse_tdt = vitesse_t + unpluphi_delta_t * acceleration_t - phi_delta_t * acceleration_prec; // 2_0 ---- calcul du champ de position à t+dt X_tdt = X_t + delta_t * vitesse_t + undemiplusgrandphi_deltat2 * acceleration_t - grandphi_deltat2 * acceleration_prec; } else // sinon c'est le résultat corrigé {// 1_1 ---- calcul du champ de vitesse vitesse_tdt = vitesse_t + un_moins_gamma_delta_t * acceleration_t + gamma_delta_t * acceleration_tdt; // 2_0 ---- calcul du champ de position à t+dt X_tdt = X_t + delta_t * vitesse_t + undemi_moins_beta_deltat2 * acceleration_t + beta_deltat2 * acceleration_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); }; // 2_1 ---- 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); // 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()) // 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 if (!SecondMembreEnergContact(lesContacts,Ass1,vcontact,aff_incr)) arret = true;break; }; // 3 ---- 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;}; // 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 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 ; // 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) list * listCondLine=NULL; if (pa.ContactType()) listCondLine = &(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 //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) // idem pour le contact conduisant à des conditions linéaires modif_repere = modif_repere || lesCondLim->CoLinCHrepere_ext(*mat_masse,(vglobaal),*listCondLine,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,*vglobaal,Ass1.Nb_cas_assemb(),cas_combi_ddl); // sur les Vi // lesCondLim->ImposeConLimtdt(lesMail,lesRef,*vglobaal,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 && (!prediction) // on ne regarde la convergence qu'à la phase finale ) { 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 << " " <-- 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() >= 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()); // 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); }; if (!prediction) { // 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) {// 4_1 ---- passage des valeurs calculées aux niveaux des maillages // passage des accélérations calculées aux niveaux des maillages // si prediction: -> prédiction d'accélération, qui sera utilisée au second passage lesMail->Vect_glob_vers_local(TEMPS_tdt,GAMMA1,acceleration_tdt,GAMMA1); if (!prediction) { acceleration_prec = acceleration_t; // sauvegarde de l'accélération précédente // 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 lesContacts->Actualisation(0); // si on n'a plus de projection // 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) }; // 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(); // arret du 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 // test de fin de calcul effectue dans charge via : charge->Fin() icharge++; Transfert_ParaGlob_COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL(icharge); brestart = false; // dans le cas où l'on était en restart, on passe l'indicateur en cas courant }; }; // --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);} ; tempsCalEquilibre.Arret_du_comptage(); // temps cpu }; // dernière passe void AlgoriDynaExpli_zhai::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_ZHAI); // 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); }; // Résolution du problème mécanique en explicite dynamique sans contact void AlgoriDynaExpli_zhai::Calcul_Equilibre(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_ZHAI); // transfert info // 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 // 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); // ---- on se place dans le cadre de l'algorithme proposé par zhai // dans le cas où l'on calcul des contraintes et/ou déformation et/ou un estimateur d'erreur // à chaque incrément, initialisation Tableau tenuXVG(3);tenuXVG(1)=X1;tenuXVG(2)=V1;tenuXVG(3)=GAMMA1; bool 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);}; // 000 ---- récup du pas de temps, proposé par l'utilisateur, initialisation et vérif / pas critique this->Gestion_pas_de_temps(true,lesMail,1); // 1 signifie qu'il y a initialisation // 00 ---- 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 int 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 Assemblage Ass1(lesMail->InitNouveauCasAssemblage(1)); 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 Assemblage Ass2(lesMail->InitNouveauCasAssemblage(1)); 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 Assemblage Ass3(lesMail->InitNouveauCasAssemblage(1)); 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 Tableau t_assemb(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 const int 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(); Vecteur X_Bl(ttsi),V_Bl(ttsi),G_Bl(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 // 6 vecteur pour une manipulation globale des positions vitesses et accélérations // vecteur qui globalise toutes les positions de l'ensemble des noeuds 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) ; // vecteur qui globalise toutes les accélérations Vecteur acceleration_prec(nbddl_X) ; // calcul des énergies 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 string sindic="calcul a t+dt: "; // indique que les calcul d'énergies sont à t et non a t+dt // initialisation du compteur d'increments de charge icharge = 1; // 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 // --- 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 }; //--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); // 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()) { // 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-------- // choix de la matrice de masse, qui est en fait celle qui correspond au ddl Xi Mat_abstraite* mat_masse=NULL;Mat_abstraite* mat_masse_sauve=NULL; // ici le numéro d'assemblage est celui de X car on projette bien sur des vitesses virtuelles c-a-d ddl X*. 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); Mat_abstraite& matrice_mas = *mat_masse;Mat_abstraite& matrice_mas_sauve = *mat_masse_sauve; // choix de la résolution if (matrice_mas.Type_matrice() == DIAGONALE) // dans le cas d'une matrice diagonale on force la résolution directe quelque soit l'entrée matrice_mas.Change_Choix_resolution(DIRECT_DIAGONAL,pa.Type_preconditionnement()); else matrice_mas.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,matrice_mas,diversStockage,lesRef,X1,lesFonctionsnD); // on sauvegarde la matrice masse matrice_mas_sauve = matrice_mas; // 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é Mat_abstraite * mat_C_pt=NULL; Vecteur forces_vis_num(0); // 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(matrice_mas,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());} }; }; Mat_abstraite & mat_C = *mat_C_pt; // 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,matrice_mas,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 matrice_mas.Preparation_resol(); OrdreVisu::EnumTypeIncre 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 tempsInitialisation.Arret_du_comptage(); // temps cpu tempsCalEquilibre.Mise_en_route_du_comptage(); // temps cpu // 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 double max_delta_X=0.; // le maxi du delta X double max_var_delta_X=0.; // idem d'une itération à l'autre 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é && !arret) { 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 // double puis_precedente = 0.; // 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); // 2 signifie cas courant 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); // 2 signifie cas courant // 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 // initialisation des coordonnees et des ddl a tdt en fonctions des // ddl imposes et de l'increment du chargement 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); lesCondLim->MiseAJour_condilineaire_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 // 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 matrice_mas = matrice_mas_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,matrice_mas ,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 matrice_mas.Preparation_resol(); } // 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 // 1_0 ---- 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 ::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 temsp }; // 1_1 ---- calcul du champ de vitesse vitesse_tdt = vitesse_t + unpluphi_delta_t * acceleration_t - phi_delta_t * acceleration_prec; // 2_0 ---- calcul du champ de position à t+dt X_tdt = X_t + delta_t * vitesse_t + undemiplusgrandphi_deltat2 * acceleration_t - grandphi_deltat2 * acceleration_prec; // -- 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); } // 2_1 ---- 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); // 3 ---- calcul des puissances internes et externe // 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;}; // 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 // calcul des maxi des puissances internes 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 ; // 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(matrice_mas_sauve,mat_C_pt,delta_X,false,vitesse_tdt); // init de C éventuelle vglobaal -= mat_C.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()); // 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); // 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); // 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 << " " <-- calcul des accélérations residu_final = vglobaal; // sauvegarde pour le post-traitement // resolution tempsResolSystemLineaire.Mise_en_route_du_comptage(); // temps cpu matrice_mas.Simple_Resol_systID_2 (vglobaal,acceleration_tdt,pa.Tolerance() ,pa.Nb_iter_nondirecte(),pa.Nb_vect_restart()); tempsResolSystemLineaire.Arret_du_comptage(); // temps cpu // calcul du maxi de variation de ddl //// maxDeltaDdl = sol->Max_val_abs(inSol); // sortie d'info sur l'increment concernant les variations de ddl ////InfoIncrementDdl(lesMail,inSol,maxDeltaDdl); // passage des accélérations calculées aux niveaux des maillages acceleration_prec = acceleration_t; // sauvegarde de l'accélération précédente 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(); // on valide l'activité des conditions limites et condition linéaires lesCondLim->Validation_blocage (lesRef,charge->Temps_courant()); // 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,matrice_mas_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,matrice_mas_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; // cas du calcul des énergies, passage des grandeurs de tdt à t Algori::TdtversT(); //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(); // arret du 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 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() icharge++; Transfert_ParaGlob_COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL(icharge); } // on remet à jour le nombre d'incréments qui ont été effectués: icharge--; Transfert_ParaGlob_COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL(icharge); // fin des calculs tempsCalEquilibre.Arret_du_comptage(); // temps cpu // passage finale dans le cas d'une visualisation au fil du calcul 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())); }; // lecture des paramètres du calcul void AlgoriDynaExpli_zhai::lecture_Parametres(UtilLecture& entreePrinc) { // dimensionnement paraTypeCalcul.Change_taille(4); Transfert_ParaGlob_ALGO_GLOBAL_ACTUEL(DYNA_EXP_ZHAI); // transfert info MotCle motCle; // ref aux mots cle deja_lue_entete_parametre = 1; // a priori pas de lecture d'entête bool lecture_effective = false; // 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 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; *(entreePrinc.entree) >> st1 ; if (st1 != "phi_minus=") { cout << "\n erreur en lecture du parametre de l'algorithme explicite " << "\n on attendait le mot : phi_minus= , au lieu de " << st1 << "\n AlgoriDynaExpli_zhai::lecture_Parametres( ... "; Sortie(1); } // lecture du parametre *(entreePrinc.entree) >> paraTypeCalcul(1); *(entreePrinc.entree) >> st1 ; if (st1 != "grand_phi=") { cout << "\n erreur en lecture du parametre de l'algorithme explicite " << "\n on attendait le mot : grand_phi= , au lieu de " << st1 << "\n AlgoriDynaExpli_zhai::lecture_Parametres( ... "; Sortie(1); } // lecture du parametre *(entreePrinc.entree) >> paraTypeCalcul(2); *(entreePrinc.entree) >> st1 ; if (st1 != "gamma=") { cout << "\n erreur en lecture du parametre de l'algorithme explicite " << "\n on attendait le mot : gamma= , au lieu de " << st1 << "\n AlgoriDynaExpli_zhai::lecture_Parametres( ... "; Sortie(1); } // lecture du parametre *(entreePrinc.entree) >> paraTypeCalcul(3); *(entreePrinc.entree) >> st1 ; if (st1 != "beta=") { cout << "\n erreur en lecture du parametre de l'algorithme explicite " << "\n on attendait le mot : beta= , au lieu de " << st1 << "\n AlgoriDynaExpli_zhai::lecture_Parametres( ... "; Sortie(1); } // lecture du parametre *(entreePrinc.entree) >> paraTypeCalcul(4); lecture_effective = true; } 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.) // et Newmark par défaut (cf. théorie zhai) { paraTypeCalcul(1) = 0.5;paraTypeCalcul(2) = 0.5; paraTypeCalcul(3) = 0.5;paraTypeCalcul(4) = 1./12.; } // on relie les paramètre à phi_minus et à grand_phi phi_minus = & paraTypeCalcul(1); grand_phi = & paraTypeCalcul(2); gamma_pt = & paraTypeCalcul(3); beta = & paraTypeCalcul(4); // 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 écrot uniquement les données variables (supposées comme telles) void AlgoriDynaExpli_zhai::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_minus= " << paraTypeCalcul(1) << " " << "grand_phi=" << paraTypeCalcul(2)<< " " << "gamma=" << paraTypeCalcul(3)<< " " << "beta=" << paraTypeCalcul(4); }; // (*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 AlgoriDynaExpli_zhai::Lecture_Base_info_Parametre(UtilLecture& entreePrinc,const int& cas,bool choix) {if (cas == 1) { // dimensionnement paraTypeCalcul.Change_taille(4); 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_minus=") { cout << "\n erreur en lecture du parametre de l'algorithme explicite " << "\n on attendait le mot : phi_minus= , au lieu de " << toto << "\n AlgoriDynaExpli_zhai::Lecture_Base_info_Parametre( ... "; Sortie(1); } *(ent) >> paraTypeCalcul(1) ; *(ent) >> toto ; if (toto != "grand_phi=") { cout << "\n erreur en lecture du parametre de l'algorithme explicite " << "\n on attendait le mot : grand_phi , au lieu de " << toto << "\n AlgoriDynaExpli_zhai::Lecture_Base_info_Parametre( ... "; Sortie(1); } *(ent) >> paraTypeCalcul(2) ; *(ent) >> toto ; if (toto != "gamma=") { cout << "\n erreur en lecture du parametre de l'algorithme explicite " << "\n on attendait le mot : gamma= , au lieu de " << toto << "\n AlgoriDynaExpli_zhai::lecture_Parametres( ... "; Sortie(1); }; // lecture du parametre *(ent) >> paraTypeCalcul(3); *(ent) >> toto ; if (toto != "beta=") { cout << "\n erreur en lecture du parametre de l'algorithme explicite " << "\n on attendait le mot : beta= , au lieu de " << toto << "\n AlgoriDynaExpli_zhai::lecture_Parametres( ... "; Sortie(1); }; // lecture du parametre *(ent) >> paraTypeCalcul(4); } else { // cas où la lecture est impossible, on attribut les valeurs par défaut paraTypeCalcul(1) = 0.5;paraTypeCalcul(2) = 0.5; paraTypeCalcul(3) = 0.5;paraTypeCalcul(4) = 1./12.; }; // on relie les paramètre à phi_minus et à grand_phi phi_minus = & paraTypeCalcul(1); grand_phi = & paraTypeCalcul(2); gamma_pt = & paraTypeCalcul(3); beta = & paraTypeCalcul(4); }; //-- fin du cas = 1 // rien si le cas est différent de 1 }; // création d'un fichier de commande: cas des paramètres spécifiques void AlgoriDynaExpli_zhai::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 explicite avec la modele de zhai |" << "\n#---------------------------------------------------------------------------------" << "\n" << "\n PARA_TYPE_DE_CALCUL" << "\n # ........................................................." << "\n # / facteur de pilotage phi_minus , grand_phi, gamma et beta /" << "\n #........................................................" << "\n phi_minus= 0.5 grand_phi= 0.5 gamma= 0.5 beta= 0.083333"; // 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 // 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 AlgoriDynaExpli_zhai::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 switch (cas) { case 1 : { // ---- récup du pas de temps, proposé par l'utilisateur delta_t = pa.Deltat(); double delta_t_old = delta_t; // ---- 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 ; if (((*phi_minus) > 1./2.) && ((*grand_phi) < (*phi_minus))) { delta_t_essai = l_sur_c/2.*sqrt((2.*(*phi_minus)-1.)/((*phi_minus)-(*grand_phi)) /(2.*(*phi_minus)+1.));} else if (((*phi_minus) > 1./2.) && ((*grand_phi) >= (*phi_minus))) { delta_t_essai = l_sur_c/2.*sqrt(2./(2.*(*grand_phi)-(*phi_minus)));} else if (((*phi_minus) == 1./2.) && ((*grand_phi) >= 1./2.)) { delta_t_essai = l_sur_c *sqrt(1./(4.*(*grand_phi)-1.));} else if (((*phi_minus) > 0.) && ((*phi_minus) < 1./2.) && ((*grand_phi) > (*phi_minus)) && ((*grand_phi) < (*phi_minus)+ 1./4.)) { delta_t_essai = MiN( l_sur_c /2.*sqrt(2./(2.*(*grand_phi)-(*phi_minus))), l_sur_c *sqrt(1./(4.*((*grand_phi)-(*phi_minus))+1.)));} else if (((*phi_minus) > 0.) && ((*phi_minus) < 1./2.) && ((*grand_phi) >= (*phi_minus)+1./4.)) { delta_t_essai = l_sur_c /2.*sqrt(2./(2.*(*grand_phi)-(*phi_minus)));} else { delta_t_essai = l_sur_c; if (ParaGlob::NiveauImpression() > 0) cout << "\n *** attention les valeurs de phi et grand_phi sont tels que le pas critique" << "\n de temps de peut pas être determine, on retiend celui des differences centrees"; }; // 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) 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; delta_t = delta_t_essai; // modification de l'increment de temps dans les paramètres switch (pa.Modif_Deltat(delta_t)) {case -1: { 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) cout << "\n >>>> rectification du pas de temps " < 0) cout << "\n >>>> rectification du pas de temps " <Longueur_arrete_mini_sur_c(TEMPS_t); double delta_t_essai = l_sur_c ;double ancien_pas=delta_t; // on intervient que si le nouveau pas de temps critique est inférieur au pas de temps en cours if ( delta_t_essai < delta_t) {// 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: { 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() > 0)) cout << "\n --->>>> modif increment de temps de " << ancien_pas << " a " << delta_t; } break; } //-- fin du cas == 2 default : {cout << "\nErreur : valeur incorrecte du cas = " << cas << "\n"; cout << "AlgoriDynaExpli_zhai::Gestion_pas_de_temps(... \n"; Sortie(1); } }; // fin du switch // ---- mise à jour éventuelle des variables simplificatrices pour le temps if ((cas==1) || ( modif_deltat)) {unsurdeltat = 1./delta_t; deltat2 = delta_t * delta_t; undemiplusgrandphi_deltat2 = (1./2.+(*grand_phi)) * deltat2; grandphi_deltat2 = (*grand_phi)* deltat2; unpluphi_delta_t = ((1.+(*phi_minus)) * delta_t); phi_delta_t = ((*phi_minus)* delta_t); //un_moins_gamma_delta_t,gamma_delta_t,undemi_moins_beta_deltat2,beta_deltat2 gamma_delta_t = (*gamma_pt)*delta_t; un_moins_gamma_delta_t = (1.-(*gamma_pt))*delta_t; undemi_moins_beta_deltat2 = (0.5-(*beta))*deltat2; beta_deltat2 = (*beta)*deltat2; }; // retour return modif_deltat; }; //---- 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 AlgoriDynaExpli_zhai::ActionInteractiveAlgo() { cout << "\n commande? "; return false; };