// FICHIER : Hysteresis_bulk.cc // CLASSE : Hysteresis_bulk // 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 "Debug.h" # include using namespace std; //introduces namespace std #include #include #include "Sortie.h" #include "ConstMath.h" #include "CharUtil.h" #include "Hysteresis_bulk.h" #include "ExceptionsLoiComp.h" #include "Enum_TypeQuelconque.h" #include "TypeQuelconqueParticulier.h" // ========== fonctions pour la classe de sauvegarde des résultats ========= // constructeur par défaut Hysteresis_bulk::SaveResulHysteresis_bulk::SaveResulHysteresis_bulk() : MPr_t(),MPr_tdt(),MPr_R() ,fonction_aide_t(0.),fonction_aide_tdt(0.),wprime_t(1.),wprime_tdt(1.) ,modif(0),MPr_R_t_a_tdt(),nb_coincidence(0) ,ip2(),indic_coin(),fct_aide_t_a_tdt(),fct_aide() ,map_type_quelconque() {// la première valeur de la fonction d'aide est l'infini fct_aide.push_front(ConstMath::tresgrand);MPr_R.push_front(0.); // // def du conteneur de grandeurs quelconques, initialisée à 0 // Grandeur_scalaire_double grand_courant(0.); // {TypeQuelconque typQ1(PRESSION_HYST_REF,SIG11,grand_courant); // map_type_quelconque[PRESSION_HYST_REF]=(typQ1); // }; // {TypeQuelconque typQ1(PRESSION_HYST_REF_M1,SIG11,grand_courant); // map_type_quelconque[PRESSION_HYST_REF_M1]=(typQ1); // }; // {TypeQuelconque typQ1(PRESSION_HYST_T,SIG11,grand_courant); // map_type_quelconque[PRESSION_HYST_T]=(typQ1); // }; }; // constructeur de copie Hysteresis_bulk::SaveResulHysteresis_bulk::SaveResulHysteresis_bulk(const SaveResulHysteresis_bulk& sav): MPr_t(sav.MPr_t),MPr_tdt(sav.MPr_tdt) ,MPr_R(sav.MPr_R) ,fonction_aide_t(sav.fonction_aide_t),fonction_aide_tdt(sav.fonction_aide_tdt) ,wprime_t(sav.wprime_t),wprime_tdt(sav.wprime_tdt) ,modif(sav.modif),MPr_R_t_a_tdt(sav.MPr_R_t_a_tdt) ,nb_coincidence(sav.nb_coincidence),fct_aide_t_a_tdt(sav.fct_aide_t_a_tdt) ,ip2(),indic_coin(sav.indic_coin),fct_aide(sav.fct_aide) ,map_type_quelconque(sav.map_type_quelconque) { if (fct_aide.size() >= 2) // pointage de ip2 {ip2= fct_aide.end(); ip2--;ip2--; } }; // a priori il n'y a rien à faire Hysteresis_bulk::SaveResulHysteresis_bulk::~SaveResulHysteresis_bulk() // destructeur {}; // affectation Loi_comp_abstraite::SaveResul & Hysteresis_bulk::SaveResulHysteresis_bulk::operator = ( const Loi_comp_abstraite::SaveResul & a) { Hysteresis_bulk::SaveResulHysteresis_bulk& sav = *((Hysteresis_bulk::SaveResulHysteresis_bulk*) &a); // recopie de toutes les grandeurs MPr_t = sav.MPr_t; MPr_tdt = sav.MPr_tdt; MPr_R_t_a_tdt = sav.MPr_R_t_a_tdt; MPr_R = sav.MPr_R; fonction_aide_t = sav.fonction_aide_t; fonction_aide_tdt = sav.fonction_aide_tdt; fct_aide = sav.fct_aide; wprime_t = sav.wprime_t; wprime_tdt = sav.wprime_tdt; if (fct_aide.size() >= 2) // pointage de ip2 {ip2= fct_aide.end(); ip2--;ip2--; }; modif = sav.modif; nb_coincidence = sav.nb_coincidence; fct_aide_t_a_tdt = sav.fct_aide_t_a_tdt; indic_coin = sav.indic_coin; indicateurs_resolution = sav.indicateurs_resolution; indicateurs_resolution_t = sav.indicateurs_resolution_t; map_type_quelconque = sav.map_type_quelconque; return *this; }; //------- lecture écriture dans base info ------- // cas donne le niveau de la récupération // = 1 : on récupère tout // = 2 : on récupère uniquement les données variables (supposées comme telles) void Hysteresis_bulk::SaveResulHysteresis_bulk::Lecture_base_info (ifstream& ent,const int ) { // ici toutes les données sont toujours a priori variables string toto; ent >> toto >> MPr_t ;MPr_tdt = MPr_t; ent >> toto >> fonction_aide_t; fonction_aide_tdt=fonction_aide_t; ent >> toto >> wprime_t; wprime_tdt=wprime_t; ent >> toto; // ---- lecture fonction d'aide #ifdef MISE_AU_POINT if (toto != "list_fct_aide") { cout << "\n erreur dans la lecture de la list_io fonction d'aide, on ne trouve pas le mot cle: debut_List_IO= " << "\n Hysteresis_bulk::SaveResulHysteresis_bulk::Lecture_base_info (... "; Sortie(1); } #endif int taille;ent >> toto >> taille; // lecture de la taille double un_elem; // on vide la liste actuelle fct_aide.clear(); for (int i=1;i<=taille;i++) // puis on lit { ent >> un_elem; // lecture fct_aide.push_back(un_elem); // enregistrement } if (fct_aide.size() >= 2) // pointage de ip2 {ip2= fct_aide.end(); ip2--;ip2--; } ent >> toto; // ---- lecture des points d'inversion #ifdef MISE_AU_POINT if (toto != "list_pt_inver") { cout << "\n erreur dans la lecture de la list_io: ref contrainte, on ne trouve pas le mot cle: debut_List_IO= " << "\n Hysteresis_bulk::SaveResulHysteresis_bulk::Lecture_base_info (... "; Sortie(1); } #endif ent >> toto >> taille; // lecture de la taille double elem; // on vide la liste actuelle MPr_R.clear(); for (int j=1;j<=taille;j++) // puis on lit { ent >> elem; // lecture MPr_R.push_back(elem); // enregistrement } //re-initialisation des variables de travail Init_debut_calcul(); }; // cas donne le niveau de sauvegarde // = 1 : on sauvegarde tout // = 2 : on sauvegarde uniquement les données variables //(supposées comme telles) void Hysteresis_bulk::SaveResulHysteresis_bulk::Ecriture_base_info(ofstream& sort,const int ) { // ici toutes les données sont toujours a priori variables sort << " MPr_t " << MPr_t; sort << " fct_aide_t " << fonction_aide_t; sort << " wprime_t " << wprime_t; // les listes: fonction d'aide et point d'inversion sort << " list_fct_aide " << "taille= " << fct_aide.size() << " "; List_io ::const_iterator iter_courant,iter_fin = fct_aide.end(); for (iter_courant=fct_aide.begin();iter_courant!=iter_fin;iter_courant++) { sort << setprecision(ParaGlob::NbdigdoCA()) << (*iter_courant) ; sort << " "; } sort << " list_pt_inver " << "taille= " << MPr_R.size() << " "; List_io ::const_iterator jter_courant,jter_fin = MPr_R.end(); for (jter_courant=MPr_R.begin();jter_courant!=jter_fin;jter_courant++) { sort << setprecision(ParaGlob::NbdigdoCA()) << (*jter_courant) ; sort << " "; } }; // mise à jour des informations transitoires void Hysteresis_bulk::SaveResulHysteresis_bulk::TdtversT() { // partie générale MPr_t = MPr_tdt;fonction_aide_t=fonction_aide_tdt; wprime_t = wprime_tdt; // partie spécifique List_io ::iterator posi0=MPr_R.end(); posi0--;// position minimal qui doit rester List_io ::iterator iafact_posi0=fct_aide.end(); iafact_posi0--; switch (modif) { case 0: // rien ne change sur le pas en terme de coincidence et d'inversion { break; } case 1: // cas où l'on a une ou plusieurs coincidences { // a chaque coincidence il faut supprimer 2 pt de ref, sauf si l'on revient // sur la courbe de première charge, et que l'on a fait une symétrie trac -> comp par exemple // dans ce cas on ne supprime qu'un point: en clair il faut toujours qu'il reste un point List_io ::iterator itens;List_io ::iterator ifct; List_io ::iterator icoin = indic_coin.begin(); for (int i=1;i<=nb_coincidence;i++,icoin++) { itens = MPr_R.begin(); ifct = fct_aide.begin(); List_io ::iterator iafact_posi1 = ifct; iafact_posi1++; #ifdef MISE_AU_POINT List_io ::iterator posi1 = itens; posi1++; // la nouvelle valeur if (itens == posi0) { cout << "\n erreur: SaveResulHysteresis_bulk::TdtversT()"; cout << " on cherche a appliquer une coincidence or le pointeur itens est deja en fin de liste " << endl ; }; //le test qui suit n'est valable que pour la dernière coïncidence if ((posi1 == posi0)&&(wprime_tdt != 1.)&&(i==nb_coincidence)&&((*icoin)==1)) { cout << "\n erreur: SaveResulHysteresis_bulk::TdtversT()"; cout << " on est de retour sur la courbe principale et wprime_tdt n'est pas egale a 1 !! " << endl ; }; #endif if (itens != posi0) MPr_R.erase(itens); if (ifct != iafact_posi0) fct_aide.erase(ifct); // arrivée ici on a dépilé une fois. // comme c'est une boucle on ne sait pas si arrivée ici on est sur une boucle secondaire ou non // on regarde s'il faut dépiler 2 fois if ((*icoin)==2) {// on refait les tests et on dépile itens = MPr_R.begin(); ifct = fct_aide.begin(); #ifdef MISE_AU_POINT posi1 = itens; posi1++; // la nouvelle valeur if (itens == posi0) { cout << "\n erreur: SaveResulHysteresis_bulk::TdtversT()"; cout << " on cherche a appliquer une coincidence (la seconde) or le pointeur itens est deja en fin de liste " << endl ; }; //le test qui suit n'est valable que pour la dernière coïncidence if ((posi1 == posi0)&&(wprime_tdt != 1.)&&(i==nb_coincidence)) { cout << "\n erreur: SaveResulHysteresis_bulk::TdtversT()"; cout << " on est de retour sur la courbe principale et wprime_tdt n'est pas egale a 1 !! " << endl ; }; #endif if (itens != posi0) MPr_R.erase(itens); // normalement le test est inutile, sauf dans le cas sans mise au point if (ifct != iafact_posi0) fct_aide.erase(ifct); // on vérifie le wprime_tdt au dernier passage if ((i==nb_coincidence)&&(MPr_R.begin()==posi0)&&(wprime_tdt != 1.)) { cout << "\n erreur: SaveResulHysteresis_bulk::TdtversT()"; cout << " apres coincidence, on est de retour sur la courbe principale et wprime_tdt n'est pas egale a 1 !! " << endl ; }; }; // // on dépile que si c'est une coincidence secondaire ou, dans le cas d'une coincidence avec // // la première charge, on dépile que si il reste 2 pt d'inversion (cas tout traction ou // // tout compression // if ((wprime_tdt != 1.) || // ((wprime_tdt == 1.)&&(MPr_R.size()!=1))) // {itens = MPr_R.begin(); // if (itens != MPr_R.end()) MPr_R.erase(itens); // ifct = fct_aide.begin(); // if (ifct != fct_aide.end()) fct_aide.erase(ifct); // } } break; } case 2: // cas où l'on a une ou plusieurs inversions { int nb_ref = fct_aide.size(); // récup du nombre de pt de ref actuel bool init_ip2 = false; if (nb_ref == 1) init_ip2=true; // pour gérer ip2 int nb_ref_new = fct_aide_t_a_tdt.size(); // récup du nombre de nouveau pt // on ajoute les différents points d'inversions List_io ::iterator iitens = MPr_R_t_a_tdt.begin(); List_io ::iterator iifct = fct_aide_t_a_tdt.begin(); for (int ii=1;ii<=nb_ref_new;ii++,iitens++,iifct++) { MPr_R.push_front(*iitens); // ----- débug (a virer ) ------ // // on vérifie que le niveau de la fonction d'aide est acceptable // if (*iifct > *fct_aide.begin()) // {cout << "\n erreur 3 " << endl;Sortie(1);} // ----- fin débug (a virer ) ------ fct_aide.push_front(*iifct); if (init_ip2) // uniquement pour mémoriser le W_a de la première charge { // cas où l'on doit initialiser ip2 (car ici le nombre de pt = 2) ip2 = fct_aide.begin(); init_ip2 = false; }; } break; } case 3: // cas où l'on a une ou plusieurs coincidence(s) et inversion (s) { // on itère sur indic_coin pour savoir si l'ordre des inversions et coincidence List_io ::iterator iindic,iindic_fin = indic_coin.end(); // on ajoute les différents points d'inversions List_io ::iterator iitens = MPr_R_t_a_tdt.begin(); List_io ::iterator iifct = fct_aide_t_a_tdt.begin(); // a chaque coincidence il faut supprimer 2 pt de ref List_io ::iterator itens;List_io ::iterator ifct; for (iindic= indic_coin.begin();iindic != iindic_fin;iindic++) { if (*iindic ) // cas d'une coincidence {itens = MPr_R.begin(); ifct = fct_aide.begin(); List_io ::iterator iafact_posi1 = ifct; iafact_posi1++; #ifdef MISE_AU_POINT List_io ::iterator posi1 = itens; posi1++; // la nouvelle valeur if (itens == posi0) { cout << "\n erreur: SaveResulHysteresis_bulk::TdtversT()"; cout << " on cherche a appliquer une coincidence (cas3) or le pointeur itens est deja en fin de liste " << endl ; }; #endif if (itens != posi0) MPr_R.erase(itens); if (ifct != iafact_posi0) fct_aide.erase(ifct); // arrivée ici on a dépilé une fois. // comme c'est une boucle on ne sait pas si arrivée ici on est sur une boucle secondaire ou non // on regarde s'il faut dépiler 2 fois if ((*iindic)==2) {// on refait les tests et on dépile itens = MPr_R.begin(); ifct = fct_aide.begin(); #ifdef MISE_AU_POINT posi1 = itens; posi1++; // la nouvelle valeur if (itens == posi0) { cout << "\n erreur: SaveResulHysteresis_bulk::TdtversT()"; cout << " on cherche a appliquer une coincidence (la seconde) or le pointeur itens est deja en fin de liste " << endl ; }; #endif if (itens != posi0) MPr_R.erase(itens); // normalement le test est inutil, sauf dans le cas sans mise au point if (ifct != iafact_posi0) fct_aide.erase(ifct); // on vérifie le wprime_tdt au dernier passage }; } else {// cas d'une inversion MPr_R.push_front(*iitens); // ----- débug (a virer ) ------ // // on vérifie que le niveau de la fonction d'aide est acceptable // if (*iifct > *fct_aide.begin()) // {cout << "\n erreur 4 " << endl;Sortie(1);} // ----- fin débug (a virer ) ------ fct_aide.push_front(*iifct); iitens++;iifct++; // incrémentation pour la fois suivante }; }; // mise à jour d'ip2 if ( fct_aide_t_a_tdt.size() > 1) { ip2 = fct_aide.end(); ip2--;ip2--;}; break; } default: cout << "\n erreur dans la mise a jour des informations transitoire: modif= " << modif << "\n Hysteresis_bulk::SaveResulHysteresis_bulk::TversTdt() "; Sortie(1); }; // on met à jour éventuellement les indicateurs de résolution if (indicateurs_resolution.Taille()) { indicateurs_resolution_t = indicateurs_resolution; // on met les indicateurs à 0 indicateurs_resolution.Change_taille(5,0.); }; // mise à jour de la liste des grandeurs quelconques internes Mise_a_jour_map_type_quelconque(); //re-initialisation des variables de travail Init_debut_calcul(); }; void Hysteresis_bulk::SaveResulHysteresis_bulk::TversTdt() { // partie générale MPr_tdt = MPr_t;fonction_aide_tdt=fonction_aide_t; wprime_tdt=wprime_t; // mise à jour de la liste des grandeurs quelconques internes Mise_a_jour_map_type_quelconque(); //re-initialisation des variables de travail Init_debut_calcul(); }; // affichage à l'écran des infos void Hysteresis_bulk::SaveResulHysteresis_bulk::Affiche() const { cout << "\n SaveResulHysteresis_bulk: " ; cout << "\n MPr_t= " << MPr_t << " MPr_tdt= " << MPr_tdt << " fonction_aide_t= " << fonction_aide_t << " fonction_aide_tdt= " << fonction_aide_tdt << " wprime_t= " << wprime_t << " wprime_tdt= " << wprime_tdt ; cout << "\n list MPr_R_t_a_tdt: " << MPr_R_t_a_tdt; cout << "\n nb_coincidence= " << nb_coincidence; cout << "\n list fct_aide_t_a_tdt: " << fct_aide_t_a_tdt; cout << "\n list indic_coin: " << indic_coin; cout << "\n list MPr_R: " << MPr_R; cout << "\n list fct_aide: " << fct_aide << " "; }; // initialise les informations de travail concernant le pas de temps en cours void Hysteresis_bulk::SaveResulHysteresis_bulk::Init_debut_calcul() { //initialisation des variables de travail modif = 0; nb_coincidence=0; indic_coin.clear(); MPr_R_t_a_tdt.clear(); fct_aide_t_a_tdt.clear(); }; // mise à jour de la liste des grandeurs quelconques internes void Hysteresis_bulk::SaveResulHysteresis_bulk::Mise_a_jour_map_type_quelconque() { map < EnumTypeQuelconque , TypeQuelconque, std::less < EnumTypeQuelconque> >::iterator il ,ilfin=map_type_quelconque.end(); for (il=map_type_quelconque.begin();il != ilfin;il++) {EnumTypeQuelconque enu = (*il).first; switch (enu) {case PRESSION_HYST_REF: // tout d'abord PRESSION_HYST_REF = la dernière pression de référence { Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) map_type_quelconque[PRESSION_HYST_REF].Grandeur_pointee()); if (MPr_R.size() != 0) // donc on a une contrainte de ref { *(tyTQ.ConteneurDouble()) = - *(MPr_R.begin());} else // sinon on met à 0 { *(tyTQ.ConteneurDouble()) = 0.;} break; }; case PRESSION_HYST_REF_M1: // maintenant PRESSION_HYST_REF_M1 c-a-d la pression de référence à R-1 { Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) map_type_quelconque[PRESSION_HYST_REF_M1].Grandeur_pointee()); if (MPr_R.size() > 1) // donc on a deux contrainte de ref { List_io ::iterator it = MPr_R.begin(); it++; // ce qui correspond à R-1 *(tyTQ.ConteneurDouble()) = - (*it);} else // sinon on met à 0 { *(tyTQ.ConteneurDouble()) = 0.;} break; }; case PRESSION_HYST_T: // et enfin PRESSION_HYST_T c-a-d la pression à T { Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) map_type_quelconque[PRESSION_HYST_T].Grandeur_pointee()); *(tyTQ.ConteneurDouble()) = - MPr_t; break; }; case FCT_AIDE: // FCT_AIDE { Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) map_type_quelconque[FCT_AIDE].Grandeur_pointee()); *(tyTQ.ConteneurDouble()) = fonction_aide_tdt; break; }; case PRESSION_HYST:// PRESSION_HYST { Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) map_type_quelconque[PRESSION_HYST].Grandeur_pointee()); *(tyTQ.ConteneurDouble()) = - MPr_tdt; break; }; case NB_INVERSION: // NB_INVERSION { Grandeur_scalaire_entier& tyTQ= *((Grandeur_scalaire_entier*) map_type_quelconque[NB_INVERSION].Grandeur_pointee()); *(tyTQ.ConteneurEntier()) = MPr_R.size();; break; }; default:; // on ne fait rien sinon }; // puis tous les indicateurs qui sont éventuellement sauvegardé if (indicateurs_resolution.Taille()) {switch (enu) { case NB_INCRE_TOTAL_RESIDU: // NB_INCRE_TOTAL_RESIDU { Grandeur_scalaire_entier& tyTQ= *((Grandeur_scalaire_entier*) map_type_quelconque[NB_INCRE_TOTAL_RESIDU].Grandeur_pointee()); *(tyTQ.ConteneurEntier()) = indicateurs_resolution_t(1); break; }; case NB_ITER_TOTAL_RESIDU:// NB_ITER_TOTAL_RESIDU { Grandeur_scalaire_entier& tyTQ= *((Grandeur_scalaire_entier*) map_type_quelconque[NB_ITER_TOTAL_RESIDU].Grandeur_pointee()); *(tyTQ.ConteneurEntier()) = indicateurs_resolution_t(2); }; case NB_APPEL_FCT:// NB_APPEL_FCT { Grandeur_scalaire_entier& tyTQ= *((Grandeur_scalaire_entier*) map_type_quelconque[NB_APPEL_FCT].Grandeur_pointee()); *(tyTQ.ConteneurEntier()) = indicateurs_resolution_t(3); }; case NB_STEP: // NB_STEP { Grandeur_scalaire_entier& tyTQ= *((Grandeur_scalaire_entier*) map_type_quelconque[NB_STEP].Grandeur_pointee()); *(tyTQ.ConteneurEntier()) = indicateurs_resolution_t(4); }; case ERREUR_RK:// ERREUR_RK { Grandeur_scalaire_entier& tyTQ= *((Grandeur_scalaire_entier*) map_type_quelconque[ERREUR_RK].Grandeur_pointee()); *(tyTQ.ConteneurEntier()) = indicateurs_resolution_t(5); }; default:; }; }; } }; // ========== fin des fonctions pour la classe de sauvegarde des résultats ========= Hysteresis_bulk::Hysteresis_bulk () : // Constructeur par defaut Loi_comp_abstraite(HYSTERESIS_BULK,CAT_MECANIQUE,3),xnp(ConstMath::tresgrand) ,xmu(ConstMath::tresgrand),Qzero(ConstMath::tresgrand) ,xnp_temperature(NULL),Qzero_temperature(NULL),xmu_temperature(NULL) ,tolerance_residu(1.e-3),tolerance_residu_rel(1.e-5),nb_boucle_maxi(100),nb_sous_increment(4) ,maxi_delta_var_sig_sur_iter_pour_Newton(-1) ,tolerance_coincidence() ,MPr_t___tdt(),MPr_i___(),MPr_R() ,delta_MPr_Rat(),delta_MPr_Ratdt() ,delta_MPr_tatdt(),residuBH(),delta_V(0),delta__alpha_V(0.),wprime() ,residu(1),derResidu(1,1),alg_zero(),alg_edp() ,aucun_pt_inversion(true) ,rdelta_MPr_Ratdt() ,rdelta_MPr_tatdt() ,cas_kutta(5),erreurAbsolue(1.e-3),erreurRelative(1.e-5) ,MPr_point_(1),MPr_point(0.),MPr_tau(),MPr_tau_vect(1),delta_MPr_R_a_tau(0.) ,type_resolution_equa_constitutive(1),nbMaxiAppel(120) ,nb_maxInvCoinSurUnPas(4) ,depassement_Q0(1.1) // ---- affichage ,sortie_post(0) // -- variables de travail pour opérateur tangent ,I_x_I_HHHH(),I_xbarre_I_HHHH(),I_x_eps_HHHH(),Ixbarre_eps_HHHH() {// initialisation des paramètres de la résolution de newton alg_zero.Modif_prec_res_abs(tolerance_residu); alg_zero.Modif_prec_res_rel(tolerance_residu_rel); alg_zero.Modif_iter_max(nb_boucle_maxi); alg_zero.Modif_nbMaxiIncre(nb_sous_increment); tolerance_coincidence = tolerance_residu; // idem avec kutta alg_edp.Modif_nbMaxiAppel(nbMaxiAppel); // --- on remplit avec les grandeurs succeptible d'être utilisées // acces à listdeTouslesQuelc_dispo_localement list & list_tousQuelc = ListdeTouslesQuelc_dispo_localement(); list_tousQuelc.push_back(NB_INVERSION); list_tousQuelc.push_back(PRESSION_HYST_REF); list_tousQuelc.push_back(PRESSION_HYST); list_tousQuelc.push_back(NB_INCRE_TOTAL_RESIDU); list_tousQuelc.push_back(NB_ITER_TOTAL_RESIDU); list_tousQuelc.push_back(NB_APPEL_FCT); list_tousQuelc.push_back(NB_STEP); list_tousQuelc.push_back(ERREUR_RK); list_tousQuelc.push_back(FCT_AIDE); list_tousQuelc.push_back(PRESSION_HYST_REF_M1); list_tousQuelc.push_back(PRESSION_HYST_T); // on supprime les doublons localement list_tousQuelc.sort(); // on ordonne la liste list_tousQuelc.unique(); // suppression des doublons }; // Constructeur de copie Hysteresis_bulk::Hysteresis_bulk (const Hysteresis_bulk& loi) : Loi_comp_abstraite(loi),xnp(loi.xnp),xmu(loi.xmu),Qzero(loi.Qzero) ,xnp_temperature(loi.xnp_temperature),Qzero_temperature(loi.xnp_temperature) ,xmu_temperature(loi.xnp_temperature) ,tolerance_residu(loi.tolerance_residu),tolerance_residu_rel(loi.tolerance_residu_rel) ,nb_boucle_maxi(loi.nb_boucle_maxi) ,nb_sous_increment(loi.nb_sous_increment) ,maxi_delta_var_sig_sur_iter_pour_Newton(loi.maxi_delta_var_sig_sur_iter_pour_Newton) ,tolerance_coincidence(loi.tolerance_coincidence) ,MPr_t___tdt(loi.MPr_t___tdt),MPr_i___(loi.MPr_i___) ,MPr_R(loi.MPr_R) ,delta_MPr_Rat(loi.delta_MPr_Rat) ,delta_MPr_Ratdt(loi.delta_MPr_Ratdt) ,delta_MPr_tatdt(loi.delta_MPr_tatdt) ,residuBH(loi.residuBH),delta_V(loi.delta_V) ,delta__alpha_V(loi.delta__alpha_V),wprime(loi.wprime) ,residu(1),derResidu(1,1),alg_zero(),alg_edp() ,aucun_pt_inversion(loi.aucun_pt_inversion) ,rdelta_MPr_Ratdt(),rdelta_MPr_tatdt() ,cas_kutta(loi.cas_kutta),erreurAbsolue(loi.erreurAbsolue),erreurRelative(loi.erreurRelative) ,MPr_point_(1),MPr_point(0.),MPr_tau(),MPr_tau_vect(1),delta_MPr_R_a_tau(0.) ,type_resolution_equa_constitutive(loi.type_resolution_equa_constitutive) ,nbMaxiAppel(loi.nbMaxiAppel) ,nb_maxInvCoinSurUnPas(loi.nb_maxInvCoinSurUnPas) ,depassement_Q0(loi.depassement_Q0) // ---- affichage ,sortie_post(loi.sortie_post) // -- variables de travail pour opérateur tangent ,I_x_I_HHHH(),I_xbarre_I_HHHH(),I_x_eps_HHHH(),Ixbarre_eps_HHHH() {// initialisation des paramètres de la résolution de newton alg_zero.Modif_prec_res_abs(tolerance_residu); alg_zero.Modif_prec_res_rel(tolerance_residu_rel); alg_zero.Modif_iter_max(nb_boucle_maxi); alg_zero.Modif_nbMaxiIncre(nb_sous_increment); // idem avec kutta alg_edp.Modif_nbMaxiAppel(nbMaxiAppel); // pour les pointeurs de courbes, on regarde s'il s'agit d'une courbe locale ou d'une courbe globale if (xnp_temperature != NULL) if (xnp_temperature->NomCourbe() == "_") xnp_temperature = Courbe1D::New_Courbe1D(*(loi.xnp_temperature)); if (Qzero_temperature != NULL) if (Qzero_temperature->NomCourbe() == "_") Qzero_temperature = Courbe1D::New_Courbe1D(*(loi.Qzero_temperature));; if (xmu_temperature != NULL) if (xmu_temperature->NomCourbe() == "_") xmu_temperature = Courbe1D::New_Courbe1D(*(loi.xmu_temperature));; // --- on remplit avec les grandeurs succeptible d'être utilisées // acces à listdeTouslesQuelc_dispo_localement list & list_tousQuelc = ListdeTouslesQuelc_dispo_localement(); list_tousQuelc.push_back(NB_INVERSION); list_tousQuelc.push_back(PRESSION_HYST_REF); list_tousQuelc.push_back(PRESSION_HYST); list_tousQuelc.push_back(NB_INCRE_TOTAL_RESIDU); list_tousQuelc.push_back(NB_ITER_TOTAL_RESIDU); list_tousQuelc.push_back(NB_APPEL_FCT); list_tousQuelc.push_back(NB_STEP); list_tousQuelc.push_back(ERREUR_RK); list_tousQuelc.push_back(FCT_AIDE); list_tousQuelc.push_back(PRESSION_HYST_REF_M1); list_tousQuelc.push_back(PRESSION_HYST_T); // on supprime les doublons localement list_tousQuelc.sort(); // on ordonne la liste list_tousQuelc.unique(); // suppression des doublons }; Hysteresis_bulk::~Hysteresis_bulk () // Destructeur { if (xnp_temperature != NULL) if (xnp_temperature->NomCourbe() == "_") delete xnp_temperature; if (Qzero_temperature != NULL) if (Qzero_temperature->NomCourbe() == "_") delete Qzero_temperature; if (xmu_temperature != NULL) if (xmu_temperature->NomCourbe() == "_") delete xmu_temperature; }; Hysteresis_bulk::SaveResul * Hysteresis_bulk::New_et_Initialise() { // création des stockages internes éventuels en fonction de sortie_post SaveResulHysteresis_bulk* pt = new SaveResulHysteresis_bulk(); // insertion éventuelle de conteneurs de grandeurs quelconque Insertion_conteneur_dans_save_result(pt); // retour return pt; }; // Lecture des donnees de la classe sur fichier void Hysteresis_bulk::LectureDonneesParticulieres (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D ,LesFonctions_nD& lesFonctionsnD) { // ---- lecture du paramètre de prager string nom; *(entreePrinc->entree) >> nom; if(nom != "np=") { cout << "\n erreur en lecture du parametre de prager " << " on attendait la chaine : np= "; cout << "\n Hysteresis_bulk::LectureDonneesParticulieres " << "(UtilLecture * entreePrinc) " << endl ; throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; // on regarde si le paramètre de prager est thermo dépendant if(strstr(entreePrinc->tablcar,"np_thermo_dependant_")!=0) { thermo_dependant=true; *(entreePrinc->entree) >> nom; if (nom != "np_thermo_dependant_") { cout << "\n erreur en lecture de la thermodependance de np, on aurait du lire le mot cle np_thermo_dependant_" << " suivi du nom d'une courbe de charge ou de la courbe elle meme "; entreePrinc->MessageBuffer("**erreur2 Hysteresis_bulk::LectureDonneesParticulieres (...**"); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); } // lecture de la loi d'évolution du paramètre de prager en fonction de la température *(entreePrinc->entree) >> nom; // on regarde si la courbe existe, si oui on récupère la référence if (lesCourbes1D.Existe(nom)) { xnp_temperature = lesCourbes1D.Trouve(nom); } else { // sinon il faut la lire maintenant string non_courbe("_"); xnp_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str())); // lecture de la courbe xnp_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc); }; entreePrinc->NouvelleDonnee(); // prepa du flot de lecture } else // sinon on lit directement le paramètre de prager {*(entreePrinc->entree) >> xnp; if (xnp <= 0) {cout << "\n erreur en lecture du parametre de prager np= " << xnp << " il doit etre superieur a 0 ! "; cout << "\n Hysteresis_bulk::LectureDonneesParticulieres " << "(UtilLecture * entreePrinc) " << endl ; throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; }; // ---- lecture de mu : coefficient de lame ---- *(entreePrinc->entree) >> nom; if(nom != "mu=") { cout << "\n erreur en lecture du coefficient mu de lame " << " on attendait la chaine : mu= "; cout << "\n Hysteresis_bulk::LectureDonneesParticulieres " << "(UtilLecture * entreePrinc) " << endl ; throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; // on regarde si le coefficient de lame est thermo dépendant if(strstr(entreePrinc->tablcar,"mu_thermo_dependant_")!=0) { thermo_dependant=true; *(entreePrinc->entree) >> nom; if (nom != "mu_thermo_dependant_") { cout << "\n erreur en lecture de la thermodependance de mu, on aurait du lire le mot cle mu_thermo_dependant_" << " suivi du nom d'une courbe de charge ou de la courbe elle meme "; entreePrinc->MessageBuffer("**erreur3 Hysteresis_bulk::LectureDonneesParticulieres (...**"); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; // lecture de la loi d'évolution du coefficient mu en fonction de la température *(entreePrinc->entree) >> nom; // on regarde si la courbe existe, si oui on récupère la référence if (lesCourbes1D.Existe(nom)) { xmu_temperature = lesCourbes1D.Trouve(nom); } else { // sinon il faut la lire maintenant string non_courbe("_"); xmu_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str())); // lecture de la courbe xmu_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc); }; entreePrinc->NouvelleDonnee(); // prepa du flot de lecture } else // sinon on lit directement le paramètre mu {*(entreePrinc->entree) >> xmu;}; // ---- lecture de la limite de plasticité infini ----- *(entreePrinc->entree) >> nom; if(nom != "Qzero=") { cout << "\n erreur en lecture du coefficient mu de lame " << " on attendait la chaine : Qzero= "; cout << "\n Hysteresis_bulk::LectureDonneesParticulieres " << "(UtilLecture * entreePrinc) " << endl ; throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; // on regarde si Qzero est thermo dépendant if(strstr(entreePrinc->tablcar,"Qzero_thermo_dependant_")!=0) { thermo_dependant=true; *(entreePrinc->entree) >> nom; if (nom != "Qzero_thermo_dependant_") { cout << "\n erreur en lecture de la thermodependance de Qzero, on aurait du lire le mot cle Qzero_thermo_dependant_" << " suivi du nom d'une courbe de charge ou de la courbe elle meme "; entreePrinc->MessageBuffer("**erreur4 Hysteresis_bulk::LectureDonneesParticulieres (...**"); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; // lecture de la loi d'évolution du coefficient Qzero en fonction de la température *(entreePrinc->entree) >> nom; // on regarde si la courbe existe, si oui on récupère la référence if (lesCourbes1D.Existe(nom)) { Qzero_temperature = lesCourbes1D.Trouve(nom); } else { // sinon il faut la lire maintenant string non_courbe("_"); Qzero_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str())); // lecture de la courbe Qzero_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc); }; entreePrinc->NouvelleDonnee(); // prepa du flot de lecture } else // sinon on lit directement le paramètre Qzero {*(entreePrinc->entree) >> Qzero;}; // --- lecture éventuelle des paramètres de réglage ---- // de l'algo de résolution de l'équation d'avancement temporel // --- lecture éventuelle des paramètres de réglage ---- // de l'algo de résolution de l'équation d'avancement temporel if(strstr(entreePrinc->tablcar,"avec_parametres_de_reglage_")!=0) {//type_resolution_equa_constitutive = 0; // valeur voulant dire qu'il n'y a pas eu de lecture entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg // on lit tant que l'on ne rencontre pas la ligne contenant "fin_parametres_reglage_Hysteresis_" // ou un nouveau mot clé global auquel cas il y a pb !! MotCle motCle; // ref aux mots cle while (strstr(entreePrinc->tablcar,"fin_parametres_reglage_Hysteresis_")==0) { // si on a un mot clé global dans la ligne courante c-a-d dans tablcar --> erreur if ( motCle.SimotCle(entreePrinc->tablcar)) { cout << "\n erreur de lecture des parametre de reglage d'hysteresis: on n'a pas trouve le mot cle " << " fin_parametres_reglage_Hysteresis_ et par contre la ligne courante contient un mot cle global "; entreePrinc->MessageBuffer("** erreur5 des parametres de reglage de la loi de comportement d'hyteresis**"); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; // lecture d'un mot clé *(entreePrinc->entree) >> nom; if ((entreePrinc->entree)->rdstate() == 0) {} // lecture normale #ifdef ENLINUX else if ((entreePrinc->entree)->fail()) // on a atteind la fin de la ligne et on appelle un nouvel enregistrement { entreePrinc->NouvelleDonnee(); // lecture d'un nouvelle enregistrement *(entreePrinc->entree) >>nom; } #else else if ((entreePrinc->entree)->eof()) // la lecture est bonne mais on a atteind la fin de la ligne { if(nom != "fin_parametres_reglage_Hysteresis_") {entreePrinc->NouvelleDonnee(); *(entreePrinc->entree) >> nom;}; } #endif else // cas d'une erreur de lecture { cout << "\n erreur de lecture inconnue "; entreePrinc->MessageBuffer("** erreur5 des parametres de reglage de la loi de comportement d'hyteresis**"); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; // cas du type de résolution if(nom == "type_de_resolution_") {*(entreePrinc->entree) >> type_resolution_equa_constitutive; if ((type_resolution_equa_constitutive < 1) || (type_resolution_equa_constitutive> 3)) { cout << "\n erreur en lecture du type de resolution, on attendait un nombre compris entre 1 et 3 " << " et on a lu : " << type_resolution_equa_constitutive << " ce cas n'est pas implante (cf. doc) " << "\n Hysteresis_bulk::LectureDonneesParticulieres(UtilLecture * entreePrinc) " << endl ; throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; } // // type de calcul de la matrice tangente // else if (nom == "type_calcul_comportement_tangent_") // {*(entreePrinc->entree) >> type_calcul_comportement_tangent; // if ((type_calcul_comportement_tangent < 1) || (type_calcul_comportement_tangent> 2)) // { cout << "\n erreur en lecture du type de calcul du comportement tangent," // << " on attendait un nombre compris entre 1 et 2 " // << " et on a lu : " << type_calcul_comportement_tangent << " ce cas n'est pas implante (cf. doc) " // << "\n Hysteresis_bulk::LectureDonneesParticulieres(UtilLecture * entreePrinc) " << endl ; // throw (UtilLecture::ErrNouvelleDonnee(-1)); // Sortie(1); // }; // } // nombre d'itération maxi else if (nom == "nb_iteration_maxi_") {*(entreePrinc->entree) >> nb_boucle_maxi; alg_zero.Modif_iter_max(nb_boucle_maxi); } // nombre de dichotomie maxi else if (nom == "nb_dichotomie_maxi_") { *(entreePrinc->entree) >> nb_sous_increment; alg_zero.Modif_nbMaxiIncre(nb_sous_increment); } // tolérance absolue sur le résidu else if (nom == "tolerance_residu_") {*(entreePrinc->entree) >> tolerance_residu; alg_zero.Modif_prec_res_abs(tolerance_residu); } // tolérance relative sur le résidu else if (nom == "tolerance_residu_rel_") {*(entreePrinc->entree) >> tolerance_residu_rel; alg_zero.Modif_prec_res_rel(tolerance_residu_rel); } // le maxi de variation que l'on tolère d'une itération à l'autre else if (nom == "maxi_delta_var_sig_sur_iter_pour_Newton_") { *(entreePrinc->entree) >> maxi_delta_var_sig_sur_iter_pour_Newton; } // type d'algo de kutta else if (nom == "cas_kutta_") {*(entreePrinc->entree) >> cas_kutta; } // l'erreur absolue else if (nom == "erreurAbsolue_") {*(entreePrinc->entree) >> erreurAbsolue; } // l'erreur relative else if (nom == "erreurRelative_") {*(entreePrinc->entree) >> erreurRelative; } // le nombre d'appel maxi de la fonction f else if (nom == "nbMaxiAppel_") {*(entreePrinc->entree) >> nbMaxiAppel; alg_edp.Modif_nbMaxiAppel(nbMaxiAppel); } // tolérance sur la coïncidence else if (nom == "tolerance_coincidence_") {*(entreePrinc->entree) >> tolerance_coincidence; } // nombre maxi d'inversion sur le pas else if (nom == "nb_maxInvCoinSurUnPas_") {*(entreePrinc->entree) >> nb_maxInvCoinSurUnPas; } // else if (nom == "force_phi_zero_si_negatif_") // {*(entreePrinc->entree) >> force_phi_zero_si_negatif; // } else if (nom == "depassement_Q0_") {*(entreePrinc->entree) >> depassement_Q0; } // } // else if (nom == "mini_QepsilonBH_") // {*(entreePrinc->entree) >> mini_QepsilonBH; // if (mini_QepsilonBH < 0.) // { cout << "\n erreur en lecture sur mini_QepsilonBH, " // << "\n il doit etre positif on a lu : " << mini_QepsilonBH // << "\n Hysteresis_bulk::LectureDonneesParticulieres(UtilLecture * entreePrinc) " << endl ; // throw (UtilLecture::ErrNouvelleDonnee(-1)); // Sortie(1); // }; // } // // prise en compte des variations de w' et w'point dans le calcul de l'évolution tangente // else if (nom == "avecVarWprimeS_") // {*(entreePrinc->entree) >> avecVarWprimeS; // } // forcer un affichage particulier pour les méthodes else if (nom == "permet_affichage_") {Lecture_permet_affichage(entreePrinc,lesFonctionsnD); alg_edp.Change_niveau_affichage(Permet_affichage()); } // forcer un stockage des indicateurs de la résolution else if (nom == "sortie_post_") {*(entreePrinc->entree) >> sortie_post; } // sinon ce n'est pas un mot clé connu, on le signale else if(nom != "fin_parametres_reglage_Hysteresis_") { cout << "\n erreur en lecture d'un parametre, le mot cle est inconnu " << " on a lu : " << nom << endl; if (ParaGlob::NiveauImpression()>3) cout << "\n Hysteresis_bulk::LectureDonneesParticulieres(UtilLecture * entreePrinc) " << endl ; throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); } }; //-- fin du while }; //-- fin de la lecture des paramètres de réglage // appel au niveau de la classe mère Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire (*entreePrinc,lesFonctionsnD); // le niveau d'affichage alg_zero.Modif_affichage(Permet_affichage()); }; // affichage de la loi void Hysteresis_bulk::Affiche() const { cout << " \n loi_de_comportement Hysteresis_bulk " ; if ( xnp_temperature != NULL) { cout << " np thermo dependant " << " courbe np=f(T): " << xnp_temperature->NomCourbe() <<" ";} else { cout << " np= " << xnp ;}; if ( Qzero_temperature != NULL) { cout << " Qzero thermo dependant " << " courbe np=f(T): " << Qzero_temperature->NomCourbe() <<" ";} else { cout << " Qzero= " << Qzero ;}; if ( xmu_temperature != NULL) { cout << " xmu thermo dependant " << " courbe xmu=f(T): " << xmu_temperature->NomCourbe() <<" ";} else { cout << " mu= " << xmu ;}; if (type_resolution_equa_constitutive == 1) { cout << " \n type de resolution: linearisation, nbmax iteration= " << nb_boucle_maxi << " nbmax_dichotonmie= " << nb_sous_increment << " tolerance_residu_absolue= " << tolerance_residu << " tolerance_residu_relative= " << tolerance_residu_rel << " maxi_delta_var_sig_sur_iter_pour_Newton= " << maxi_delta_var_sig_sur_iter_pour_Newton; } else { cout << " \n type de resolution: Runge_Kutta " << cas_kutta << " erreurAbsolue " << erreurAbsolue << " erreurRelative= " << erreurRelative ; }; cout << " tolerance_coincidence= " << tolerance_coincidence << " nb_maxInvCoinSurUnPas= " << nb_maxInvCoinSurUnPas << " depassement_Q0= " << depassement_Q0 << " "; // niveau d'affichage Affiche_niveau_affichage(); cout << " sortie_post "<< sortie_post << " "; // appel de la classe mère Loi_comp_abstraite::Affiche_don_classe_abstraite(); }; // affichage et definition interactive des commandes particulières à chaques lois void Hysteresis_bulk::Info_commande_LoisDeComp(UtilLecture& entreePrinc) { ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier cout << "\n definition standart (rep o) ou exemples exhaustifs (rep n'importe quoi) ? "; string rep = "_"; // procédure de lecture avec prise en charge d'un retour chariot rep = lect_return_defaut(true,"o"); Qzero = 200.; xnp = 1.; xmu = 1000; // init à des exemples de valeurs sort << "\n#-------------------------------------------------------------------" << "\n# ....... loi_de_comportement d'hysteresis spherique ........ |" << "\n# para de prager(>=0) : mu : limite de plasticite |" << "\n#-------------------------------------------------------------------" << "\n np= " << setprecision(6) << xnp << " mu= " << setprecision(6) << xmu << " Qzero= "<< setprecision(6) << Qzero ; if ((rep != "o") && (rep != "O" ) && (rep != "0") ) { sort << "\n#-----------------------------------" << "\n# il est possible de definir des parametres thermo-dependants (1 ou 2 ou 3 parametres)" << "\n# par exemple pour les trois parametres on ecrit: " << "\n#-----------------------------------" << "\n# np= np_thermo_dependant_ courbe1 " << "\n# mu= xmu_temperature_ courbe4 " << "\n# Qzero= Qzero_temperature_ courbe3 " << "\n#-----------------------------------" << "\n# noter qu'apres la definition de chaque courbe, on change de ligne, d'une maniere inverse " << "\n# si la valeur du parametre est fixe, on poursuit sur la meme ligne. " << "\n# par exemple supposons np fixe, mu et Qzero thermo-dependant, on ecrit: " << "\n#-----------------------------------" << "\n# np= 2. mu= xmu_temperature_ courbe4 " << "\n# Qzero= Qzero_temperature_ courbe3 " << "\n#-----------------------------------" << "\n# il est egalement possible (mais pas obligatoire) de definir les parametres de reglage " << "\n# de la resolution. Dans ce cas, a la suite de la limite maxi de plasticite on indique " << "\n# le mot cle: avec_parametres_de_reglage_ " << "\n# on peut avoir un ou plusieur couple parametre-grandeur sur chaque ligne " << "\n# par contre la derniere ligne doit comporter uniquement le mot cle: " << "\n# fin_parametres_reglage_Hysteresis_ " << "\n# les differents parametres sont: " << "\n# " << "\n# le type de resolution: par linearisation (1) (valeur par defaut) ou integration directe via Runge_Kutta (2)" << "\n# le nombre d'iteration (pour (1)), le nombre de dichotomie(pour (1)), " << "\n# la tolerance abolue et relative sur le residu(pour (1))" << "\n# la norme maximale du delta sigma qui est permise a chaque iteration de Newton " << "\n# par defaut = -1 , si on veut une autre valeur: exe: " << "\n# maxi_delta_var_sig_sur_iter_pour_Newton_ 2 # donc 2 MPa a chaque iter " << "\n# si on donne une valeur negative (val par defaut), il n'y a pas de limite " << "\n# " << "\n# le type de Runge_kutta (pour (2)): 3 pour une methode imbriquee 2-3, 4 pour 3-4, 5 pour 4-5, " << "\n# l'erreur absolue sur la contrainte finale(pour (2)), l'erreur relative sur la contrainte finale(pour (2)) " << "\n# le nombre maxi d'appel de la fonction derivee (sigma_point) (pour (2)) " << "\n# la tolerance sur la detection des coincidences," << "\n# NB: tous les parametres n'ont pas a etre presents, il en faut seulement 1 au minimum " << "\n#-------------- limitation de l'intensite du deviateur ------------------" << "\n# " << "\n# - normalement l'intensite de QPression doit etre inferieure a la saturation Q_0 " << "\n# le parametre depassement_Q0 donne la valeur toleree de depassement, au-dessus de laquelle " << "\n# la contrainte n'est plus recevable " << "\n# ex: depassement_Q0_ 1.2 # signifie que l'on tolere 20% de deplacement" << "\n# " << "\n# " << "\n# -------------- affichage des erreurs et des warning ---------- " << "\n# - l'affichage normale est fonction du parametre global d'affichage gerer par le niveau d'affichage" << "\n# cependant pour des raisons par exemple de mise au point, il est possible de permettre l'affichage " << "\n# a un niveau particulier (mot cle : permet_affichage_ suivi d'un nombre entier) en plus de l'affichage normal. " << "\n# l'affichage s'effectuera donc en fonction de l'affichage normale et de l'affichage particulier." << "\n# Le fonctionnement de l'affichage particulier suit les mêmes règles que l'affichage globale" << "\n# soit permet_affichage_ est nulle (cas par defaut), dans ce cas l'affichage est fonction du niveau global" << "\n# soit permet_affichage_ vaut n par exemple, dans ce cas l'affichage est fonction uniquement de n " << "\n# " << "\n# ex: permet_affichage_ 5 " << "\n# --------------- acces aux indicateurs de la resolution par -------" << "\n# A chaque resolution, il est possible de stocker les indicateurs: nombre d'iteration, d'increment " << "\n# precision etc. les indicateurs sont renseignes en fonction du type de resolution " << "\n# le mot cle est sortie_post_ , par defaut il vaut 0, dans ce cas aucun indicateur n'est stoke" << "\n# s'il est different de 0, on peut acceder aux indicateurs en post-traitement" << "\n# seules les indicateurs en cours sont disponibles, il n'y a pas de stockage sur plusieurs increment " << "\n# " << "\n# ex: sortie_post_ 1 " << "\n# " << "\n# \n Exemple de declaration pour une resolution de l'equation linearisee: " << "\n#-------------------------------------------------------------------------" << "\n# ....... loi_de_comportement d'hysteresis spherique ........ |" << "\n# para de prager : mu : limite maxi de plasticite |" << "\n#-------------------------------------------------------------------------" << "\n# np= " << setprecision(6) << xnp << " mu= " << setprecision(6) << xmu << "\n# Qzero= "<< setprecision(6) << Qzero << " avec_parametres_de_reglage_ " << "\n# nb_iteration_maxi_ 20 nb_dichotomie_maxi_ 20 tolerance_residu_ 1.e-3 tolerance_residu_rel_ 1.e-5 " << "\n# tolerance_coincidence_ 1.e-3 " << "\n# \n Exemple de declaration pour une resolution avec une methode de Runge_Kutta: " << "\n#-------------------------------------------------------------------------" << "\n# ....... loi_de_comportement d'hysteresis spherique ........ |" << "\n# para de prager : mu : limite maxi de plasticite |" << "\n#-------------------------------------------------------------------------" << "\n# np= " << setprecision(6) << xnp << " mu= " << setprecision(6) << xmu << "\n# Qzero= "<< setprecision(6) << Qzero << " avec_parametres_de_reglage_ " << "\n# type_de_resolution_ 2 " << "\n# cas_kutta_ 5 " << "\n# erreurAbsolue_ 1.e-3 erreurRelative_ 1.e-5" << "\n# nbMaxiAppel_ 100" << "\n# tolerance_coincidence_ 1.e-3 nb_maxInvCoinSurUnPas_ 20 " << "\n# fin_parametres_reglage_Hysteresis_ "; }; sort << endl; // appel de la classe mère Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc); }; // test si la loi est complete int Hysteresis_bulk::TestComplet() { int ret = LoiAbstraiteGeneral::TestComplet(); if ((xnp_temperature == NULL) && (xnp == ConstMath::tresgrand)) { cout << " \n le parametre de prager n'est pas defini pour la loi " << Nom_comp(id_comp) << '\n'; ret = 0; } if ((xmu_temperature) && (xmu == ConstMath::tresgrand)) { cout << " \n la pente initiale mu n'est pas defini pour la loi " << Nom_comp(id_comp) << '\n'; ret = 0; } if ((Qzero_temperature) && (Qzero == ConstMath::tresgrand)) { cout << " \n la limite de plasticite n'est pas defini pour la loi " << Nom_comp(id_comp) << '\n'; ret = 0; } return ret; }; // calcul d'un module d'young équivalent à la loi, ceci pour un // chargement nul double Hysteresis_bulk::Module_young_equivalent(Enum_dure ,const Deformation & ,SaveResul * ) { // par défaut on prend un mu de 0.3 et on utilise l'équivalent d'une loi élastique // K = E/(2(1+nu)) -> E=2.6 * K if (xmu_temperature != NULL) xmu = xmu_temperature->Valeur(*temperature); return (2.6 * xmu); }; // récupération d'un module de compressibilité équivalent à la loi pour un chargement nul // il s'agit ici de la relation -pression = sigma_trace/3. = module de compressibilité * I_eps double Hysteresis_bulk::Module_compressibilite_equivalent(Enum_dure ,const Deformation & def) {if (xmu_temperature != NULL) xmu = xmu_temperature->Valeur(*temperature); return xmu; }; // récupération des grandeurs particulière (hors ddl ) // correspondant à liTQ // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière void Hysteresis_bulk::Grandeur_particuliere (bool ,List_io& liTQ,Loi_comp_abstraite::SaveResul * saveDon,list& decal) const { // ici on est en 3D et les grandeurs sont par principe en absolue, donc la variable absolue ne sert pas // on passe en revue la liste List_io::iterator itq,itqfin=liTQ.end(); list::iterator idecal=decal.begin(); for (itq=liTQ.begin();itq!=itqfin;itq++,idecal++) {TypeQuelconque& tipParticu = (*itq); // pour simplifier if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur switch (tipParticu.EnuTypeQuelconque().EnumTQ()) {case NB_INVERSION: // a) ----- cas du nombre d'inversion actuelle { SaveResulHysteresis_bulk & save_resul = *((SaveResulHysteresis_bulk*) saveDon); Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier tyTQ(1+(*idecal))=save_resul.MPr_R.size(); // cout << "\n save_resul.MPr_R.size() = "<& liTQ) const {// ici on est en 3D et les grandeurs sont par principe en absolue, donc la variable absolue ne sert pas Tableau tab_1(1); Tab_Grandeur_scalaire_double grand_courant(tab_1); // def d'un type quelconque représentatif à chaque grandeur // a priori ces grandeurs sont défini aux points d'intégration identique à la contrainte par exemple // enu_ddl_type_pt est définit dans la loi Abtraite générale // on n'ajoute que si sortie_post est vraie, sinon aucune grandeur n'est sauvegardé, donc on ne peut // plus y accèder //on regarde si ce type d'info existe déjà: si oui on augmente la taille du tableau, si non on crée // a) $$$ cas du nombre d'inversion actuelle {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; for (itq=liTQ.begin();itq!=itqfin;itq++) if ((*itq).EnuTypeQuelconque() == NB_INVERSION) { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier int taille = tyTQ.Taille()+1; tyTQ.Change_taille(taille); nexistePas = false; }; if (nexistePas) {TypeQuelconque typQ1(NB_INVERSION,enu_ddl_type_pt,grand_courant); liTQ.push_back(typQ1); }; }; // c) $$$ cas de la pression de référence actuelle {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; for (itq=liTQ.begin();itq!=itqfin;itq++) if ((*itq).EnuTypeQuelconque() == PRESSION_HYST_REF) {Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier int taille = tyTQ.Taille()+1; tyTQ.Change_taille(taille); nexistePas = false; }; if (nexistePas) {// def d'un type quelconque représentatif TypeQuelconque typQ(PRESSION_HYST_REF,SIG11,grand_courant); liTQ.push_back(typQ); }; }; // $$$ cas de la pression courante actuelle {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; for (itq=liTQ.begin();itq!=itqfin;itq++) if ((*itq).EnuTypeQuelconque() == PRESSION_HYST) {Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier int taille = tyTQ.Taille()+1; tyTQ.Change_taille(taille); nexistePas = false; }; if (nexistePas) {// def d'un type quelconque représentatif TypeQuelconque typQ(PRESSION_HYST,SIG11,grand_courant); liTQ.push_back(typQ); }; }; // ----- fonction d'aide {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; for (itq=liTQ.begin();itq!=itqfin;itq++) if ((*itq).EnuTypeQuelconque() == FCT_AIDE) {Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier int taille = tyTQ.Taille()+1; tyTQ.Change_taille(taille); nexistePas = false; }; if (nexistePas) {TypeQuelconque typQ1(FCT_AIDE,SIG11,grand_courant); liTQ.push_back(typQ1); }; }; // ---- la suite dépend de l'indicateur : sortie_post if (sortie_post) { // j) ----- NB_INCRE_TOTAL_RESIDU {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; for (itq=liTQ.begin();itq!=itqfin;itq++) if ((*itq).EnuTypeQuelconque() == NB_INCRE_TOTAL_RESIDU) {Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier int taille = tyTQ.Taille()+1; tyTQ.Change_taille(taille); nexistePas = false; }; if (nexistePas) {TypeQuelconque typQ1(NB_INCRE_TOTAL_RESIDU,SIG11,grand_courant); liTQ.push_back(typQ1); }; }; // k) ----- NB_ITER_TOTAL_RESIDU {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; for (itq=liTQ.begin();itq!=itqfin;itq++) if ((*itq).EnuTypeQuelconque() == NB_ITER_TOTAL_RESIDU) {Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier int taille = tyTQ.Taille()+1; tyTQ.Change_taille(taille); nexistePas = false; }; if (nexistePas) {TypeQuelconque typQ1(NB_ITER_TOTAL_RESIDU,SIG11,grand_courant); liTQ.push_back(typQ1); }; }; // l) ----- NB_APPEL_FCT {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; for (itq=liTQ.begin();itq!=itqfin;itq++) if ((*itq).EnuTypeQuelconque() == NB_APPEL_FCT) {Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier int taille = tyTQ.Taille()+1; tyTQ.Change_taille(taille); nexistePas = false; }; if (nexistePas) {TypeQuelconque typQ1(NB_APPEL_FCT,SIG11,grand_courant); liTQ.push_back(typQ1); }; }; // m) ----- NB_STEP {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; for (itq=liTQ.begin();itq!=itqfin;itq++) if ((*itq).EnuTypeQuelconque() == NB_STEP) {Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier int taille = tyTQ.Taille()+1; tyTQ.Change_taille(taille); nexistePas = false; }; if (nexistePas) {TypeQuelconque typQ1(NB_STEP,SIG11,grand_courant); liTQ.push_back(typQ1); }; }; // n) ----- ERREUR_RK {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; for (itq=liTQ.begin();itq!=itqfin;itq++) if ((*itq).EnuTypeQuelconque() == ERREUR_RK) {Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier int taille = tyTQ.Taille()+1; tyTQ.Change_taille(taille); nexistePas = false; }; if (nexistePas) {TypeQuelconque typQ1(ERREUR_RK,SIG11,grand_courant); liTQ.push_back(typQ1); }; }; }; // fin du cas ou sortie_post est actif, c-a-d que l'on veut des infos sur les indicateurs // de résolution }; // activation du stockage de grandeurs quelconques qui pourront ensuite être récupéré // via le conteneur SaveResul, si la grandeur n'existe pas ici, aucune action void Hysteresis_bulk::Activation_stockage_grandeurs_quelconques(list & listEnuQuelc) { // récup de la liste de stockage list & listlocale = ListQuelc_mis_en_acces_localement(); // on parcours la liste des grandeurs à activer // et on remplit la liste locale list ::iterator il, ilfin = listEnuQuelc.end(); for (il = listEnuQuelc.begin();il != ilfin; il++) // for (EnumTypeQuelconque enu : listEnuQuelc) // on ne remplit que s'il s'agit d'une grandeur qui peut-être accessible {switch (*il) {case NB_INVERSION : case PRESSION_HYST_REF: case PRESSION_HYST: case FCT_AIDE: case NB_INCRE_TOTAL_RESIDU : case NB_ITER_TOTAL_RESIDU: case NB_APPEL_FCT: case NB_STEP: case ERREUR_RK : case PRESSION_HYST_T: case PRESSION_HYST_REF_M1: listlocale.push_back(*(il)); break; default: ; // pour les autres cas on ne fait rien }; }; // on supprime les doublons localement listlocale.sort(); // on ordonne la liste listlocale.unique(); // suppression des doublons }; // insertion des conteneurs ad hoc concernant le stockage de grandeurs quelconques // passée en paramètre, dans le save result: ces conteneurs doivent être valides // c-a-d faire partie de listdeTouslesQuelc_dispo_localement void Hysteresis_bulk::Insertion_conteneur_dans_save_result(SaveResul * sr) { // récup de la liste de stockage list & listlocale = ListQuelc_mis_en_acces_localement(); // on spécialise saveresult SaveResulHysteresis_bulk & save_resul = *((SaveResulHysteresis_bulk*) sr); // -- autre stockage éventuel en fonction des grandeurs quelconques demandées par d'autres lois // on va regarder s'il faut activer ou non , pour le cas // de la récupération de grandeur quelconque au moment de l'utilisation de la loi List_io ::iterator jk,jkfin = listlocale.end(); for (jk=listlocale.begin();jk != jkfin;jk++) {EnumTypeQuelconque enu = *jk; switch (enu) {case PRESSION_HYST_REF: case PRESSION_HYST: case FCT_AIDE: case PRESSION_HYST_REF_M1: case PRESSION_HYST_T: { // on crée le conteneur ad hoc pour le passage d'info // def d'un conteneur de grandeurs quelconques, initialisée à 0 Grandeur_scalaire_double grand_courant(0.); TypeQuelconque typQ1(enu,EPS11,grand_courant); save_resul.map_type_quelconque[enu]=(typQ1); break; } case NB_INVERSION : { // on crée le conteneur ad hoc pour le passage d'info // def d'un conteneur de grandeurs quelconques, initialisée à 0 Grandeur_scalaire_entier grand_courant(0); TypeQuelconque typQ1(enu,EPS11,grand_courant); save_resul.map_type_quelconque[enu]=(typQ1); break; } case NB_INCRE_TOTAL_RESIDU : case NB_ITER_TOTAL_RESIDU: case NB_APPEL_FCT: case NB_STEP: case ERREUR_RK : // cas on ces grandeurs sont demandées {if (sortie_post) // si on dispose des infos {// on crée le conteneur ad hoc pour le passage d'info // def d'un conteneur de grandeurs quelconques, initialisée à 0 Grandeur_scalaire_entier grand_courant(0); TypeQuelconque typQ1(enu,EPS11,grand_courant); save_resul.map_type_quelconque[enu]=(typQ1); }; // sinon ce n'est pas dispo donc on ne fait rien break; } default: cout << "\n *** erreur on demande l'acces a : " << NomTypeQuelconque(enu) << " or celui-ci n'est pas dispo pour la loi "; this->Affiche(); cout << " revoir la mise en donnees ! " << endl; Sortie(1); }; }; }; // ========== codage des METHODES VIRTUELLES protegees:================ // calcul des contraintes void Hysteresis_bulk::Calcul_SigmaHH (TenseurHH & sigHH_t,TenseurBB& ,DdlElement & tab_ddl, TenseurBB & ,TenseurHH & ,BaseB& ,BaseH& ,TenseurBB& epsBB_, TenseurBB& delta_epsBB_,TenseurBB & gijBB_, TenseurHH & gijHH_,Tableau & d_gijBB_,double& ,double& , TenseurHH & sigHH_ ,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement ,const Met_abstraite::Expli_t_tdt& ex) { #ifdef MISE_AU_POINT if (epsBB_.Dimension() != 3) { cout << "\nErreur : la dimension devrait etre 3 !\n"; cout << " Hysteresis_bulk::Calcul_SigmaHH\n"; Sortie(1); }; if (tab_ddl.NbDdl() != d_gijBB_.Taille()) { cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_ !\n"; cout << " Hysteresis_bulk::Calcul_SigmaHH\n"; Sortie(1); }; #endif const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_); // passage en dim 3 const Tenseur3BB & delta_epsBB = *((Tenseur3BB*) &delta_epsBB_); // passage en dim 3 const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_); // " " " " const Tenseur3BB & gijBB = *((Tenseur3BB*) &gijBB_); // " " " " Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // " " " " // initialisation des variables de travail Tenseur3HH & sigHH_i = *((Tenseur3HH*) &sigHH_t); // " " " " SaveResulHysteresis_bulk & save_resul = *((SaveResulHysteresis_bulk*) saveResul); // initialisation des variables de travail save_resul.Init_debut_calcul(); // initialisation éventuelle des variables thermo-dépendantes Init_thermo_dependance(); // calcul de delta_V : incrément de variation relative de volume double V_tdt = (*ex.jacobien_tdt)/(*ex.jacobien_0); double V_t = (*ex.jacobien_t)/(*ex.jacobien_0); delta_V = V_tdt - V_t; // dans le cas de la prise en compte de la dilatation il faut changer V et delta_V if (dilatation) // on considère que la variation de volume thermique = la trace du tenseur de def thermique {double var_vol_thermique = (gijHH * (*epsBB_therm)).Trace(); V_tdt -= var_vol_thermique; // recup de l'incrément de temps double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant(); // on calcul un incrément en utilisant la vitesse de déformation thermique double delta_var_vol_thermique = (gijHH * (*DepsBB_therm)).Trace()*deltat; delta_V -= delta_var_vol_thermique; }; // === 1 == calcul de l'avancement temporel sur 1 pas, // ici il n'y a pas de différence entre orthonormeee ou pas au niveau du type de calcul int cas = 1; // donc cas avec dérivée de Jauman en mixte Avancement_temporel(gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ); LibereTenseur(); }; // calcul des contraintes a t+dt et de ses variations void Hysteresis_bulk::Calcul_DsigmaHH_tdt (TenseurHH & sigHH_t,TenseurBB& ,DdlElement & tab_ddl ,BaseB& ,TenseurBB & ,TenseurHH & ,BaseB& ,Tableau & ,BaseH& ,Tableau & ,TenseurBB & epsBB_tdt,Tableau & d_epsBB,TenseurBB & delta_epsBB_ ,TenseurBB & gijBB_tdt ,TenseurHH & gijHH_tdt ,Tableau & d_gijBB_tdt ,Tableau & d_gijHH_tdt,double& ,double& ,Vecteur& ,TenseurHH& sigHH_tdt,Tableau & d_sigHH ,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite ,double& module_cisaillement,const Met_abstraite::Impli& ex) { #ifdef MISE_AU_POINT if (epsBB_tdt.Dimension() != 3) { cout << "\nErreur : la dimension devrait etre 3 !\n"; cout << " Hysteresis_bulk::Calcul_DsigmaHH_tdt\n"; Sortie(1); }; if (tab_ddl.NbDdl() != d_gijBB_tdt.Taille()) { cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_tdt !\n"; cout << " Hysteresis_bulk::Calcul_DsigmaHH_tdt\n"; Sortie(1); }; #endif const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3 const Tenseur3BB & delta_epsBB = *((Tenseur3BB*) &delta_epsBB_); // passage en dim 3 const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); // " " " " const Tenseur3BB & gijBB = *((Tenseur3BB*) &gijBB_tdt); // " " " " Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // " " " " Tenseur3HH & sigHH_i = *((Tenseur3HH*) &sigHH_t); // " " " " SaveResulHysteresis_bulk & save_resul = *((SaveResulHysteresis_bulk*) saveResul); // initialisation des variables de travail save_resul.Init_debut_calcul(); // initialisation éventuelle des variables thermo-dépendantes Init_thermo_dependance(); // calcul de delta_V : incrément de variation relative de volume double un_sur_jaco_0 = 1./(*ex.jacobien_0); double V_tdt = (*ex.jacobien_tdt) * un_sur_jaco_0; double V_t = (*ex.jacobien_t) * un_sur_jaco_0; delta_V = V_tdt - V_t; // dans le cas de la prise en compte de la dilatation il faut changer V et delta_V if (dilatation) // on considère que la variation de volume thermique = la trace du tenseur de def thermique {double var_vol_thermique = (gijHH * (*epsBB_therm)).Trace(); V_tdt -= var_vol_thermique; // recup de l'incrément de temps double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant(); // on calcul un incrément en utilisant la vitesse de déformation thermique double delta_var_vol_thermique = (gijHH * (*DepsBB_therm)).Trace()*deltat; delta_V -= delta_var_vol_thermique; }; //// --- debug //cout << "\n debug Hysteresis_bulk::Calcul_DsigmaHH_tdt " // << "V_tdt= " << V_tdt << " V_t= " << V_t << " delta_V= "<< delta_V << endl; // //// --- fin debug // === 1 == calcul de l'avancement temporel sur 1 pas, int cas = 1; // cas normal avec dérivée de Jauman try // ce qui permet de donner les numéros d'éléments et de pti { Avancement_temporel(gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ); } catch (ErrNonConvergence_loiDeComportement excep) { // pour débugger ce qui se passe // initialisation des variables de travail save_resul.Init_debut_calcul(); // initialisation éventuelle des variables thermo-dépendantes Init_thermo_dependance(); //------- debug ------ #ifdef MISE_AU_POINT { cout << "\n debug Hysteresis_bulk::Calcul_DsigmaHH_tdt "; Avancement_temporel(gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ); }; #endif //------- fin debug ------ // on régénère une exception throw ErrNonConvergence_loiDeComportement(); } // ==== 2 === calcul de l'opérateur tangent ============= bool en_base_orthonormee = false;double T_d_Mpres_d_V; TenseurHHHH* d_sigma_deps=NULL; // ici comme d_sigma_deps = NULL, il n'y aura pas de calcul explicite de d_sigma_deps // mais seulement le calcul de T_d_Mpres_d_V Hysteresis_bulk::Dsig_depsilon(T_d_Mpres_d_V,en_base_orthonormee,gijHH,d_sigma_deps,V_tdt); // ==== 3 === calcul de la variation de contrainte / au ddl // nombre de ddl pour la variation du tenseur des contraintes int nbddl = d_gijBB_tdt.Taille(); Tenseur3BH depsBH_dll; for (int i = 1; i<= nbddl; i++) { // on fait uniquement une égalité d'adresse de manière à ne pas utiliser // le constructeur d'ou la profusion d'* et de () Tenseur3HH & dsigHH = *((Tenseur3HH*) (d_sigHH(i))); // passage en dim 3 const Tenseur3HH & dgijHH = *((Tenseur3HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture // compte tenu des relations permettant le calcul de delta_V // on a : d_delta_V = d_jacobien_tdt / (*ex.jacobien_0) double d_delta_V = (*ex.d_jacobien_tdt)(i) * un_sur_jaco_0; // d'autre part sachant que : sigHH = gijHH * (- pression_t___tdt); // et que T_d_Mpres_d_V = d_Mpression_t___tdt dsigHH = dgijHH * (MPr_t___tdt) + gijHH * (T_d_Mpres_d_V) * d_delta_V; }; LibereTenseur(); }; // calcul des contraintes et ses variations par rapport aux déformations a t+dt // en_base_orthonormee: le tenseur de contrainte en entrée est en orthonormee // le tenseur de déformation et son incrémentsont également en orthonormees // si = false: les bases transmises sont utilisées // ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a void Hysteresis_bulk::Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & sigHH_t,TenseurBB& ,TenseurBB & epsBB_tdt,TenseurBB & delta_epsBB_,double& ,double& ,TenseurHH& sigHH_tdt,TenseurHHHH& d_sigma_deps ,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement ,const Met_abstraite::Umat_cont& ex) { #ifdef MISE_AU_POINT if (epsBB_tdt.Dimension() != 3) { cout << "\nErreur : la dimension devrait etre 3 !\n"; cout << " Hysteresis_bulk::Calcul_dsigma_deps\n"; Sortie(1); }; #endif const Tenseur3HH & gijHH = *((Tenseur3HH*) ex.gijHH_tdt); // " " " " const Tenseur3BB & gijBB = *((Tenseur3BB*) ex.gijBB_tdt); // " " " " const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3 const Tenseur3BB & delta_epsBB = *((Tenseur3BB*) &delta_epsBB_); // passage en dim 3 Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // " " " " // Tenseur3HH & sigHH_i = *((Tenseur3HH*) &sigHH_t); // " " " " SaveResulHysteresis_bulk & save_resul = *((SaveResulHysteresis_bulk*) saveResul); // initialisation des variables de travail save_resul.Init_debut_calcul(); // initialisation éventuelle des variables thermo-dépendantes Init_thermo_dependance(); // calcul de delta_V : incrément de variation relative de volume double V_tdt = (*ex.jacobien_tdt)/(*ex.jacobien_0); double V_t = (*ex.jacobien_t)/(*ex.jacobien_0); delta_V = V_tdt - V_t; // dans le cas de la prise en compte de la dilatation il faut changer V et delta_V if (dilatation) // on considère que la variation de volume thermique = la trace du tenseur de def thermique {double var_vol_thermique = (gijHH * (*epsBB_therm)).Trace(); V_tdt -= var_vol_thermique; // recup de l'incrément de temps double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant(); // on calcul un incrément en utilisant la vitesse de déformation thermique double delta_var_vol_thermique = (gijHH * (*DepsBB_therm)).Trace()*deltat; delta_V -= delta_var_vol_thermique; }; // === 1 == calcul de l'avancement temporel sur 1 pas, // ici il n'y a pas de différence entre orthonormee ou pas au niveau du type de calcul int cas = 1; // donc cas avec dérivée de Jauman en mixte if (en_base_orthonormee) Avancement_temporel(IdHH3,IdBB3,cas,save_resul,sigHH,energ_t,energ); else Avancement_temporel(gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ); // ==== 2 === calcul de l'opérateur tangent ============= double T_d_pres_d_V; Hysteresis_bulk::Dsig_depsilon(T_d_pres_d_V,en_base_orthonormee,gijHH,&d_sigma_deps,V_tdt); // on libère les tenseurs intermédiaires LibereTenseur();LibereTenseurQ(); }; //----- lecture écriture de restart ----- // cas donne le niveau de la récupération // = 1 : on récupère tout // = 2 : on récupère uniquement les données variables (supposées comme telles) void Hysteresis_bulk::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D ,LesFonctions_nD& lesFonctionsnD) {string toto; if (cas == 1) { bool test;string nom; // param_prager_(np) ent >> nom >> test; if (!test) { ent >> xnp; if (xnp_temperature != NULL) {if (xnp_temperature->NomCourbe() == "_") delete xnp_temperature; xnp_temperature = NULL;}; } else { ent >> nom; xnp_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,xnp_temperature); }; // limite_plasticite_(Q0) ent >> nom >> test; if (!test) { ent >> Qzero; if (Qzero_temperature != NULL) {if (Qzero_temperature->NomCourbe() == "_") delete Qzero_temperature; Qzero_temperature = NULL;}; } else { ent >> nom; Qzero_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,Qzero_temperature); }; // param_lame_(mu) ent >> nom >> test; if (!test) { ent >> xmu; if (xmu_temperature != NULL) {if (xmu_temperature->NomCourbe() == "_") delete xmu_temperature; xmu_temperature = NULL;}; } else { ent >> nom; xmu_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,xmu_temperature); }; // les paramètres de réglage de l'algorithme ent >> nom >> tolerance_residu >> nom >> tolerance_residu_rel >> nom >> nb_boucle_maxi >> nom >> maxi_delta_var_sig_sur_iter_pour_Newton >> nom >> tolerance_coincidence; ent >> nom >> type_resolution_equa_constitutive >> nom >> cas_kutta >> nom >> erreurAbsolue >> nom >> erreurRelative; ent >> nom >> nb_maxInvCoinSurUnPas; ent >> nom >> depassement_Q0; // le niveau d'affichage Lecture_permet_affichage(ent,cas,lesFonctionsnD); ent >> nom >> sortie_post; }; // appel de la classe mère Loi_comp_abstraite::Lecture_don_base_info(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD); }; // cas donne le niveau de sauvegarde // = 1 : on sauvegarde tout // = 2 : on sauvegarde uniquement les données variables (supposées comme telles) void Hysteresis_bulk::Ecriture_base_info_loi(ofstream& sort,const int cas) { if (cas == 1) { sort << " HYSTERESIS_BULK "; sort << " param_prager_(np) "; if (xnp_temperature == NULL) { sort << false << " " << xnp << " ";} else { sort << true << " fonction_xnp_temperature "; LesCourbes1D::Ecriture_pour_base_info(sort,cas,xnp_temperature); }; sort << " limite_plasticite_(Q0) "; if (Qzero_temperature == NULL) { sort << false << " " << Qzero << " ";} else { sort << true << " fonction_Qzero_temperature "; LesCourbes1D::Ecriture_pour_base_info(sort,cas,Qzero_temperature); }; sort << " param_xmu_(mu) "; if (xmu_temperature == NULL) { sort << false << " " << xmu << " ";} else { sort << true << " fonction_xmu_temperature "; LesCourbes1D::Ecriture_pour_base_info(sort,cas,xmu_temperature); }; sort << "\n tolerance_algo_newton_equadiff:_absolue= " << tolerance_residu << " relative= "<< tolerance_residu_rel << " nb_boucle_maxi " << nb_boucle_maxi << " maxi_delta_var_sig_sur_iter_pour_Newton " << maxi_delta_var_sig_sur_iter_pour_Newton << " tolerance_coincidence " << tolerance_coincidence; sort << "\n type_resolution_equa_constitutive " << type_resolution_equa_constitutive << " cas_kutta " << cas_kutta << " erreurAbsolue " << erreurAbsolue << " erreurRelative " << erreurRelative << " "; sort << " nb_maxInvCoinSurUnPas " << nb_maxInvCoinSurUnPas << " "; sort << " depassement_Q0 " << depassement_Q0 << " "; // niveau d'affichage Affiche_niveau_affichage(sort,cas); sort << " sortie_post " << sortie_post << " "; }; // appel de la classe mère Loi_comp_abstraite::Ecriture_don_base_info(sort,cas); }; // affichage des informations pour le débug void Hysteresis_bulk::Affiche_debug(double& unSur_wBaseCarre,SaveResulHysteresis_bulk & save_resul ,List_io ::iterator& iatens,bool& pt_sur_principal ,const List_io ::iterator& iatens_princ ) { // affichage de la loi Affiche(); // départ cout << "\n -- valeurs initiales -- "; cout << "\n MPr_t= " << save_resul.MPr_t; // à t // variables actuelles cout << "\n -- variables actuelles -- "; cout << "\n MPr_tdt= " << save_resul.MPr_tdt; // valeur finale cout << "\n MPr_i___= " << MPr_i___; // MPr de début de calcul (à t au début) cout << "\n MPr_R= " << MPr_R; // MPr de Référence en cours cout << "\n delta_MPr_Rat= " << delta_MPr_Rat; // deltat MPr de R a t cout << "\n delta_MPr_Ratdt= " << delta_MPr_Ratdt; // deltat MPr de R a tdt cout << "\n delta_MPr_tatdt= " << delta_MPr_tatdt; // delta MPr de t à tdt cout << "\n unSur_wBaseCarre= " << unSur_wBaseCarre; cout << "\n nb_coincidence= " << save_resul.nb_coincidence; cout << "\n modif= " << save_resul.modif; // variables de liste cout << "\n -- variables de liste -- "; cout << "\n *iatens= " << *iatens; cout << "\n pt_sur_principal= " << pt_sur_principal; cout << "\n *iatens_princ= " << *iatens_princ; // les listes cout << "\n -- les listes -- "; cout << "\n save_resul.MPr_R_t_a_tdt= " << save_resul.MPr_R_t_a_tdt; cout << "\n save_resul.fct_aide_t_a_tdt= " << save_resul.fct_aide_t_a_tdt; cout << "\n list indic_coin: " << save_resul.indic_coin; cout << "\n list MPr_R: " << save_resul.MPr_R; cout << "\n list fct_aide: " << save_resul.fct_aide << " "; // on vide le buffer cout << endl; };