// FICHIER : Maheo_hyper.cc // CLASSE : Maheo_hyper // 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 "TypeConsTens.h" #include "CharUtil.h" #include "Maheo_hyper.h" #include "MathUtil.h" #include "Enum_TypeQuelconque.h" #include "TypeQuelconqueParticulier.h" Maheo_hyper::Maheo_hyper () : // Constructeur par defaut Hyper_W_gene_3D(MAHEO_HYPER,CAT_MECANIQUE,3) ,Qsig_rev(-ConstMath::trespetit) , mu1_rev(-ConstMath::trespetit) ,mu2_rev(-ConstMath::trespetit) ,Qsig_use(0.), mu1_use(0.),mu2_use(0.) // recalculés avant chaque utilisation ,Qsig_nD(NULL),mu1_nD(NULL),mu2_nD(NULL) ,fact_regularisation(ConstMath::pasmalpetit) ,Qsig_rev_temperature(NULL),mu1_rev_temperature(NULL),mu2_rev_temperature(NULL) ,W_d(0.),W_v(0.) ,W_d_J1(0.),W_d_J2(0.),W_v_J3(0.),W_v_J3J3(0.) ,W_d_J1_2(0.),W_d_J1_J2(0.),W_d_J2_2(0.) {nb_para_loi = 3; }; // Constructeur de copie Maheo_hyper::Maheo_hyper (const Maheo_hyper& loi) : Hyper_W_gene_3D(loi) ,Qsig_rev(loi.Qsig_rev) , mu1_rev(loi.mu1_rev) ,mu2_rev(loi.mu2_rev) ,Qsig_use(0.), mu1_use(0.),mu2_use(0.) // recalculés avant chaque utilisation ,Qsig_rev_temperature(loi.Qsig_rev_temperature) ,mu1_rev_temperature(loi.mu1_rev_temperature) ,mu2_rev_temperature(loi.mu2_rev_temperature) ,Qsig_nD(loi.Qsig_nD),mu1_nD(loi.mu1_nD),mu2_nD(loi.mu2_nD) ,fact_regularisation(loi.fact_regularisation) ,W_d(loi.W_d),W_v(loi.W_v) ,W_d_J1(loi.W_d_J1),W_d_J2(loi.W_d_J2),W_v_J3(loi.W_v_J3),W_v_J3J3(loi.W_v_J3J3) ,W_d_J1_2(loi.W_d_J1_2),W_d_J1_J2(loi.W_d_J1_J2),W_d_J2_2(loi.W_d_J2_2) { nb_para_loi = 3; // nb paramètres de la loi // on regarde s'il s'agit de courbes locales ou de courbes globales if (Qsig_rev_temperature != NULL) if (Qsig_rev_temperature->NomCourbe() == "_") Qsig_rev_temperature = Courbe1D::New_Courbe1D(*(loi.Qsig_rev_temperature)); if (mu1_rev_temperature != NULL) if (mu1_rev_temperature->NomCourbe() == "_") mu1_rev_temperature = Courbe1D::New_Courbe1D(*(loi.mu1_rev_temperature)); if (mu2_rev_temperature != NULL) if (mu2_rev_temperature->NomCourbe() == "_") mu2_rev_temperature = Courbe1D::New_Courbe1D(*(loi.mu2_rev_temperature)); // idem pour les fonctions nD if (Qsig_nD != NULL) if (Qsig_nD->NomFonction() == "_") {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi) string non_courbe("_"); Qsig_nD = Fonction_nD::New_Fonction_nD(*loi.Qsig_nD); }; if (mu1_nD != NULL) if (mu1_nD->NomFonction() == "_") {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi) string non_courbe("_"); mu1_nD = Fonction_nD::New_Fonction_nD(*loi.mu1_nD); }; if (mu2_nD != NULL) if (mu2_nD->NomFonction() == "_") {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi) string non_courbe("_"); mu2_nD = Fonction_nD::New_Fonction_nD(*loi.mu2_nD); }; }; Maheo_hyper::~Maheo_hyper () // Destructeur { if (Qsig_rev_temperature != NULL) if (Qsig_rev_temperature->NomCourbe() == "_") delete Qsig_rev_temperature; if (mu1_rev_temperature != NULL) if (mu1_rev_temperature->NomCourbe() == "_") delete mu1_rev_temperature; if (mu2_rev_temperature != NULL) if (mu2_rev_temperature->NomCourbe() == "_") delete mu2_rev_temperature; // idem pour les fonctions nD if (Qsig_nD != NULL) if (Qsig_nD->NomFonction() == "_") delete Qsig_nD; if (mu1_nD != NULL) if (mu1_nD->NomFonction() == "_") delete mu1_nD; if (mu2_nD != NULL) if (mu2_nD->NomFonction() == "_") delete mu2_nD; }; // Lecture des donnees de la classe sur fichier void Maheo_hyper::LectureDonneesParticulieres (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D ,LesFonctions_nD& lesFonctionsnD) { string toto,nom; // lecture de Qsig_rev string mot_cle = "Qsig_rev="; string nom_class_methode = "Maheo_hyper::LectureDonneesParticulieres"; entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle); // maintenant choix thermodépendant ou pas if(strstr(entreePrinc->tablcar,"Qsig_rev_thermo_dependant_")!=0) { mot_cle = "Qsig_rev_thermo_dependant_"; entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle); // lecture de la loi d'évolution Qsig_rev_temperature = Lecture_courbe(entreePrinc,lesCourbes1D); thermo_dependant=true; // on indique que la loi est thermo-dépendant } else *(entreePrinc->entree) >> Qsig_rev ; // lecture de mu1_rev mot_cle = "mu1_rev="; entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle); // maintenant choix thermodépendant ou pas if(strstr(entreePrinc->tablcar,"mu1_rev_temperature_")!=0) { mot_cle = "mu1_rev_temperature_"; entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle); // lecture de la loi d'évolution mu1_rev_temperature = Lecture_courbe(entreePrinc,lesCourbes1D); thermo_dependant=true; // on indique que la loi est thermo-dépendant } else *(entreePrinc->entree) >> mu1_rev; // lecture de mu2_rev mot_cle = "mu2_rev="; entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle); // maintenant choix thermodépendant ou pas if(strstr(entreePrinc->tablcar,"mu2_rev_temperature_")!=0) { mot_cle = "mu2_rev_temperature_"; entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle); // lecture de la loi d'évolution mu2_rev_temperature = Lecture_courbe(entreePrinc,lesCourbes1D); thermo_dependant=true; // on indique que la loi est thermo-dépendant } else {*(entreePrinc->entree) >> mu2_rev; // s'il n'y a plus rien n'a lire, il faut passer un enregistrement if((strstr(entreePrinc->tablcar,"avec_regularisation_")==0) && (strstr(entreePrinc->tablcar,"avec_nD_")==0) ) entreePrinc->NouvelleDonnee(); }; // -------- on regarde s'il y a une dépendance à une fonction nD if (strstr(entreePrinc->tablcar,"avec_nD_")!=NULL) { string nom1,nom2; // lecture des fonctions nD entreePrinc->NouvelleDonnee(); // passage d'une ligne // on lit tant que l'on ne rencontre pas la ligne contenant "fin_parametres_avec_nD_" // ou un nouveau mot clé global auquel cas il y a pb !! MotCle motCle; // ref aux mots cle while (strstr(entreePrinc->tablcar,"fin_parametres_avec_nD_")==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 dependance a une fonction nD: on n'a pas trouve le mot cle " << " fin_parametres_avec_nD_ et par contre la ligne courante contient un mot cle global "; entreePrinc->MessageBuffer ("** erreur51 des parametres de reglage de la loi de comportement Maheo_hyper**"); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; // lecture d'un mot clé *(entreePrinc->entree) >> nom1; 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) >>nom1; } #else else if ((entreePrinc->entree)->eof()) // la lecture est bonne mais on a atteind la fin de la ligne { if(nom1 != "fin_parametres_avec_nD_") {entreePrinc->NouvelleDonnee(); *(entreePrinc->entree) >> nom1;}; } #endif else // cas d'une erreur de lecture { cout << "\n erreur de lecture inconnue "; entreePrinc->MessageBuffer ("** erreur61 des parametres de reglage de la loi de comportement Maheo_hyper**"); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; if(nom1 == "Qsig_nD=") {string nom_fonct; // lecture du nom de la fonction *(entreePrinc->entree) >> nom_fonct; // maintenant on définit la fonction if (lesFonctionsnD.Existe(nom_fonct)) {Qsig_nD = lesFonctionsnD.Trouve(nom_fonct); } else {// sinon il faut la lire maintenant string non("_"); Qsig_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct)); // lecture de la fonction Qsig_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc); // maintenant on vérifie que la fonction est utilisable if (Qsig_nD->NbComposante() != 1 ) { cout << "\n erreur en lecture, la fonction " << nom_fonct << " est une fonction vectorielle a " << Qsig_nD->NbComposante() << " composante alors qu'elle devrait etre scalaire ! " << " elle n'est donc pas utilisable !! "; string message("\n**erreur031** \n"+nom_class_methode+"(..."); entreePrinc->MessageBuffer(message); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; }; // on regarde si la fonction nD intègre la température const Tableau & tab_enu = Qsig_nD->Tab_enu_etendu(); if (tab_enu.Contient(TEMP)) thermo_dependant=true; // prepa du flot de lecture entreePrinc->NouvelleDonnee(); } else if(nom1 == "mu1_nD=") {string nom_fonct; // lecture du nom de la fonction *(entreePrinc->entree) >> nom_fonct; // maintenant on définit la fonction if (lesFonctionsnD.Existe(nom_fonct)) {mu1_nD = lesFonctionsnD.Trouve(nom_fonct); } else {// sinon il faut la lire maintenant string non("_"); mu1_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct)); // lecture de la fonction mu1_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc); // maintenant on vérifie que la fonction est utilisable if (mu1_nD->NbComposante() != 1 ) { cout << "\n erreur en lecture, la fonction " << nom_fonct << " est une fonction vectorielle a " << mu1_nD->NbComposante() << " composante alors qu'elle devrait etre scalaire ! " << " elle n'est donc pas utilisable !! "; string message("\n**erreur031** \n"+nom_class_methode+"(..."); entreePrinc->MessageBuffer(message); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; }; // on regarde si la fonction nD intègre la température const Tableau & tab_enu = mu1_nD->Tab_enu_etendu(); if (tab_enu.Contient(TEMP)) thermo_dependant=true; // prepa du flot de lecture entreePrinc->NouvelleDonnee(); } else if(nom1 == "mu2_nD=") {string nom_fonct; // lecture du nom de la fonction *(entreePrinc->entree) >> nom_fonct; // maintenant on définit la fonction if (lesFonctionsnD.Existe(nom_fonct)) {mu2_nD = lesFonctionsnD.Trouve(nom_fonct); } else {// sinon il faut la lire maintenant string non("_"); mu2_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct)); // lecture de la fonction mu2_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc); // maintenant on vérifie que la fonction est utilisable if (mu2_nD->NbComposante() != 1 ) { cout << "\n erreur en lecture, la fonction " << nom_fonct << " est une fonction vectorielle a " << mu2_nD->NbComposante() << " composante alors qu'elle devrait etre scalaire ! " << " elle n'est donc pas utilisable !! "; string message("\n**erreur031** \n"+nom_class_methode+"(..."); entreePrinc->MessageBuffer(message); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; }; // on regarde si la fonction nD intègre la température const Tableau & tab_enu = mu2_nD->Tab_enu_etendu(); if (tab_enu.Contient(TEMP)) thermo_dependant=true; // prepa du flot de lecture entreePrinc->NouvelleDonnee(); } else // cas d'une erreur { cout << "\n erreur en lecture de la dependance a une fonction nD des coefficients, on aurait du lire " << " le mot cle Qsig_nD= ou mu1_nD= ou mu2_nD= et on a lu: " << nom1; entreePrinc->MessageBuffer("**erreur7 Maheo_hyper::LectureDonneesParticulieres (...**"); throw (UtilLecture::ErrNouvelleDonnee(-1)); Sortie(1); }; }; //-- fin du while entreePrinc->NouvelleDonnee(); // passage du mot clé "fin_parametres_avec_nD_" }; // cas avec une régularisation éventuelle if (strstr(entreePrinc->tablcar,"avec_regularisation_")!=NULL) { mot_cle = "avec_regularisation_"; entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle); *(entreePrinc->entree) >> fact_regularisation; }; // lecture de l'indication éventuelle du post traitement string le_mot_cle = "sortie_post_"; if (entreePrinc->Lecture_un_parametre_int(0,nom_class_methode,0,1,le_mot_cle,sortie_post)) entreePrinc->NouvelleDonnee(); // passage éventuel du mot clé sortie_post_; // --- appel au niveau de la classe mère // ici il n'y a pas de type de déformation associé // mais on prend la def standart d'almansi, pour les fonctions associées éventuelles Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire (*entreePrinc,lesFonctionsnD,true); }; // affichage de la loi void Maheo_hyper::Affiche() const { cout << " \n loi de comportement hyper elastique 3D Maheo_hyper \n"; if ( Qsig_rev_temperature != NULL) { cout << " Qsig_rev thermo dependant " << " courbe Qsig_rev=f(T): " << Qsig_rev_temperature->NomCourbe() <<" ";} else { cout << " Qsig_rev= " << Qsig_rev << " ";}; if ( mu1_rev_temperature != NULL) { cout << " mu1_rev thermo dependant " << " courbe mu1_rev=f(T): " << mu1_rev_temperature->NomCourbe() <<" ";} else { cout << " mu1_rev= " << mu1_rev << " ";}; if ( mu2_rev_temperature != NULL) { cout << " mu2_rev thermo dependant " << " courbe mu2_rev=f(T): " << mu2_rev_temperature->NomCourbe() <<" ";} else { cout << " mu2_rev= " << mu2_rev << " ";}; // les dépendances à une fonction nD if (Qsig_nD) {cout << " fonction nD Qsig_nD: "; if (Qsig_nD->NomFonction() != "_") cout << Qsig_nD->NomFonction(); else Qsig_nD->Affiche(); }; if (mu1_nD) {cout << " fonction nD mu1_nD: "; if (mu1_nD->NomFonction() != "_") cout << mu1_nD->NomFonction(); else mu1_nD->Affiche(); }; if (mu2_nD) {cout << " fonction nD mu2_nD: "; if (mu2_nD->NomFonction() != "_") cout << mu2_nD->NomFonction(); else mu2_nD->Affiche(); }; // appel de la classe mère Loi_comp_abstraite::Affiche_don_classe_abstraite(); }; // affichage et definition interactive des commandes particulières à chaques lois void Maheo_hyper::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"); sort << "\n# ----------------------------------------------------------------------------------" << "\n# |... loi de comportement hyper elastique 3D Maheo_hyper .. |" << "\n# ----------------------------------------------------------------------------------" << "\n\n# exemple de definition de loi" << "\n Qsig_rev= 2 mu1_rev= 0.0167 mu2_rev= 0.45 " << "\n# .. fin de la definition de la loi Maheo_hyper \n" << endl; if ((rep != "o") && (rep != "O" ) && (rep != "0") ) { sort << "\n# " << "\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# Qsig_rev= Qsig_rev_thermo_dependant_ courbe1 " << "\n# mu1_rev= mu1_rev_temperature_ courbe2 " << "\n# mu2_rev= mu2_rev_temperature_ courbe3 " << "\n# ou encore " << "\n# Qsig_rev= Qsig_rev_temperature_ courbe2 " << "\n# mu1_rev= 0.0167 mu2_rev= mu2_rev_temperature_ courbe4 " << "\n#------------------------------------------------------------------------------------" << "\n# noter qu'apres la definition de chaque courbe, on change de ligne, a l'inverse " << "\n# si la valeur du parametre est fixe, on poursuit sur la meme ligne. " << "\n#" << "\n#--------------------- dependance explicite a une fonction nD -----------------------------------" << "\n# il est egalement possible de definir une dependance des parametres a une fonction nD " << "\n# dans ce cas, a la suite de la dependance explicite a la temperature via une courbr 1D " << "\n# ou a la suite du dernier parametre fixe, on met le mot cle: avec_nD_ " << "\n# puis sur les lignes qui suivent on met des fonctions nD multiplicatives des parametres " << "\n# initiaux, puis un mot cle signalant la fin des fonctions: fin_parametres_avec_nD_ " << "\n# suit un exemple de declaration (il est possible d'ommettre certaine fonction ! " << "\n# la presence de fct nD est facultative), de meme il est possible soit d'utiliser une fonction " << "\n# deja defini soit de donner directement la fonction nD (comme pour les autres lois) " << "\n# Enfin, apres chaque definition de fonction il faut passer a une nouvelle ligne " << "\n# " << "\n# exemple de definition de loi" << "\n# Qsig_rev= 1 mu1_rev= 1 mu2_rev= 1 avec_nD_ " << "\n# Qsig_nD= fct_nD_1 " << "\n# mu1_nD= fct_nD_2 " << "\n# mu2_nD= fct_nD_3 " << "\n# fin_parametres_avec_nD_ " << "\n# " << "\n# Remarques: " << "\n# - les fonction nD sont des fonctions multiplicatives " << "\n# des parametres initiaux (ou de ceux obtenus apres dependance explicite a la temperature " << "\n# ( bien noter que la loi obtenue peut-etre quelconque, en particulier plus du tout " << "\n# hyperelastique, tout depend des choix des fcts nD !) " << "\n# " << "\n#----------------------------------- avec_regularisation_ -------------------------------" << "\n# il est possible d'indiquer un facteur de regularisation qui permet d'eviter " << "\n# des problemes de NaN, de type division par 0 par exemple " << "\n# 1/a est remplace par 1/(a+fact_regularisation), par defaut fact_regularisation = 1.e-12 " << "\n# pour indiquer un facteur de regulation non nul on indique en dernier parametre " << "\n# le mot cle avec_regularisation_ suivi du facteur voulu " << "\n# ex: " << "\n# avec_regularisation_ 1.e-12 " << "\n# ce mot cle doit se situer avant eventuelement le mot cle sortie_post_ " << "\n#" << "\n#----------------------------------- sortie_post_ -------------------------------------" << "\n# il est possible de recuperer differentes grandeurs de travail par exemple " << "\n# l'intensite du potentiel, comme ces grandeurs sont calculees au moment de la resolution " << "\n# et ne sont pas stockees. On peut egalement recuperer la valeur des coefficients de la loi " << "\n# (utile s'ils varient) sous forme d'un tableau de triplet : Qsig_rev C01 C10 mu1_rev mu2_use " << "\n# Il faut indiquer le mot cle sortie_post_ suivi de 1 (par defaut = 0) " << "\n# ensuite au moment de la constitution du .CVisu on aura acces aux grandeurs de travail " << "\n# ex: " << "\n# sortie_post_ 1 " << "\n# ce mot cle est le dernier des parametres specifiques de la loi il doit se situe " << "\n# a la fin de la derniere ligne de donnees " << "\n#" << "\n#------------------------------------------------------------------------------------" << endl; }; // appel de la classe Hyper_W_gene_3D Hyper_W_gene_3D::Info_commande_LoisDeComp_hyper3D(entreePrinc); // appel de la classe mère Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc); }; // test si la loi est complete int Maheo_hyper::TestComplet() { int ret = LoiAbstraiteGeneral::TestComplet(); if ((Qsig_rev_temperature == NULL) && (Qsig_rev == (-ConstMath::trespetit))) { cout << " \n le coefficient Qsig_rev de la loi Maheo_hyper n'est pas defini pour la loi " << Nom_comp(id_comp) << '\n' << endl; ret = 0; }; if ((mu1_rev_temperature == NULL) && (mu1_rev == (-ConstMath::trespetit))) { cout << " \n le coefficient mu1_rev de la loi Maheo_hyper n'est pas defini pour la loi " << Nom_comp(id_comp) << '\n' << endl; ret = 0; }; if ((mu2_rev_temperature == NULL) && (mu2_rev == (-ConstMath::trespetit))) { cout << " \n le coefficient mu2_rev de la loi Maheo_hyper n'est pas defini pour la loi " << Nom_comp(id_comp) << '\n' << endl; ret = 0; }; // if (ret == 0) {this-> Affiche(); ret = 0; }; return ret; }; //----- 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 Maheo_hyper::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D ,LesFonctions_nD& lesFonctionsnD) { string nom; bool test; if (cas == 1) { // === Qsig_rev ent >> nom >> test; // vérification #ifdef MISE_AU_POINT if (nom != "Qsig_rev") { cout << "\n erreur en lecture du parametres Qsig_rev de la loi de mooney rivlin 3D" << " on devait lire Qsig_rev= avant le second parametre " << " et on a lue: " << nom << " " << "\n Maheo_hyper::Lecture_base_info_loi(..." << endl; Sortie(1); } #endif if (!test) { ent >> Qsig_rev; if (Qsig_rev_temperature != NULL) {if (Qsig_rev_temperature->NomCourbe() == "_") delete Qsig_rev_temperature; Qsig_rev_temperature = NULL;}; } else { ent >> nom; Qsig_rev_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,Qsig_rev_temperature); }; // on regarde si il y a une fonction nD multiplicative ent >> nom; if (nom == "Qsig_nD:") Qsig_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,Qsig_nD); // === mu1_rev ent >> nom >> test; // vérification #ifdef MISE_AU_POINT if (nom != "mu1_rev") { cout << "\n erreur en lecture du parametres mu1_rev de la loi de mooney rivlin 3D" << " on devait lire mu1_rev= avant le second parametre " << " et on a lue: " << nom << " " << "\n Maheo_hyper::Lecture_base_info_loi(..." << endl; Sortie(1); } #endif if (!test) { ent >> mu1_rev; if (mu1_rev_temperature != NULL) {if (mu1_rev_temperature->NomCourbe() == "_") delete mu1_rev_temperature; mu1_rev_temperature = NULL;}; } else { ent >> nom; mu1_rev_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,mu1_rev_temperature); }; // on regarde si il y a une fonction nD multiplicative ent >> nom; if (nom == "mu1_nD:") mu1_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,mu1_nD); // === mu2_rev ent >> nom >> test; // vérification #ifdef MISE_AU_POINT if (nom != "mu2_rev") { cout << "\n erreur en lecture du parametres mu2_rev de la loi de mooney rivlin 3D" << " on devait lire mu2_rev= avant le second parametre " << " et on a lue: " << nom << " " << "\n Maheo_hyper::Lecture_base_info_loi(..." << endl; Sortie(1); } #endif if (!test) { ent >> mu2_rev; if (mu2_rev_temperature != NULL) {if (mu2_rev_temperature->NomCourbe() == "_") delete mu2_rev_temperature; mu2_rev_temperature = NULL;}; } else { ent >> nom; mu2_rev_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,mu2_rev_temperature); }; // on regarde si il y a une fonction nD multiplicative ent >> nom; if (nom == "mu2_nD:") mu2_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,mu2_nD); }; // 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 Maheo_hyper::Ecriture_base_info_loi(ofstream& sort,const int cas) { if (cas == 1) { // Qsig_rev sort << " Qsig_rev= "; if (Qsig_rev_temperature == NULL) { sort << false << " " << Qsig_rev << " ";} else { sort << true << " fonction_Qsig_rev_temperature "; LesCourbes1D::Ecriture_pour_base_info(sort,cas,Qsig_rev_temperature); }; if (Qsig_nD != NULL) {sort << " Qsig_nD: " << " "; LesFonctions_nD::Ecriture_pour_base_info(sort, cas,Qsig_nD); } else {sort << " non_Qsig_nD "; }; // mu1_rev sort << " mu1_rev= "; if (mu1_rev_temperature == NULL) { sort << false << " " << mu1_rev << " ";} else { sort << true << " fonction_mu1_rev_temperature "; LesCourbes1D::Ecriture_pour_base_info(sort,cas,mu1_rev_temperature); }; if (mu1_nD != NULL) {sort << " 2 mu1_nD: " << " "; LesFonctions_nD::Ecriture_pour_base_info(sort, cas,mu1_nD); } else {sort << " non_mu1_nD "; }; // mu2_rev sort << " mu2_rev= "; if (mu2_rev_temperature == NULL) { sort << false << " " << mu2_rev << " ";} else { sort << true << " fonction_mu2_rev_temperature "; LesCourbes1D::Ecriture_pour_base_info(sort,cas,mu2_rev_temperature); }; if (mu2_nD != NULL) {sort << " 2 mu1_nD: " << " "; LesFonctions_nD::Ecriture_pour_base_info(sort, cas,mu1_nD); } else {sort << " non_mu2_nD "; }; }; // appel de la classe mère Loi_comp_abstraite::Ecriture_don_base_info(sort,cas); }; // calcul d'un module d'young équivalent à la loi: pour l'instant à faire !! double Maheo_hyper::Module_young_equivalent(Enum_dure temps,const Deformation & ,SaveResul * ) {// comme la loi n'a pas de partie volumique, le module d'Young équivalent est nul return 0.; }; // calcul d'un module compressibilité équivalent à la loi: double Maheo_hyper::Module_compressibilite_equivalent(Enum_dure temps,const Deformation & ,SaveResul * ) {// comme la loi n'a pas de partie volumique, le module de compressibilité est nul return 0.; }; // ========== codage des METHODES VIRTUELLES protegees:================ // calcul des contraintes a t+dt void Maheo_hyper::Calcul_SigmaHH (TenseurHH& ,TenseurBB& ,DdlElement & tab_ddl, TenseurBB & ,TenseurHH & ,BaseB& ,BaseH& ,TenseurBB & epsBB_, TenseurBB & , TenseurBB & gijBB_,TenseurHH & gijHH_,Tableau & d_gijBB_, double& jacobien_0,double& jacobien,TenseurHH & sigHH_ ,EnergieMeca & energ,const EnergieMeca & ,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 << " Maheo_hyper::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 << " Maheo_hyper::Calcul_SigmaHH\n"; Sortie(1); }; #endif Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // passage explicite en tenseur dim 3 // init coef matériaux Qsig_use = Qsig_rev; mu1_use = mu1_rev; mu2_use = mu2_rev; // cas de la thermo dépendance, on calcul les grandeurs if (Qsig_rev_temperature != NULL) Qsig_use = Qsig_rev_temperature->Valeur(*temperature); if (mu1_rev_temperature != NULL) mu1_use = mu1_rev_temperature->Valeur(*temperature); if (mu2_rev_temperature != NULL) mu2_use = mu2_rev_temperature->Valeur(*temperature); // récup du conteneur spécifique du point, pour sauvegarde éventuelle SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul); // dans le cas de l'utilisation de fct nD if ( (Qsig_nD != NULL) || (mu1_nD != NULL) || (mu2_nD != NULL) ) { // opération de transmission de la métrique const Met_abstraite::Impli* ex_impli = NULL; const Met_abstraite::Expli_t_tdt* ex_expli_tdt = &ex; const Met_abstraite::Umat_cont* ex_expli = NULL; if (Qsig_nD != NULL) // là il faut calculer la fonction nD { // on utilise la méthode générique de loi abstraite list list_save; // inter pour l'appel de la fonction list_save.push_back(saveResul); Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee (Qsig_nD,1 // une seule valeur attendue en retour ,ex_impli,ex_expli_tdt,ex_expli ,NULL ,NULL ,&list_save ); // on récupère le premier élément du tableau uniquement Qsig_use *= tab_val(1); }; if (mu1_nD != NULL) // là il faut calculer la fonction nD { // on utilise la méthode générique de loi abstraite list list_save; // inter pour l'appel de la fonction list_save.push_back(saveResul); Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee (mu1_nD,1 // une seule valeur attendue en retour ,ex_impli,ex_expli_tdt,ex_expli ,NULL ,NULL ,&list_save ); // on récupère le premier élément du tableau uniquement mu1_use *= tab_val(1); }; if (mu2_nD != NULL) // là il faut calculer la fonction nD { // on utilise la méthode générique de loi abstraite list list_save; // inter pour l'appel de la fonction list_save.push_back(saveResul); Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee (mu2_nD,1 // une seule valeur attendue en retour ,ex_impli,ex_expli_tdt,ex_expli ,NULL ,NULL ,&list_save ); // on récupère le premier élément du tableau uniquement mu2_use *= tab_val(1); }; }; if (sortie_post) {// sauvegarde des paramètres matériau (*save_resulHyper_W.para_loi)(1) = Qsig_use; (*save_resulHyper_W.para_loi)(2) = mu1_use; (*save_resulHyper_W.para_loi)(3) = mu2_use; }; // calcul des invariants et de leurs variations premières (méthode de Hyper_W_gene_3D) Invariants_et_var1(*(ex.gijBB_0),*(ex.gijHH_0),gijBB_,gijHH_,jacobien_0,jacobien); // calcul du potentiel et de ses dérivées premières / aux invariants J_r Potentiel_et_var(); // stockage éventuel pour du post-traitement if (sortie_post) { save_resulHyper_W.invP->potentiel= W_d + W_v; }; // calcul du tenseur des contraintes: on tiend compte du fait que le potentiel ne dépend que de J1 sigHH = (W_d_J1/V) * d_J_r_epsBB_HH(1); // sigHH = (W_d_J1/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2) // + (W_v_J3/V) * d_J_r_epsBB_HH(3); // calcul du module de compressibilité module_compressibilite = 2. * W_v_J3; // ici 0 // pour le module de cisaillement, pour l'instant je ne fais rien !! à voir *** module_cisaillement = 0.; // traitement des énergies energ.Inita(0.); energ.ChangeEnergieElastique((W_d)/V); // energ.ChangeEnergieElastique((W_d+W_v)/V); LibereTenseur(); }; // calcul des contraintes a t+dt et de ses variations void Maheo_hyper::Calcul_DsigmaHH_tdt (TenseurHH& ,TenseurBB& ,DdlElement & tab_ddl ,BaseB& ,TenseurBB & ,TenseurHH & ,BaseB& ,Tableau & ,BaseH& ,Tableau & ,TenseurBB & epsBB_tdt,Tableau & d_epsBB,TenseurBB & ,TenseurBB & gijBB_tdt,TenseurHH & gijHH_tdt_ ,Tableau & d_gijBB_tdt ,Tableau & ,double& jacobien_0,double& jacobien ,Vecteur& ,TenseurHH& sigHH_tdt,Tableau & d_sigHH ,EnergieMeca & energ,const EnergieMeca & ,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 << " Maheo_hyper::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 << " Maheo_hyper::Calcul_SDsigmaHH_tdt\n"; Sortie(1); }; #endif Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // passage en dim 3 explicite Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); // passage en dim 3 explicite // init coef matériaux Qsig_use = Qsig_rev; mu1_use = mu1_rev; mu2_use = mu2_rev; // cas de la thermo dépendance, on calcul les grandeurs if (Qsig_rev_temperature != NULL) Qsig_use = Qsig_rev_temperature->Valeur(*temperature); if (mu1_rev_temperature != NULL) mu1_use = mu1_rev_temperature->Valeur(*temperature); if (mu2_rev_temperature != NULL) mu2_use = mu2_rev_temperature->Valeur(*temperature); // récup du conteneur spécifique du point, pour sauvegarde éventuelle SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul); // dans le cas de l'utilisation de fct nD if ( (Qsig_nD != NULL) || (mu1_nD != NULL) || (mu2_nD != NULL) ) { // opération de transmission de la métrique const Met_abstraite::Impli* ex_impli = &ex; const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL; const Met_abstraite::Umat_cont* ex_expli = NULL; if (Qsig_nD != NULL) // là il faut calculer la fonction nD { // on utilise la méthode générique de loi abstraite list list_save; // inter pour l'appel de la fonction list_save.push_back(saveResul); Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee (Qsig_nD,1 // une seule valeur attendue en retour ,ex_impli,ex_expli_tdt,ex_expli ,NULL ,NULL ,&list_save ); // on récupère le premier élément du tableau uniquement Qsig_use *= tab_val(1); }; if (mu1_nD != NULL) // là il faut calculer la fonction nD { // on utilise la méthode générique de loi abstraite list list_save; // inter pour l'appel de la fonction list_save.push_back(saveResul); Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee (mu1_nD,1 // une seule valeur attendue en retour ,ex_impli,ex_expli_tdt,ex_expli ,NULL ,NULL ,&list_save ); // on récupère le premier élément du tableau uniquement mu1_use *= tab_val(1); }; if (mu2_nD != NULL) // là il faut calculer la fonction nD { // on utilise la méthode générique de loi abstraite list list_save; // inter pour l'appel de la fonction list_save.push_back(saveResul); Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee (mu2_nD,1 // une seule valeur attendue en retour ,ex_impli,ex_expli_tdt,ex_expli ,NULL ,NULL ,&list_save ); // on récupère le premier élément du tableau uniquement mu2_use *= tab_val(1); }; }; if (sortie_post) {// sauvegarde des paramètres matériau (*save_resulHyper_W.para_loi)(1) = Qsig_use; (*save_resulHyper_W.para_loi)(2) = mu1_use; (*save_resulHyper_W.para_loi)(3) = mu2_use; }; // calcul des invariants et de leurs variations premières et secondes Invariants_et_var2(*(ex.gijBB_0),*(ex.gijHH_0),gijBB_tdt,gijHH_tdt,jacobien_0,jacobien); // calcul du potentiel et de ses dérivées premières et secondes / aux invariants J_r Potentiel_et_var2(); // stockage éventuel pour du post-traitement if (sortie_post) { // récup du conteneur spécifique du point, pour sauvegarde éventuelle SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul); save_resulHyper_W.invP->potentiel= W_d + W_v; }; // calcul du tenseur des contraintes: on tiend compte du fait que le potentiel ne dépend que de J1 double unSurV=1./V; sigHH = (W_d_J1*unSurV) * d_J_r_epsBB_HH(1); // sigHH = (W_d_J1*unSurV) * d_J_r_epsBB_HH(1) + (W_d_J2*unSurV) * d_J_r_epsBB_HH(2) // + (W_v_J3*unSurV) * d_J_r_epsBB_HH(3); // cout << "\n sigHH " << sigHH // << "\n d_J_r_epsBB_HH1 " << d_J_r_epsBB_HH(1) // << "\n d_J_r_epsBB_HH2 " << d_J_r_epsBB_HH(2) // << "\n d_J_r_epsBB_HH3 " << d_J_r_epsBB_HH(3) // << "\n (W_d_J1*unSurV) " << (W_d_J1*unSurV) // << " (W_d_J2*unSurV) " << (W_d_J2*unSurV) // << " (W_d_J3*unSurV) " << (W_v_J3*unSurV); // calcul de la variation seconde du potentiel par rapport à epsij epskl // on tiend compte du fait que le potentiel ne dépend que de J1 Tenseur3HHHH d2W_d2epsHHHH = // tout d'abord les dérivées secondes du potentiel déviatoire W_d_J1_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1)) // puis les dérivées premières du potentiel déviatoire + W_d_J1 * d_J_1_eps2BB_HHHH; // = // tout d'abord les dérivées secondes du potentiel déviatoire // W_d_J1_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1)) // + W_d_J1_J2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(2)) // + W_d_J2_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(2),d_J_r_epsBB_HH(2)) // // puis les dérivées premières du potentiel déviatoire // + W_d_J1 * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH // // enfin les dérivées seconde et première du potentiel sphérique // + W_v_J3 * d_J_3_eps2BB_HHHH // + W_v_J3J3 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3)); // calcul de la variation du tenseur des contraintes par rapports aux déformations // on tient compte du fait que V*sigHH = d W/ d epsij Tenseur3HH interHH = -sigHH ; //* (-0.5*unSurV*unSurV); // Tenseur3HHHH dSigdepsHHHH(1,interHH,d_J_r_epsBB_HH(3)); Tenseur3HHHH dSigdepsHHHH(1,interHH,gijHH_tdt); // dSigdepsHHHH += (unSurV) * d2W_d2epsHHHH; Tenseur3HHHH interHHHH((unSurV) * d2W_d2epsHHHH); // cas des tenseurs généraux dSigdepsHHHH += interHHHH; // cas des tenseurs généraux //--------------------------------------------------------------------------------- // vérif numérique de l'opérateur tangent // Cal_dsigma_deps_num (*(ex.gijBB_0),*(ex.gijHH_0),gijBB_tdt,gijHH_tdt,jacobien_0,jacobien,dSigdepsHHHH); //--------------------------------------------------------------------------------- // calcul des variations / aux ddl int nbddl = d_gijBB_tdt.Taille(); for (int i = 1; i<= nbddl; i++) { // on fait uniquement une égalité d'adresse et de ne pas utiliser // le constructeur d'ou la profusion d'* et de () Tenseur3HH & dsigHH = *((Tenseur3HH*) (d_sigHH(i))); // passage en dim 3 const Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); // " dsigHH = dSigdepsHHHH && depsBB; }; // calcul du module de compressibilité module_compressibilite = 2. * W_v_J3;// ici 0 // pour le module de cisaillement, pour l'instant je ne fais rien !! à voir *** module_cisaillement = 0.; // traitement des énergies energ.Inita(0.); energ.ChangeEnergieElastique((W_d+W_v)/V); LibereTenseur(); LibereTenseurQ(); }; // 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 Maheo_hyper::Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & ,TenseurBB& ,TenseurBB & epsBB_tdt,TenseurBB &, double& jacobien_0,double& jacobien ,TenseurHH& sigHH_tdt,TenseurHHHH& d_sigma_deps_ ,EnergieMeca & energ,const EnergieMeca & ,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 << " Maheo_hyper::Calcul_dsigma_deps\n"; Sortie(1); }; #endif Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // // passage en dim 3 explicite Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) ex.gijHH_tdt); // passage en dim 3 explicite // init coef matériaux Qsig_use = Qsig_rev; mu1_use = mu1_rev; mu2_use = mu2_rev; // cas de la thermo dépendance, on calcul les grandeurs if (Qsig_rev_temperature != NULL) Qsig_use = Qsig_rev_temperature->Valeur(*temperature); if (mu1_rev_temperature != NULL) mu1_use = mu1_rev_temperature->Valeur(*temperature); if (mu2_rev_temperature != NULL) mu2_use = mu2_rev_temperature->Valeur(*temperature); // récup du conteneur spécifique du point, pour sauvegarde éventuelle SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul); // dans le cas de l'utilisation de fct nD if ( (Qsig_nD != NULL) || (mu1_nD != NULL) || (mu2_nD != NULL) ) { // opération de transmission de la métrique const Met_abstraite::Impli* ex_impli = NULL; const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL; const Met_abstraite::Umat_cont* ex_expli = &ex; if (Qsig_nD != NULL) // là il faut calculer la fonction nD { // on utilise la méthode générique de loi abstraite list list_save; // inter pour l'appel de la fonction list_save.push_back(saveResul); Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee (Qsig_nD,1 // une seule valeur attendue en retour ,ex_impli,ex_expli_tdt,ex_expli ,NULL ,NULL ,&list_save ); // on récupère le premier élément du tableau uniquement Qsig_use *= tab_val(1); }; if (mu1_nD != NULL) // là il faut calculer la fonction nD { // on utilise la méthode générique de loi abstraite list list_save; // inter pour l'appel de la fonction list_save.push_back(saveResul); Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee (mu1_nD,1 // une seule valeur attendue en retour ,ex_impli,ex_expli_tdt,ex_expli ,NULL ,NULL ,&list_save ); // on récupère le premier élément du tableau uniquement mu1_use *= tab_val(1); }; if (mu2_nD != NULL) // là il faut calculer la fonction nD { // on utilise la méthode générique de loi abstraite list list_save; // inter pour l'appel de la fonction list_save.push_back(saveResul); Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee (mu2_nD,1 // une seule valeur attendue en retour ,ex_impli,ex_expli_tdt,ex_expli ,NULL ,NULL ,&list_save ); // on récupère le premier élément du tableau uniquement mu2_use *= tab_val(1); }; }; if (sortie_post) {// sauvegarde des paramètres matériau (*save_resulHyper_W.para_loi)(1) = Qsig_use; (*save_resulHyper_W.para_loi)(2) = mu1_use; (*save_resulHyper_W.para_loi)(3) = mu2_use; }; // calcul des invariants et de leurs variations premières et secondes Invariants_et_var2(*(ex.gijBB_0),*(ex.gijHH_0),*(ex.gijBB_tdt),gijHH_tdt,jacobien_0,jacobien); // calcul du potentiel et de ses dérivées premières et secondes / aux invariants J_r Potentiel_et_var2(); // stockage éventuel pour du post-traitement if (sortie_post) { // récup du conteneur spécifique du point, pour sauvegarde éventuelle SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul); save_resulHyper_W.invP->potentiel= W_d + W_v; }; // calcul du tenseur des contraintes, on travaille ici dans le repère matériel finale correspondant // aux coordonnées initiales X_(0)^a, on obtient donc un tenseur dans la base naturelle finale double unSurV=1./V; // on tiend compte du fait que le potentiel ne dépend que de J1 Tenseur3HH sig_localeHH = (W_d_J1*unSurV) * d_J_r_epsBB_HH(1); // Tenseur3HH sig_localeHH = (W_d_J1*unSurV) * d_J_r_epsBB_HH(1) + (W_d_J2*unSurV) * d_J_r_epsBB_HH(2) // + (W_v_J3*unSurV) * d_J_r_epsBB_HH(3); // passage éventuelle dans la base I_a if (en_base_orthonormee) {sig_localeHH.BaseAbsolue(sigHH,*(ex.giB_tdt));} else {sigHH = sig_localeHH;}; // sinon la base locale est la bonne // calcul de la variation seconde du potentiel par rapport à epsij epskl // on tiend compte du fait que le potentiel ne dépend que de J1 Tenseur3HHHH d2W_d2epsHHHH = // tout d'abord les dérivées secondes du potentiel déviatoire W_d_J1_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1)) // puis les dérivées premières du potentiel déviatoire + W_d_J1 * d_J_1_eps2BB_HHHH; // = // tout d'abord les dérivées secondes du potentiel déviatoire // W_d_J1_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1)) // + W_d_J1_J2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(2)) // + W_d_J2_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(2),d_J_r_epsBB_HH(2)) // // puis les dérivées premières du potentiel déviatoire // + W_d_J1 * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH // // enfin les dérivées seconde et première du potentiel sphérique // + W_v_J3 * d_J_3_eps2BB_HHHH // + W_v_J3J3 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3)); // calcul de la variation du tenseur des contraintes par rapports aux déformations // on tient compte du fait que V*sigHH = d W/ d epsij Tenseur3HH interHH = -sig_localeHH ; //* (-0.5*unSurV*unSurV); // calcul de la variation du tenseur des contraintes par rapports aux déformations // on tient compte du fait que V*sigHH = d W/ d epsij Tenseur3HHHH dSigdepsHHHH(1,interHH,gijHH_tdt); // dSigdepsHHHH += (unSurV) * d2W_d2epsHHHH; Tenseur3HHHH interHHHH((unSurV) * d2W_d2epsHHHH); // cas des tenseurs généraux dSigdepsHHHH += interHHHH; // cas des tenseurs généraux // transfert des informations: on pas d'un tenseur de 81 composantes à 36 // avec des symétries par rapport aux deux premiers indices et par rapport aux deux derniers /// Tenseur3HHHH d_sigma_depsHHHH; d_sigma_depsHHHH.TransfertDunTenseurGeneral(dSigdepsHHHH.Symetrise1et2_3et4()); // calcul de la première partie de l'opérateur tangent (correspond au changement de repère // gi_tdt -> Ia de l'opérateur calculer précédemment Tenseur3HHHH & d_sigma_depsFinHHHH = *((Tenseur3HHHH*) &d_sigma_deps_); // pour accés directe // passage éventuelle dans la base I_a if (en_base_orthonormee) {dSigdepsHHHH.ChangeBase(d_sigma_depsFinHHHH,*(ex.giB_tdt));} else {d_sigma_depsFinHHHH = dSigdepsHHHH;}; // traitement des énergies energ.Inita(0.); energ.ChangeEnergieElastique((W_d)/V); // energ.ChangeEnergieElastique((W_d+W_v)/V); LibereTenseur(); LibereTenseurQ(); }; //---------------------- méthodes privées ------------------------------- // calcul du potentiel et de ses dérivées premières / aux invariants J_r void Maheo_hyper::Potentiel_et_var() { // calcul de grandeurs intermédiaires #ifdef MISE_AU_POINT if ((J_r(1)-3.) < 0.) { cout << "\nErreur : (J_r(1)-3.)= "<<(J_r(1)-3.) <<" est < 0. !" << " ce n'est pas possible avec les formules actuelles "; cout << "\n Maheo_hyper::Calcul_dsigma_deps\n"; Sortie(1); }; if ((mu1_use/Qsig_use) < 0.) { cout << "\nErreur : mu1_rev/Qsig_rev= "<Valeur(*temperature); if (mu1_rev_temperature != NULL) mu1_use = mu1_rev_temperature->Valeur(*temperature); if (mu2_rev_temperature != NULL) mu2_use = mu2_rev_temperature->Valeur(*temperature); // cas des contraintes et de ses variations analytiques // Tenseur3HHHH dSigdepsHHHH; // le tenseur contenant les dérivées analytiques Tenseur3HH SigmaHH_deb; Cal_sigmaEtDer_pour_num(gijBB_0_,gijHH_0_,gijBB_tdt_,gijHH_tdt_ ,jacobien_0,jacobien,SigmaHH_deb,dSigdepsHHHH); // dimensionnement pour la matrice numérique Tenseur3HHHH dSigdepsHHHH_num; // on va boucler sur les composantes de gijBB for (int i=1;i<=3;i++) for (int j=1;j<=3;j++) { gijBBtdt_N = gijBB_tdt; gijBBtdt_N.Coor(i,j) += delta; // en fait dans l'opération précédente on a modifier les termes (i,j) et (j,i) // car le tenseur est symétrique // on a donc en variation numérique la somme des deux dérivées // on définit un coeff multiplicatif qui vaut 1 ou 0.5 double coef=1.; if (i != j) coef = 0.5; gijHHtdt_N = gijBBtdt_N.Inverse(); double jacobien_N=sqrt(gijBBtdt_N.Det()); // cas des contraintes Tenseur3HH SigmaHH_N; Cal_sigma_pour_num(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N ,jacobien_0,jacobien_N,SigmaHH_N); // calcul des dérivées numériques et comparaisons for (int k=1;k<=3;k++) for (int l=1;l<=3;l++) { // double derSigNum = coef * 2.*(SigmaHH_N(k,l) - SigmaHH_deb(k,l) )*unSurDelta; dSigdepsHHHH_num.Change(k,l,i,j,derSigNum); double derSigAna = dSigdepsHHHH(k,l,i,j);//0.5*(dSigdepsHHHH(k,l,i,j) + dSigdepsHHHH(k,l,j,i)); bool erreur = false; if (diffpourcent(derSigNum,derSigAna,MaX(Dabs(derSigNum),Dabs(derSigAna)),0.1)) if (MaX(Dabs(derSigNum),Dabs(derSigAna)) > 200.) {if (MiN(Dabs(derSigNum),Dabs(derSigAna)) == 0.) {if ( MaX(Dabs(derSigNum),Dabs(derSigAna)) > 50.*delta) erreur = true;} else erreur = true; }; // erreur = false; // a virer if (erreur) { // calcul des dérivées d'éléments intermédiaires pour voir // cout << "\n erreur dans le calcul analytique de l'operateur tangent " << "\n derSigNum= " << derSigNum << " derSigAna= " << derSigAna << " klij= " << k << " " << l << " " << i << " " << j << " SigmaHH_N(k,l)= " << SigmaHH_N(k,l); cout << "\n Maheo_hyper::Calcul_derivee_numerique(.."; cout << "\n un caractere "; // --- pour le débug ---- // calcul des invariants et de leurs variations premières en numérique Invariants_et_var1_deb(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N ,jacobien_0,jacobien_N); // calcul des invariants et de leurs variations premières et secondes Invariants_et_var2_deb(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N ,jacobien_0,jacobien_N); string toto; toto= lect_chaine(); }; }; }; // passage des dérivées numériques aux dérivées finales dSigdepsHHHH= dSigdepsHHHH_num; }; // calcul de la contrainte avec le minimum de variable de passage, utilisé pour le numérique void Maheo_hyper::Cal_sigma_pour_num(const TenseurBB & gijBB_0,const TenseurHH & gijHH_0 ,const TenseurBB & gijBB_tdt,const TenseurHH & gijHH_tdt ,const double& jacobien_0,const double& jacobien,TenseurHH & sigHH_) { Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // passage en dim 3 explicite // calcul des invariants et de leurs variations premières (méthode de Hyper_W_gene_3D) // Invariants_et_var1(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien); // pour vérif on appelle var2, mais c'est à virer Invariants_et_var1(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien); // calcul du potentiel et de ses dérivées premières / aux invariants J_r Potentiel_et_var(); // calcul du tenseur des contraintes sigHH = (W_d_J1/V) * d_J_r_epsBB_HH(1); // sigHH = (W_d_J1/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2) // + (W_v_J3/V) * d_J_r_epsBB_HH(3); }; // idem avec la variation void Maheo_hyper::Cal_sigmaEtDer_pour_num(const TenseurBB & gijBB_0,const TenseurHH & gijHH_0 ,const TenseurBB & gijBB_tdt,const TenseurHH & gijHH_tdt ,const double& jacobien_0,const double& jacobien ,TenseurHH & sigHH_,Tenseur3HHHH& dSigdepsHHHH) { Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // passage en dim 3 explicite // calcul des invariants et de leurs variations premières et seconde (méthode de Hyper_W_gene_3D) Invariants_et_var2(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien); // calcul du potentiel et de ses dérivées premières / aux invariants J_r Potentiel_et_var2(); // calcul du tenseur des contraintes double unSurV=1./V; sigHH = (W_d_J1*unSurV) * d_J_r_epsBB_HH(1); // sigHH = (W_d_J1*unSurV) * d_J_r_epsBB_HH(1) + (W_d_J2*unSurV) * d_J_r_epsBB_HH(2) // + (W_v_J3*unSurV) * d_J_r_epsBB_HH(3); // calcul de la variation seconde du potentiel par rapport à epsij epskl Tenseur3HHHH d2W_d2epsHHHH = W_d_J1 * d_J_1_eps2BB_HHHH; // = W_d_J1 * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH // + W_v_J3 * d_J_3_eps2BB_HHHH // + W_v_J3J3 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3)); // cout << "\n d2W_d2epsHHHH(1,1,1,1)= " << d2W_d2epsHHHH(1,1,1,1); // calcul de la variation du tenseur des contraintes par rapports aux déformations // on tient compte du fait que V*sigHH = d W/ d epsij Tenseur3HH interHH = -sigHH ; //* (-0.5*unSurV*unSurV); Tenseur3HHHH d_igdepsHHHH(1,interHH,gijHH_tdt); d_igdepsHHHH += (unSurV) * d2W_d2epsHHHH; dSigdepsHHHH = d_igdepsHHHH; };