// 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-2021 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 "ImpliNonDynaCont.h" #include "MathUtil.h" /*#ifdef ENLINUX_STREAM #include "Profiler.h" #endif #define USE_PROFILE*/ // CONSTRUCTEURS : ImpliNonDynaCont::ImpliNonDynaCont () : // par defaut Algori() ,vglobin(),vglobex(),vglobaal(),vcontact(),forces_vis_num() { }; // constructeur en fonction du type de calcul // il y a ici lecture des parametres attaches au type ImpliNonDynaCont::ImpliNonDynaCont (const bool avec_typeDeCal ,const list & soustype ,const list & avec_soustypeDeCal ,UtilLecture& entreePrinc) : Algori(NON_DYNA_CONT,avec_typeDeCal,soustype,avec_soustypeDeCal,entreePrinc) ,vglobin(),vglobex(),vglobaal(),vcontact(),forces_vis_num() { // lecture des paramètres attachés au type de calcul (ici aucun) switch (entreePrinc.Lec_ent_info()) { case 0 : {// pour signaler à Algori qu'il n'y a pas eu de lecture de paramètre deja_lue_entete_parametre=0; // puis appel de la méthode de la classe mère Algori::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 ImpliNonDynaCont::ImpliNonDynaCont (const ImpliNonDynaCont& algo): Algori(algo) ,vglobin(),vglobex(),vglobaal(),vcontact(),forces_vis_num() {}; // destructeur ImpliNonDynaCont::~ImpliNonDynaCont () { }; // execution de l'algorithme dans le cas non dynamique, implicit, sans contact void ImpliNonDynaCont::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) { // 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 { // on ne continue que si on n'a pas dépasser le nombre d'incréments maxi ou le temps maxi if (! (charge->Fin(icharge,true) ) ) Calcul_Equilibre(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 { // 1- on fait le calcul d'équilibre // on ne continue que si on n'a pas dépasser le nombre d'incréments maxi ou le temps maxi if (! (charge->Fin(icharge,true) ) ) Calcul_Equilibre(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); }; }; } 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); }; }; } } }; // execution de l'algorithme dans le cas non dynamique avec prise en // compte du contact void ImpliNonDynaCont::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(NON_DYNA_CONT); // transfert info // on défini globalement que l'on a pas de combinaison des ddl c'est à dire que le seul ddl est X int cas_combi_ddl=0; // 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 le temps fin stricte, vrai par défaut en implicite charge->Change_temps_fin_non_stricte(0); // dans le cas où l'on calcul des contraintes et/ou déformation et/ou un estimateur d'erreur // à chaque incrément, initialisation 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);}; // on définit un nouveau cas d'assemblage // à travers la définition d'une instance de la classe assemblage Nb_assemb cas_a = lesMail->InitNouveauCasAssemblage(1); // récup du numéro Assemblage Ass(cas_a); // stockage numéro // 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()); // mise a zero de tous les ddl actifs et creation des tableaux a t+dt // les ddl de position ne sont pas mis a zero ! ils sont initialises // 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 actifs avec les conditions initiales lesCondLim->Initial(lesMail,lesRef,lesCourbes1D,lesFonctionsnD,true,cas_combi_ddl); // activation des ddl de position lesMail->Active_un_type_ddl_particulier(X1); lesMail->MiseAJourPointeurAssemblage(Ass.Nb_cas_assemb());// mise a jour des pointeurs d'assemblage int nbddl = lesMail->NbTotalDdlActifs(); // nb total de ddl // choix initial de la matrice de raideur du système linéaire Tableau < Mat_abstraite*> tab_mato(1);tab_mato(1)=NULL; tab_mato = Choix_matriciel(nbddl,tab_mato,lesMail,lesRef,Ass.Nb_cas_assemb(),lesCondLim,lescontacts); Mat_abstraite* mato=tab_mato(1); Mat_abstraite& matglob = *mato; // choix de la résolution matglob.Change_Choix_resolution(pa.Type_resolution(),pa.Type_preconditionnement()); // vérification d'une singularité éventuelle de la matrice de raideur VerifSingulariteRaideurMeca(nbddl,*lesMail); // def vecteurs globaux vglobin.Change_taille(nbddl); // puissance interne vglobex.Change_taille(nbddl); // puissance externe vglobaal.Change_taille(nbddl,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); // 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_tdt.Change_taille(nbddl); delta_X.Change_taille(nbddl); var_delta_X.Change_taille(nbddl); // vecteur qui globalise toutes les vitesses de l'ensemble des noeuds vitesse_t.Change_taille(nbddl);vitesse_tdt.Change_taille(nbddl); // vecteur qui globalise toutes les accélérations acceleration_t.Change_taille(nbddl);acceleration_tdt.Change_taille(nbddl) ; // calcul des énergies 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); F_ext_t.Change_taille(nbddl); // forces généralisées int et ext au pas précédent F_int_tdt.Change_taille(nbddl); F_ext_tdt.Change_taille(nbddl); // forces généralisées int et ext au pas actuel residu_final.Change_taille(nbddl); // pour la sauvegarde du résidu pour le post-traitement // définition d'un pointeur de fonction d'assemblage en fonction du type voulu // définition d'un pointeur de fonction d'assemblage (symétrique ou non) // utilisé dans la prise en compte des efforts extérieurs void (Assemblage::* assembMat) // le pointeur (Mat_abstraite & matglob,const Mat_abstraite & matloc, const DdlElement& tab_ddl,const Tableau&tab_noeud) = & Assemblage::AssembMatSym; if (pa.Symetrie_matrice()) // cas d'un assemblage symétrique assembMat = & Assemblage::AssembMatSym; else // cas d'un assemblage non symétrique assembMat = & Assemblage::AssembMatnonSym; // calcul de la largeur de bande relative aux maillages individuellement int demi,total; lesMail->Largeur_Bande(demi,total,Ass.Nb_cas_assemb()); total = std::min(total,nbddl); // boucle sur les increments de charge icharge = 1; //--** // un premier increment pour demmarer le processus //--** charge->Avance(); // premier increment de charge ????????? différent des autres algos // --- 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 int cal_front = lesMail->CreeElemFront(); 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; brestart = true; // 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()); // 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-------- // 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); }; OrdreVisu::EnumTypeIncre type_incre = OrdreVisu::PREMIER_INCRE; // pour la visualisation au fil du calcul // Tableau *> tesc; /* #ifdef USE_PROFILE // préparation du profiler if (!ProfilerInit(collectDetailed, bestTimeBase,4000, 50)) { #endif*/ /* #ifdef USE_PROFILE int erro; // sortie des infos pour le profiler erro= ProfilerDump("\pHerezh.prof"); ProfilerTerm(); } #endif*/ // Sortie(1); // def d'un indicateur pour le premier contact de chaque increment bool premiercontact = true; int niveau_substitution = 1; // par défaut on utilise la matrice de raideur mato = tab_mato(1) // mise à jour au cas où Algori::MiseAJourAlgoMere(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp,diversStockage ,charge,lesCondLim,lescontacts,resultats); tempsInitialisation.Arret_du_comptage(); // temps cpu tempsCalEquilibre.Mise_en_route_du_comptage(); // temps cpu // 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; // pour arrêt du calcul au niveau du pilotage 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é ) { // def de variables 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 double last_var_ddl_max=0; // maxi var ddl juste après résolutuion, que l'on sauvegarde d'une it à l'autre // 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,Ass.Nb_cas_assemb(),*lesRef); varExpor->RenseigneVarUtilisateur(*lesMail,*lesRef); lesMail->CalStatistique(); // calcul éventuel de statistiques Transfert_ParaGlob_NORME_CONVERGENCE(ConstMath::grand);// on met une norme grande par défaut au début if (premiercontact) // modification de la charge et du pas de temps qu'au premier passage // mais pas après un changement de statut { Pilotage_du_temps(charge,arret); // appel du Pilotage if (arret) break; }// pilotage -> arret du calcul // **** pour l'instant on ne met pas à jour les boites d'encombrement des différenents élément, de frontières // **** mais il faudra y penser surtout lorsque les déformations et déplacement deviendront important !!!!! bool aff_incr=pa.Vrai_commande_sortie(icharge,temps_derniere_sauvegarde); // pour simplifier if (premiercontact) { // affichage de l'increment de charge 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 et des conditions linéaires imposées bool change_statut = false; // init des changements de statut, ne sert pas ici lesCondLim->MiseAJour_tdt (pa.Multiplicateur(),lesMail,charge->Increment_de_Temps(),lesRef,charge->Temps_courant() ,lesCourbes1D,lesFonctionsnD,charge->MultCharge(),change_statut,cas_combi_ddl); // definition des elements de contact eventuels // - imposition des conditions de non penetration lescontacts->DefElemCont(lesMail->Max_var_dep_t_a_tdt()); } else cout << "\n=============================================================================" << "\n ....... re-analyse du contact ........ " << "\nINCREMENT DE CHARGE : " << icharge << " intensite " << charge->IntensiteCharge() << " t= " << charge->Temps_courant() << " dt= " << ParaGlob::Variables_de_temps().IncreTempsCourant() << "\n============================================================================"; lesLoisDeComp->MiseAJour_umat_nbincr(icharge); // init pour les lois Umat éventuelles // mise à jour éventuelle de la matrice de raideur en fonction du contact Mise_a_jour_Choix_matriciel_contact(tab_mato,Ass.Nb_cas_assemb(),lescontacts,niveau_substitution); mato=tab_mato(1); // boucle de convergence sur un increment Vecteur * sol; // pointeur du vecteur solution // boucle d'equilibre sur un pas // on ne réexamine pas le decollement ou de nouveaux contact a chaque iteration, // car il peut y avoir beaucoup d'oscillation, du au comportement matériel int compteur; // déclaré à l'extérieur de la boucle car utilisé après la boucle for (compteur = 0; compteur<= pa.Iterations();compteur++) //---//\\//\\// début de la boucle sur les itérations d'équilibres //\\//\\// {// initialisation de la matrice et du second membre matglob.Initialise (0.); vglobin.Zero(); vglobex.Zero(); // 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,Ass.Nb_cas_assemb(),*lesRef); varExpor->RenseigneVarUtilisateur(*lesMail,*lesRef); lesMail->CalStatistique(); // calcul éventuel de statistiques Transfert_ParaGlob_COMPTEUR_ITERATION_ALGO_GLOBAL(compteur); 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 // mise en place des conditions linéaires sur les ddl // les 3 lignes qui suivent ont été rajouté à la suite de la modification de la prise en compte // des conditions linéaires (version >= 6.425) maintenant il faudra tester, en particulier // peut-être faut-il prendre en compte le change_statut dans la matrice ??? // dans un premier temps se reporter à AlgoriNonDyna::CalEquilibre bool change_statut = false; lesCondLim->MiseAJour_condilineaire_tdt (pa.Multiplicateur(),lesMail,charge->Increment_de_Temps(),lesRef,charge->Temps_courant() ,lesCourbes1D,lesFonctionsnD,charge->MultCharge(),change_statut,cas_combi_ddl); // affichage ou non de l'itération bool aff_iteration = (pa.freq_affich_iter() > 0) ? (aff_incr && (compteur % pa.freq_affich_iter()==0) &&(compteur!=0)) : false ; // utilisation d'un comportement tangent simplifié si nécessaire if (compteur <= pa.Init_comp_tangent_simple() ) lesLoisDeComp->Loi_simplifie(true); else lesLoisDeComp->Loi_simplifie(false); lesLoisDeComp->MiseAJour_umat_nbiter(compteur); // init pour les lois Umat éventuelles // actualisation des conditions de contact qui peuvent changer // la largeur de bande, quand un noeud glisse d'une facette sur une voisine if (compteur != 0) { lescontacts->Actualisation(); // mise à jour éventuelle de la matrice de raideur en fonction du contact Mise_a_jour_Choix_matriciel_contact(tab_mato,Ass.Nb_cas_assemb(),lescontacts,niveau_substitution); mato=tab_mato(1); } // mise en place du chargement impose sur le second membre // et éventuellement sur la raideur en fonction de sur_raideur // mise en place du chargement impose, c-a-d calcul de la puissance externe // si pb on sort de la boucle if (!(charge->ChargeSMembreRaideur_Im_mecaSolid (Ass,lesMail,lesRef,vglobex,matglob,assembMat,pa,lesCourbes1D,lesFonctionsnD))) { Change_PhaseDeConvergence(-10);break;}; maxPuissExt = vglobex.Max_val_abs(); F_ext_tdt = vglobex; // sauvegarde des forces généralisées extérieures // - definition des conditions limites de contact const Tableau & tabCondLine = lescontacts->ConditionLin(Ass.Nb_cas_assemb()); // -- appel du calcul de la raideur et du second membre, énergies // dans le cas d'un calcul inexploitable arrêt de la boucle if (!RaidSmEner(lesMail,Ass,vglobin,matglob)) 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 ( vglobal est identique a vglobin) 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 bool decol; lescontacts->CalculReaction(vglobin,decol,Ass.Nb_cas_assemb(),aff_iteration); // initialisation des sauvegardes sur matrice et second membre lesCondLim->InitSauve(Ass.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,Ass.Nb_cas_assemb(),cas_combi_ddl); // mise en place de condition externe lineaires /* string toti; cout << "\n matrice et second membre avant mise en place du contact"; matglob.Affiche1(1,4,1,1,4,1); vglobal.Affiche();*/ // cout <<" \n une valeur continuer ? "; cin >> toti; // expression de la raideur et du second membre dans un nouveau repere lesCondLim->CoLinCHrepere_ext(matglob,vglobaal,tabCondLine,Ass.Nb_cas_assemb(),vglob_stat); /* cout << "\n matrice et second membre dans le nouveau repere"; matglob.Affiche1(1,4,1,1,4,1); vglobal.Affiche(); */ // cout <<" \n une valeur continuer ? "; cin >> toti; // sauvegarde des reactions aux ddl bloque // en fait ne sert à rien, car les réactions maxi sont calculées dans condlin, mais comme c'est peut-être un peu spécial ici on laisse // mais dans les autres algos c'est supprimé !!!!! lesCondLim->ReacApresCHrepere(vglobin,lesMail,lesRef,Ass.Nb_cas_assemb(),cas_combi_ddl); // blocage des conditions limites autres que le contact lesCondLim->ImposeConLimtdt(lesMail,lesRef,matglob,vglobaal ,Ass.Nb_cas_assemb(),cas_combi_ddl,vglob_stat); /* cout << "\n matrice et second membre apres blocage des CL autres que ceux du contact"; matglob.Affiche1(1,4,1,1,4,1); vglobal.Affiche();*/ // cout <<" \n une valeur continuer ? "; cin >> toti; // blocage des conditions lineaires de contact lesCondLim->CoLinBlocage(matglob,vglobaal,Ass.Nb_cas_assemb(),vglob_stat); /* cout << "\n matrice et second membre apres blocage des CL de contact"; matglob.Affiche1(1,4,1,1,4,1); vglobal.Affiche(); cout <<" \n une valeur continuer ? "; cin >> toti;*/ // calcul du maxi des reactions maxReaction = lesCondLim->MaxEffort(inReaction,Ass.Nb_cas_assemb()); // --- resolution /// matglob.Affiche(1,8,1,1,6,1); // -- pour le debug /// vglobal.Affiche(); // -- pour le debug // sortie d'info sur l'increment concernant les réactions if (compteur != 0) if (aff_iteration) InfoIncrementReac(lesMail,compteur,inReaction,maxReaction,Ass.Nb_cas_assemb()); // calcul des énergies et affichage des balances // récupération (initialisation) des ddl position à t et t+dt lesMail->Vect_loc_vers_glob(TEMPS_t,X1,X_t,X1); lesMail->Vect_loc_vers_glob(TEMPS_tdt,X1,X_tdt,X1); // calcul de la variation de ddl / delta t // delta_X.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); Vecteur forces_vis_num; // ne sert pas CalEnergieAffichage(delta_X,icharge,brestart,forces_vis_num); // test de convergence sur un increment // en tenant compte du contact (non decollement) if (Convergence(aff_iteration,last_var_ddl_max,vglobaal,maxPuissExt,maxPuissInt,maxReaction,compteur,arret) && !decol) { // sortie des itérations sauf si l'on est en loi simplifiée if (lesLoisDeComp->Test_loi_simplife() ) // cas d'une loi simplifié on remet normal lesLoisDeComp->Loi_simplifie(false); else // cas normal, break; } else if (arret) {break;} // cas ou la méthode Convergence() demande l'arret // sinon on continue // ici sol en fait = vecglob qui est ecrase par la resolution // ---- resolution ---- residu_final = vglobaal; // sauvegarde pour le post-traitement tempsResolSystemLineaire.Mise_en_route_du_comptage(); // temps cpu bool erreur_resolution_syst_lineaire = false; // init try {// ici sol en fait = vecglob qui est ecrase par la resolution sol = &(matglob.Resol_systID (vglobaal,pa.Tolerance(),pa.Nb_iter_nondirecte(),pa.Nb_vect_restart())); // affichage éventuelle du vecteur solution : deltat_ddl if (ParaGlob::NiveauImpression() >= 10) { string entete = " affichage du vecteur solution (delta_ddl) "; sol->Affichage_ecran(entete); }; // retour des ddl dans les reperes generaux, dans le cas où ils ont ete modifie // par des conditions linéaires lesCondLim->RepInitiaux( *sol,Ass.Nb_cas_assemb()); //---debug //cout << "\ndelta ddl "; for(int i=1;i<=3;i++) { cout << "\nnoe:"<Affiche(); // --- pour le debug } catch (ErrResolve_system_lineaire excep) // cas d'une d'une erreur survenue pendant la résolution { erreur_resolution_syst_lineaire = true; // on prépare l'arrêt } // car il faut néanmoins effacer les marques avant l'arrêt catch (ErrSortieFinale) // cas d'une direction voulue vers la sortie // on relance l'interuption pour le niveau supérieur { ErrSortieFinale toto; throw (toto); } catch ( ... ) { erreur_resolution_syst_lineaire = true; // on prépare l'arrêt } tempsResolSystemLineaire.Arret_du_comptage(); // temps cpu // effacement du marquage de ddl bloque du au conditions lineaire imposée par l'entrée lesCondLim->EffMarque(); if (pa.ContactType()) lescontacts->EffMarque(); if (erreur_resolution_syst_lineaire) {Change_PhaseDeConvergence(-9); // on signale une divergence due à la résolution break; // arrêt si on avait détecté une erreur à la résolution }; // affichage éventuelle du vecteur solution : deltat_ddl if (ParaGlob::NiveauImpression() >= 10) { string entete = " affichage du vecteur solution (delta_ddl) apres retour dans le repere principal"; sol->Affichage_ecran(entete); }; // calcul du maxi de variation de ddl maxDeltaDdl = sol->Max_val_abs(inSol); last_var_ddl_max=maxDeltaDdl; double maxDeltatDdl_signe=(*sol)(inSol); // récupération de la grandeur signée // pilotage à chaque itération: ramène: maxDeltaDdl,sol modifiés éventuellement et insol Pilotage_chaque_iteration(sol,maxDeltaDdl,compteur,lesMail); // sortie d'info sur l'increment concernant les variations de ddl if (aff_iteration&&(ParaGlob::NiveauImpression() > 1)) InfoIncrementDdl(lesMail,inSol,maxDeltatDdl_signe,Ass.Nb_cas_assemb()); // sol->Affiche(); // -- pour le debug // actualisation des ddl actifs a t+dt lesMail->PlusDelta_tdt(*sol,Ass.Nb_cas_assemb()); //---//\\//\\// fin de la boucle sur les itérations d'équilibres //\\//\\// } // effacement du marquage de ddl bloque du au conditions lineaire imposée par l'entrée // car ici on n'est pas passé par cette endroit si l'on a convergé if (compteur<= pa.Iterations()) {lescontacts->EffMarque(); lesCondLim->EffMarque();}; // gestion de la fin des itérations if (!Pilotage_fin_iteration_implicite(compteur)) { // cas d'une non convergence // comme on incrémente pas les ddl on remet cependant les ddl // et les grandeurs actives de tdt aux valeurs de ceux de t lesMail->TversTdt(); } else { // réexamen du contact pour voir s'il n'y a pas de nouveau element de contact if (lescontacts->Nouveau(lesMail->Max_var_dep_t_a_tdt())) // !!! il faudra prevoir un compteur sur le reexamen du contact premiercontact = false; // retour sur l'increment pour un re examen // avec les nouveaux elements de contact else { premiercontact = true; //pour le prochain increment // actualisation des ddl et des grandeurs actives de t+dt vers t lesMail->TdtversT(); lesMail->Mise_a_jour_boite_encombrement_elem_front(TEMPS_tdt); // cas du calcul des énergies, passage des grandeurs de tdt à t Algori::TdtversT(); // 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(X1);}; // sauvegarde de l'incrément si nécessaire tempsCalEquilibre.Arret_du_comptage(); // arrêt 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); } } // effacement du marquage de ddl bloque du au conditions lineaire de contact // déjà fait ??????????? sans doute à supprimer lescontacts->EffMarque(); } //-- fin du while ((charge->Fin())||(icharge == 1)) // exploitation des resultats en fonctions de l'instance de Resultats // on remet à jour le nombre d'incréments qui ont convergé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'état finale 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())); };