// 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 "ElemMeca.h" #include #include "ConstMath.h" #include "Util.h" #include "Coordonnee2.h" #include "Coordonnee3.h" #include "TypeQuelconqueParticulier.h" #include "TypeConsTens.h" // affichage dans la sortie transmise, des variables duales "nom" // aux differents points d'integration // dans le cas ou nom est vide, affichage de "toute" les variables void ElemMeca::VarDualSort(ofstream& sort, Tableau& nom,int ,int cas) { if ((cas == 1) || (cas == 2)) { // cas d'une premiere initialisation Tableau tab(5); tab(1) = iM0; tab(2) = iMt; tab(3) = iMtdt; tab(4) = igiH_0; tab(5) = igijHH_0; met->PlusInitVariables(tab) ; }; // ----- def de grandeurs de travail // il y a deux pb a gérer: 1) le fait que la dimension absolue peut-être différente de la dimension des tenseurs // 2) le fait que l'on veut une sortie dans une base ad hoc ou pas int dim = lesPtIntegMecaInterne->DimTens(); int dim_sortie_tenseur = dim; // dans le cas ou l'on veut une sortie en base absolue, il faut que dim_sortie_tenseur = la dimension de la base absolue // if (absolue) // dim_sortie_tenseur = ParaGlob::Dimension(); // pour ne faire qu'un seul test ensuite // bool prevoir_change_dim_tenseur = false; // if ((absolue) && (dim != dim_sortie_tenseur)) // prevoir_change_dim_tenseur = true; // a priori ici on sort dans un repère ad hoc bool absolue=false; // def de la dimension des tenseurs // int dim = lesPtIntegMecaInterne->DimTens(); // def de tenseurs pour la sortie TenseurHH& sig0HH = *(NevezTenseurHH(dim_sortie_tenseur)) ; TenseurHH& sigHH = *(NevezTenseurHH(dim_sortie_tenseur)) ; TenseurBB& eps0BB = *(NevezTenseurBB(dim_sortie_tenseur)) ; // pour def green-lagrange TenseurBB& epsBB = *(NevezTenseurBB(dim_sortie_tenseur)) ; // pour def courante TenseurBB& epslogBB = *(NevezTenseurBB(dim_sortie_tenseur)) ; // pour def logarithmique TenseurBB& epsAlmBB = *(NevezTenseurBB(dim_sortie_tenseur)) ; // pour def d'almansi TenseurHB& sigHB = *(NevezTenseurHB(dim_sortie_tenseur)) ; TenseurHB& epsHB = *(NevezTenseurHB(dim_sortie_tenseur)) ; TenseurBB& DepsBB = *(NevezTenseurBB(dim_sortie_tenseur)) ; TenseurHB& DepsHB = *(NevezTenseurHB(dim_sortie_tenseur)) ; TenseurBB& DeltaEpsBB = *(NevezTenseurBB(dim_sortie_tenseur)) ; Enum_dure temps; int ni; // compteur globale de point d'integration for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++) // for (int ni=1;ni<= nbi;ni++) {// éléments de métrique et matrices de passage TenseurHH* gijHH;TenseurBB* gijBB; Coordonnee* Mpt0;Coordonnee* Mptfin; Mat_pleine Aa0(dim,dim),Aafin(dim,dim); if ((cas==1) || (cas==11)) // calcul d'une ortho interessance de visualisation des tenseurs // cas de tenseur 3D -> Ipa, cas 1D on prend un vecteur norme collineaire // avec g1, dans le cas 2D // la nouvelle base gamma est calculee dans def par projection de "Ipa" sur Galpha // le resultat est une matrice de passage utilisable pour le changement de base // Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a // tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer // resultat a t+dt {const Met_abstraite::InfoImp& ex = def->RemontImp(absolue,Aa0,Aafin); temps=TEMPS_tdt; Mpt0 = ex.M0; Mptfin = ex.Mtdt; gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; } else if ((cas==2) || (cas==12)) { // resultat a tdt const Met_abstraite::InfoExp_tdt& ex= def->RemontExp_tdt(absolue,Aa0,Aafin); temps=TEMPS_tdt; Mpt0 = ex.M0; Mptfin = ex.Mtdt; gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; }; // pour les formules de passage de repère il nous faut : // Ip_a = beta_a^{.j} g_j et Ip^b = gamma^b_{.j} g^j // on a: [beta_a^{.j}] = [Aa^j_{.a}]^T // et [gamma^b_{.j}] = [beta_a^{.j}]^{-1T} = [Aa^j_{.a}]^{-1} Mat_pleine gamma0(Aa0.Inverse()); Mat_pleine gammafin(Aafin.Inverse()); // on détermine également les matrices beta Mat_pleine beta0(Aa0.Transpose()); Mat_pleine betafin(Aafin.Transpose()); // calcul des tenseurs PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni); sig0HH = *(ptIntegMeca.SigHH()); sig0HH.ChBase(gamma0); // changement de base initiale sigHH = *(ptIntegMeca.SigHH()); sigHH.ChBase(gammafin); // changement de base finale eps0BB = *(ptIntegMeca.EpsBB()); eps0BB.ChBase(beta0); // changement de base initiale epsBB = *(ptIntegMeca.EpsBB());; epsBB.ChBase(betafin); // changement de base finale sigHB = *(ptIntegMeca.SigHH()) * (*gijBB); epsHB = (*gijHH) * (*(ptIntegMeca.EpsBB())); DepsBB = *(ptIntegMeca.DepsBB()); DepsBB.ChBase(betafin); // changement de base finale DepsHB = (*gijHH) * (*(ptIntegMeca.DepsBB())); DeltaEpsBB = *(ptIntegMeca.DeltaEpsBB()); DeltaEpsBB.ChBase(betafin); // changement de base finale switch (def->Type_de_deformation()) { case DEFORMATION_STANDART : // c'est à dire almansi { eps0BB = *(ptIntegMeca.EpsBB()); eps0BB.ChBase(beta0); // changement de base initiale epsAlmBB = epsBB; // si l'on veut sortir la déformation logarithmique le plus simple est de la calculer def->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE); def->Cal_deformation (temps,epslogBB); def->Change_type_de_deformation(DEFORMATION_STANDART); // on revient au type initial epslogBB.ChBase(betafin); break; } case DEFORMATION_LOGARITHMIQUE : { epslogBB=epsBB; // si l'on veut sortir la déformation d'Almansi ou de green-lagrange le plus simple est de les calculer def->Change_type_de_deformation(DEFORMATION_STANDART); def->Cal_deformation (temps,epsAlmBB); def->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE); // on revient au type initial epsAlmBB.ChBase(betafin); eps0BB=epsAlmBB; eps0BB.ChBase(beta0); break; } default: cout << "\n cas de deformation non encore implante en sortie de visualisation " << Nom_type_deformation(def->Type_de_deformation()) << " affichage donc errone des valeurs !!!"; }; int cas=0; Coordonnee valPropreSig = sigHB.ValPropre(cas); if (cas == -1) { cout << "\n warning *** erreur dans le calcul des valeurs propres de la contrainte"; if (ParaGlob::NiveauImpression() >= 7) {sigHB.Ecriture(cout); cout << "\nElemMeca::VarDualSort(...";}; cout << endl;}; Coordonnee valPropreEps = epsHB.ValPropre(cas); if (cas == -1) { cout << "\n warning *** erreur dans le calcul des valeurs propres de la deformation"; if (ParaGlob::NiveauImpression() >= 7) {epsHB.Ecriture(cout); cout << "\nElemMeca::VarDualSort(...";}; cout << endl;}; Coordonnee valPropreDeps = DepsHB.ValPropre(cas); if (cas == -1) { cout << "\n warning *** erreur dans le calcul des valeurs propres de la vitesse de deformation"; if (ParaGlob::NiveauImpression() >= 7) {DepsHB.Ecriture(cout); cout << "\nElemMeca::VarDualSort(...";}; cout << endl;}; // def de la contrainte de Mises double Mises = 0.; Coordonnee& vv = valPropreSig; // pour condenser l'écriture if (dim ==1) {Mises = Dabs(vv(1));} else if (dim ==2) {Mises = sqrt( ((vv(1)-vv(2))*(vv(1)-vv(2)) + vv(1) * vv(1) + vv(2) * vv(2)) * 0.5); } else {Mises = sqrt( ((vv(1)-vv(2))*(vv(1)-vv(2)) + (vv(1)-vv(3))*(vv(1)-vv(3)) + (vv(3)-vv(2))*(vv(3)-vv(2))) * 0.5);}; // def de la déformation duale de mises = sqrt(2/3 * eps_barre:eps_barre) //--- on utilise directement la grandeur stockée au pt d'integ, mais ici ne sert à rien, puis cette grandeur n'est pas utilisée par la suite ! // l'expression est la même que celle de mises, ormis un coeff 4/9 qui permet de passer de 3/2 à 2/3 // double defDualMises = 0.; Coordonnee& vdef = valPropreEps; // pour condenser l'écriture // if (dim ==1) {defDualMises = Dabs(vdef(1));} // else if (dim ==2) {defDualMises = sqrt( ((vdef(1)-vdef(2))*(vdef(1)-vdef(2)) // + vdef(1) * vdef(1) + vdef(2) * vdef(2)) * 2./3.); } // else {defDualMises = sqrt( ((vdef(1)-vdef(2))*(vdef(1)-vdef(2)) + (vdef(1)-vdef(3))*(vdef(1)-vdef(3)) // + (vdef(3)-vdef(2))*(vdef(3)-vdef(2))) * 2./9.);}; // donnees propre a la loi mécanique au pt d'integ Loi_comp_abstraite::SaveResul * sDon = tabSaveDon(ni); // donnees propre a la loi thermo physique au pt d'integ CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques if (loiTP != NULL) {sTP = tabSaveTP(ni);}; // au pt d'integ si elles existes // données propre à la déformation au pt d'integ Deformation::SaveDefResul * sDefDon = tabSaveDefDon(ni); // affichage sort << "\n=================== Point d\'integration [[[ " << ni<<" ]]] ==================" ; sort << "\nCoordonnees avant deformation " ; Mpt0->Affiche(sort,16); sort << "\nCoordonnees apres deformation " ; Mptfin->Affiche(sort,16); sort <<" \n et les variations " ; ((*Mptfin) - (*Mpt0)).Affiche(sort,16); if (nom.Taille() == 0) { // cas d'une sortie complete for (int indic =1; indic<=11; indic++) AffDefCont (sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,indic) ; } else for (int ij = 1;ij<= nom.Taille(); ij++) if (nom(ij) == "Green-Lagrange") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,1); else if (nom(ij) == "Almansi") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,2); else if (nom(ij) == "Cauchy_global") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,3); else if (nom(ij) == "Def_mixte_local") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,4); else if (nom(ij) == "Sigma_mixte_local") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,5); else if (nom(ij) == "Sigma_principale") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,6); else if (nom(ij) == "Def_principale") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,7); else if (nom(ij) == "Resultat_loi") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,8); else if (nom(ij) == "Vitesse_deformation") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,9); else if (nom(ij) == "Increment_deformation") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,10); else if (nom(ij) == "VitDef_principale") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,11); else if (nom(ij) == "logarithmique") AffDefCont(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB,valPropreEps,valPropreDeps,valPropreSig,Mises ,epsAlmBB,epslogBB,12); } // delete des tenseurs TenseurHH* ptHH; TenseurBB* ptBB; ptHH = &sig0HH; delete ptHH; ptHH = &sigHH; delete ptHH; ptBB = &eps0BB; delete ptBB; ptBB = &epsBB; delete ptBB; TenseurHB * ptHB = &sigHB; delete ptHB; ptHB = &epsHB; delete ptHB; // liberation des tenseurs intermediaires LibereTenseur(); }; void ElemMeca::AffDefCont(ofstream& sort,Loi_comp_abstraite::SaveResul * sDon, TenseurBB& eps0BB,TenseurBB& epsBB,TenseurBB& DepsBB,TenseurBB& DeltaEpsBB, TenseurHH& sigHH, TenseurHB& epsHB,TenseurHB& sigHB, Coordonnee& valPropreEps,Coordonnee& valPropreDeps,Coordonnee& valPropreSig, double Mises,TenseurBB& epsAlmBB,TenseurBB& epslogBB,int indic) { int dim = sigHH.Dimension(); switch (dim) { case 1: AffDefCont1D(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB, valPropreEps,valPropreDeps,valPropreSig,Mises,epsAlmBB,epslogBB,indic); break; case 2: AffDefCont2D(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB, valPropreEps,valPropreDeps,valPropreSig,Mises,epsAlmBB,epslogBB,indic); break; case 3: AffDefCont3D(sort,sDon,eps0BB,epsBB,DepsBB,DeltaEpsBB,sigHH,epsHB,sigHB, valPropreEps,valPropreDeps,valPropreSig,Mises,epsAlmBB,epslogBB,indic); break; default : {cout << "\nerreur seules les cas 1D 2D 3D sont pris considere et non " << dim; cout << "\nElemMeca::AffDefCont( etc ..." << endl; } } }; void ElemMeca::AffDefCont1D(ofstream& sort,Loi_comp_abstraite::SaveResul * saveDon, TenseurBB& eps0BB,TenseurBB& ,TenseurBB& DepsBB,TenseurBB& DeltaEpsBB, TenseurHH& sigHH, TenseurHB& epsHB,TenseurHB& sigHB, Coordonnee& valPropreEps,Coordonnee& valPropreDeps,Coordonnee& valPropreSig, double Mises,TenseurBB& epsAlmBB,TenseurBB& epslogBB,int indic) { switch (indic) {case 1 :{ sort << "\nDeformations Green-Lagrange dans repere global"; sort << "\n E11 " ; sort << setw (16) << eps0BB(1,1) ; } break; case 2 :{ sort << "\nDeformations Almansi dans repere global"; sort << "\n eps11 " ; sort << setw (16) << epsAlmBB(1,1) ; } break; case 3 :{ sort << "\nContraintes Cauchy dans repere global"; sort << "\n Sig11 " ; sort << setw (16) << sigHH(1,1) ; } break; case 4 :{ sort << "\nDeformations mixtes dans repere local"; sort << "\n epsHB11 " ; sort << setw (16) << epsHB(1,1) ; } break; case 5 :{ sort << "\nContraintes mixtes dans repere local"; sort << "\n sigHB11 " ; sort << setw (16) << sigHB(1,1) ; } break; case 6 :{ sort << "\nContraintes principales"; sort << "\n SigI " << " VonMises " ; sort << setw (16) << valPropreSig(1) << " "<< setw (16) << Mises; } break; case 7 :{ sort << "\nDeformations principales"; sort << "\n EpsI " ; sort << setw (16) << valPropreEps(1) ; } break; case 8 :{ if (saveDon != NULL) { sort << "\ngrandeurs attachees a la loi et au point "; loiComp->AfficheDataSpecif(sort,saveDon); } } break; case 9 :{ sort << "\nVitesse de Deformations dans le repere global"; sort << "\n Deps11 " ; sort << setw (16) << DepsBB(1,1) ; } break; case 10 :{ sort << "\nIncrement de Deformations d'Almansi dans repere global"; sort << "\n DeltaEps11 " ; sort << setw (16) << DeltaEpsBB(1,1) ; } break; case 11 :{ sort << "\nVitesses de Deformations principales"; sort << "\n DepsI " ; sort << setw (16) << valPropreDeps(1) ; } break; case 12 :{ sort << "\nDeformations logarithmique dans repere global"; sort << "\n epslog11 " ; sort << setw (16) << epslogBB(1,1) ; } break; default : {cout << "\nerreur seules les cas de 1 a 11 sont pris en consideration et non " << indic; cout << "\nElemMeca::AffDefCont1D( etc ..." << endl; } } }; void ElemMeca::AffDefCont2D(ofstream& sort,Loi_comp_abstraite::SaveResul * saveDon, TenseurBB& eps0BB,TenseurBB& ,TenseurBB& DepsBB,TenseurBB& DeltaEpsBB, TenseurHH& sigHH, TenseurHB& epsHB,TenseurHB& sigHB, Coordonnee& valPropreEps,Coordonnee& valPropreDeps,Coordonnee& valPropreSig, double Mises,TenseurBB& epsAlmBB,TenseurBB& epslogBB, int indic) { string type_repere=" repere global"; if (ParaGlob::Dimension() == 3) // lorsque que l'on sort des tenseurs d'ordre 2 en dimension 3, on utilise un repère particulier local orthonormé type_repere=" repere orthonorme local"; switch (indic) {case 1 :{ sort << "\nDeformations Green-Lagrange dans" << type_repere; sort << "\n E11 " << " E22 " << " E12 " << " E21 \n"; sort << setw (16) << eps0BB(1,1) << " "<< setw (16) << eps0BB(2,2) << " " << setw (16) << eps0BB(1,2) << " "<< setw (16) << eps0BB(2,1) << " "; } break; case 2 :{ sort << "\nDeformations Almansi dans" << type_repere; sort << "\n eps11 " << " eps22 " << " eps12 " << " eps21 \n"; sort << setw (16) << epsAlmBB(1,1) << " "<< setw (16) << epsAlmBB(2,2) << " " << setw (16) << epsAlmBB(1,2) << " "<< setw (16) << epsAlmBB(2,1) << " "; } break; case 12 :{ sort << "\nDeformations logarithmique dans" << type_repere; sort << "\n epslog11 " << " epslog22 " << " epslog12 " << " epslog21 \n"; sort << setw (16) << epslogBB(1,1) << " "<< setw (16) << epslogBB(2,2) << " " << setw (16) << epslogBB(1,2) << " "<< setw (16) << epslogBB(2,1) << " "; } break; case 3 :{ sort << "\nContraintes Cauchy dans" << type_repere; sort << "\n Sig11 " << " Sig22 " << " Sig12 " << " Sig21 \n"; sort << setw (16) << sigHH(1,1) << " "<< setw (16) << sigHH(2,2) << " " << setw (16) << sigHH(1,2) << " "<< setw (16) << sigHH(2,1) << " "; } break; case 4 :{ sort << "\nDeformations mixtes dans repere local"; sort << "\n epsHB11 " << " epsHB22 " << " epsHB12 " << " epsHB21 \n"; sort << setw (16) << epsHB(1,1) << " "<< setw (16) << epsHB(2,2) << " " << setw (16) << epsHB(1,2) << " "<< setw (16) << epsHB(2,1) << " "; } break; case 5 :{ sort << "\nContraintes mixtes dans repere local"; sort << "\n sigHB11 " << " sigHB22 " << " sigHB12 " << " sigHB21 \n"; sort << setw (16) << sigHB(1,1) << " "<< setw (16) << sigHB(2,2) << " " << setw (16) << sigHB(1,2) << " "<< setw (16) << sigHB(2,1) << " "; } break; case 6 :{ sort << "\nContraintes principales"; sort << "\n SigI " << " SigII " << " VonMises \n" ; sort << setw (16) << valPropreSig(1) << " "<< setw (16) << valPropreSig(2) << " " << setw (16) << Mises; } break; case 7 :{ sort << "\nDeformations principales"; sort << "\n EpsI " << " EpsII \n" ; sort << setw (16) << valPropreEps(1) << " "<< setw (16) << valPropreEps(2) ; } break; case 8 :{ if (saveDon != NULL) { sort << "\ngrandeurs attachees a la loi et au point "; loiComp->AfficheDataSpecif(sort,saveDon); } } break; case 9 :{ sort << "\nVitesse de Deformations dans le" << type_repere; sort << "\n Deps11 " << " Deps22 " << " Deps12 " << " Deps21 \n"; sort << setw (16) << DepsBB(1,1) << " "<< setw (16) << DepsBB(2,2) << " " << setw (16) << DepsBB(1,2) << " "<< setw (16) << DepsBB(2,1) << " "; } break; case 10 :{ sort << "\nIncrement de Deformations d'Almansi dans" << type_repere; sort << "\n DeltaEpsBB11 " << " DeltaEpsBB22 " << " DeltaEpsBB12 " << " DeltaEpsBB21 \n"; sort << setw (16) << DeltaEpsBB(1,1) << " "<< setw (16) << DeltaEpsBB(2,2) << " " << setw (16) << DeltaEpsBB(1,2) << " "<< setw (16) << DeltaEpsBB(2,1) << " "; } break; case 11 :{ sort << "\nVitesses de Deformations principales"; sort << "\n DepsI " << " DepsII \n" ; sort << setw (16) << valPropreDeps(1) << " "<< setw (16) << valPropreDeps(2) ; } break; default : {cout << "\nerreur seules les cas de 1 a 11 sont pris en consideration et non " << indic; cout << "\nElemMeca::AffDefCont2D( etc ..." << endl; } } }; void ElemMeca::AffDefCont3D(ofstream& sort,Loi_comp_abstraite::SaveResul * saveDon, TenseurBB& eps0BB,TenseurBB& ,TenseurBB& DepsBB,TenseurBB& DeltaEpsBB, TenseurHH& sigHH, TenseurHB& epsHB,TenseurHB& sigHB, Coordonnee& valPropreEps,Coordonnee& valPropreDeps,Coordonnee& valPropreSig, double Mises,TenseurBB& epsAlmBB,TenseurBB& epslogBB,int indic) { switch (indic) {case 1 : { sort << "\nDeformations Green-Lagrange dans repere global"; sort << "\n E11 " << " E12 " << " E13 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << eps0BB(1,i) << " "; sort << "\n E21 " << " E22 " << " E23 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << eps0BB(2,i) << " "; sort << "\n E31 " << " E32 " << " E33 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << eps0BB(3,i) << " "; } break; case 2 : { sort << "\nDeformations Almansi dans repere global"; sort << "\n eps11 " << " eps12 " << " eps13 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << epsAlmBB(1,i) << " "; sort << "\n eps21 " << " eps22 " << " eps23 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << epsAlmBB(2,i) << " "; sort << "\n eps31 " << " eps32 " << " eps33 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << epsAlmBB(3,i) << " "; } break; case 12 : { sort << "\nDeformations logarithmique dans repere global"; sort << "\n epslog11 " << " epslog12 " << " epslog13 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << epslogBB(1,i) << " "; sort << "\n epslog21 " << " epslog22 " << " epslog23 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << epslogBB(2,i) << " "; sort << "\n epslog31 " << " epslog32 " << " epslog33 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << epslogBB(3,i) << " "; } break; case 3 : { sort << "\nContraintes Cauchy dans repere global"; sort << "\n Sig11 " << " Sig12 " << " Sig13 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << sigHH(1,i) << " "; sort << "\n Sig21 " << " Sig22 " << " Sig23 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << sigHH(2,i) << " "; sort << "\n Sig31 " << " Sig32 " << " Sig33 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << sigHH(3,i) << " "; } break; case 4 : { sort << "\nDeformations mixtes dans repere locall"; sort << "\n epsHB11 " << " epsHB12 " << " epsHB13 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << epsHB(1,i) << " "; sort << "\n epsHB21 " << " epsHB22 " << " epsHB23 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << epsHB(2,i) << " "; sort << "\n epsHB31 " << " epsHB32 " << " epsHB33 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << epsHB(3,i) << " "; } break; case 5 : { sort << "\nContraintes mixtes dans repere local"; sort << "\n sigHB11 " << " sigHB12 " << " sigHB13 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << sigHB(1,i) << " "; sort << "\n sigHB21 " << " sigHB22 " << " sigHB23 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << sigHB(2,i) << " "; sort << "\n sigHB31 " << " sigHB32 " << " sigHB33 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << sigHB(3,i) << " "; } break; case 6 : {sort << "\nContraintes principales"; sort << "\n SigI " << " SigII " << " SigIII \n"; for (int i=1;i<=3;i++) sort << setw (16) << valPropreSig(i) << " "; sort << "\n VonMises = " << Mises ; } break; case 7 : {sort << "\nDeformations principales"; sort << "\n EpsI " << " EpsII " << " EpsIII \n"; for (int i=1;i<=3;i++) sort << setw (16) << valPropreEps(i) << " "; } break; case 8 :{ if (saveDon != NULL) { sort << "\ngrandeurs attachees a la loi et au point "; loiComp->AfficheDataSpecif(sort,saveDon); } } break; case 9 : { sort << "\nVitesse de Deformations dans le repere global"; sort << "\n Deps11 " << " Deps12 " << " Deps13 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << DepsBB(1,i) << " "; sort << "\n Deps21 " << " Deps22 " << " Deps23 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << DepsBB(2,i) << " "; sort << "\n Deps31 " << " Deps32 " << " Deps33 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << DepsBB(3,i) << " "; } break; case 10 : { sort << "\nIncrement de Deformations d'Almansi dans repere global"; sort << "\n DeltaEpsBB11 " << " DeltaEpsBB12 " << " DeltaEpsBB13 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << DeltaEpsBB(1,i) << " "; sort << "\n DeltaEpsBB21 " << " DeltaEpsBB22 " << " DeltaEpsBB23 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << DeltaEpsBB(2,i) << " "; sort << "\n DeltaEpsBB31 " << " DeltaEpsBB32 " << " DeltaEpsBB33 \n" ; for (int i=1;i<=3;i++) sort << setw (16) << DeltaEpsBB(3,i) << " "; } break; case 11 : { sort << "\nVitesses de Deformations principales"; sort << "\n DepsI " << " DepsII " << " DepsIII \n"; for (int i=1;i<=3;i++) sort << setw (16) << valPropreDeps(i) << " "; } break; default : {cout << "\nerreur seules les cas de 1 a 11 sont pris en consideration et non " << indic; cout << "\nElemMeca::AffDefCont2D( etc ..." << endl; } } }; // récupération des valeurs au numéro d'ordre = iteg pour // les grandeur enu // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière // cas = -1 ou -2: cela signifie que la métrique vient juste d'être calculée au pti iteg // sinon il faut recalculer qq éléments de la métrique // NB Importante: il faut faire attention à ce que ces métriques soient identiques à celles qui ont servit // pour le calcul des tenseurs: en particulier si c'est utilisé pour calculer les grandeurs pour le chargement // il faut s'assurer que ce sont les "mêmes pti" qui servent pour la charge et pour la raideur !! Tableau ElemMeca::Valeur_multi (bool absolue, Enum_dure temps,const List_io& enu,int iteg,int cas ) { // ----- def de grandeurs de travail // il y a deux pb a gérer: 1) le fait que la dimension absolue peut-être différente de la dimension des tenseurs // 2) le fait que l'on veut une sortie dans une base ad hoc ou pas int dim = lesPtIntegMecaInterne->DimTens();int dim_sortie_tenseur = dim; // dans le cas ou l'on veut une sortie en base absolue, il faut que dim_sortie_tenseur = la dimension de la base absolue if (absolue) dim_sortie_tenseur = ParaGlob::Dimension(); // --- pour ne faire qu'un seul test ensuite bool prevoir_change_dim_tenseur = false; // initialement on faisait le test suivant, // if ((absolue) && (dim != dim_sortie_tenseur)) if (absolue) // mais en fait, quand on sort en absolu on est obligé d'utiliser un tenseur intermédiaire // car BaseAbsolue(.. modifie tenseur passé en paramètre, donc dans tous les cas de sortie absolue // il faut un tenseur intermédiaire qui a ou non une dimension différente prevoir_change_dim_tenseur = true; #ifdef MISE_AU_POINT if(iteg > lesPtIntegMecaInterne->NbPti()) { cout << "\n erreur, les informations demandees au pti "< 3 ) cout << "\n ElemMeca::Valeur_multi(.."; cout << endl; Sortie (1); }; #endif PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(iteg); // recup de l'incrément de temps double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant(); double unSurDeltat=0; if (Abs(deltat) >= ConstMath::trespetit) {unSurDeltat = 1./deltat;} else // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand { // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb if (unSurDeltat < 0) { cout << "\n le pas de temps est négatif !! "; }; unSurDeltat = ConstMath::tresgrand; }; // -- def des tenseurs locaux TenseurHB* sigHB = NULL ; TenseurHB* sig_barreHB = NULL ; TenseurHB* epsHB = NULL ; TenseurHB* eps_barreHB = NULL ; TenseurBB* DepsBB = NULL ; TenseurHB* Deps_barreHB = NULL ; TenseurHB* DepsHB = NULL ; TenseurBB* epsAlmTotalBB=NULL; // pour la déformation totale d'almansi TenseurBB* epsGLTotalBB=NULL; // pour la déformation totale de green_lagrange TenseurBB* epsLogTotalBB=NULL; // pour la déformation totale logarithmique Coordonnee* Mtdt = NULL; // coordonnées finales éventuelles du point d'intégration considéré Coordonnee* Mt = NULL; // coordonnées à t éventuelles du point d'intégration considéré Coordonnee* M0 = NULL; // coordonnées initiales éventuelles du point d'intégration considéré Coordonnee* N_surf = NULL; // coordonnée d'un vecteur normal si c'est adéquate Coordonnee* N_surf_t = NULL; // coordonnée d'un vecteur normal à t si c'est adéquate Coordonnee* N_surf_t0 = NULL; // coordonnée d'un vecteur normal à t0 si c'est adéquate Coordonnee* Vitesse = NULL; // cas des vitesses // on définie des indicateurs pour ne pas faire plusieurs fois le même calcul List_io::const_iterator ie,iefin=enu.end(); bool besoin_des_contraintes=false; bool besoin_des_deformation=false; bool besoin_des_contraintes_barre=false; bool besoin_des_deformation_barre=false; bool besoin_des_vitesses_deformation=false; bool besoin_des_vitesse_deformation_barre=false; bool besoin_des_valpropre_sigma=false; bool besoin_des_valpropre_deformation = false; bool besoin_des_valpropre_vitdef = false; bool besoin_deformation_logarithmique = false; bool besoin_deformation_greenlagrange = false; bool besoin_deformation_almansi = false; bool besoin_coordonnees = false; bool besoin_deplacements = false; bool besoin_coordonnees_t = false;bool besoin_coordonnees_t0 = false; for (ie=enu.begin(); ie!=iefin;ie++) { if (Meme_famille((*ie).Enum(),SIG11)) besoin_des_contraintes=true; if (Meme_famille((*ie).Enum(),EPS11)) besoin_des_deformation=true; if (Meme_famille((*ie).Enum(),DEPS11)) besoin_des_vitesses_deformation=true; if (Meme_famille((*ie).Enum(),X1)) besoin_coordonnees=true; if (Meme_famille((*ie).Enum(),UX)) {besoin_deplacements=true;besoin_coordonnees=true;}; int posi = (*ie).Position()-NbEnum_ddl(); switch (posi) { case 1: case 2: case 3: case 4: case 5: case 6: {besoin_deformation_greenlagrange=true; break;} case 7: case 8: case 9: case 10: case 11: case 12: {besoin_deformation_almansi=true; break;} case 28: case 29: case 30: case 31: case 32: {besoin_des_valpropre_sigma=true; break;} case 25: case 26: case 27: case 77: {besoin_des_valpropre_deformation=true;break;} case 40: case 41: case 42: {besoin_des_valpropre_vitdef=true;break;} case 49: case 50: case 51: case 52: case 53: case 54: {besoin_deformation_logarithmique=true;break;} case 55: case 56: case 57: case 58: case 59: case 60: {if ((epsAlmTotalBB == NULL) && (dilatation)) {epsAlmTotalBB = (NevezTenseurBB(dim_sortie_tenseur));}; break; } case 61: case 62: case 63: case 64: case 65: case 66: { if ((epsGLTotalBB == NULL) && (dilatation)) {epsGLTotalBB = (NevezTenseurBB(dim_sortie_tenseur));}; break; } case 67: case 68: case 69: case 70: case 71: case 72: {if ((epsLogTotalBB == NULL) && (dilatation)) {epsLogTotalBB = (NevezTenseurBB(dim_sortie_tenseur));}; break; } case 78: case 79: case 80: {besoin_des_deformation_barre=true;break;} case 81: case 82: case 83: {besoin_des_contraintes_barre=true;break;} case 84: case 85: case 86: {besoin_des_vitesse_deformation_barre=true;break;} case 114: case 115: case 116: // le vecteur normal { N_surf = new Coordonnee(ParaGlob::Dimension()); break;} case 117: case 118: case 119: // le vecteur normal à t { N_surf_t = new Coordonnee(ParaGlob::Dimension()); break;} case 120: case 121: case 122: // le vecteur normal à t0 { N_surf_t0 = new Coordonnee(ParaGlob::Dimension()); break;} case 123: case 124: case 125: // la position à t { Mt = new Coordonnee(ParaGlob::Dimension()); besoin_coordonnees_t = true; break; } case 126: case 127: case 128: // la position à t0 { M0 = new Coordonnee(ParaGlob::Dimension()); besoin_coordonnees_t0 = true; break; } default: break; }; }; // -- def de tenseurs pour la sortie TenseurHH* sigHH = NULL ; TenseurBB* eps0BB = NULL ; // pour def green-lagrange TenseurBB* epsBB = NULL ; // pour def courante TenseurBB* epslogBB = NULL ; // pour def logarithmique TenseurBB* epsAlmBB = NULL ; // pour def d'almansi TenseurBB* DeltaEpsBB = NULL ; bool besoin_matrice_chg_base = false; if (besoin_des_contraintes || besoin_des_valpropre_sigma || besoin_des_contraintes_barre) {sigHH = (NevezTenseurHH(dim_sortie_tenseur)) ; sigHB = (NevezTenseurHB(dim)) ; sig_barreHB = (NevezTenseurHB(dim)) ;besoin_matrice_chg_base=true; }; if (besoin_des_deformation || besoin_deformation_almansi || besoin_des_deformation_barre || besoin_deformation_greenlagrange || besoin_deformation_logarithmique ) {eps0BB = (NevezTenseurBB(dim_sortie_tenseur)) ; // pour def green-lagrange epsBB = (NevezTenseurBB(dim_sortie_tenseur)) ; // pour def courante epslogBB = (NevezTenseurBB(dim_sortie_tenseur)) ; // pour def logarithmique epsAlmBB = (NevezTenseurBB(dim_sortie_tenseur)) ; // pour def d'almansi epsHB = (NevezTenseurHB(dim)) ; eps_barreHB = (NevezTenseurHB(dim)) ; DeltaEpsBB = (NevezTenseurBB(dim_sortie_tenseur)) ; besoin_matrice_chg_base=true; }; if (besoin_des_valpropre_vitdef || besoin_des_vitesses_deformation ) {DepsBB = (NevezTenseurBB(dim)) ; Deps_barreHB = (NevezTenseurHB(dim)) ; DepsHB = (NevezTenseurHB(dim)) ; besoin_matrice_chg_base=true; }; Tableau tab_ret (enu.size()); // on recupère le tableau pour la lecture des coordonnées des tenseurs int nbcompo = ParaGlob::NbCompTens(); // définition des grandeurs qui sont indépendante de la boucle sur les ddl_enum_etendue def->ChangeNumInteg(iteg); // on change le numéro de point d'intégration courant if (((cas == 1) || (cas == 2))) { // cas d'une premiere initialisation Tableau tab(5); tab(1) = iM0; tab(2) = iMt; tab(3) = iMtdt; tab(4) = igiH_0; tab(5) = igijHH_0; met->PlusInitVariables(tab) ; }; // éléments de métrique et matrices de passage TenseurHH* gijHH=NULL;TenseurBB* gijBB=NULL;BaseB* giB=NULL; BaseH* giH_0=NULL;BaseH* giH=NULL; BaseB* giB_0=NULL; BaseB* giB_t=NULL; // n'est définit "que" pour certains cas Mat_pleine* Aa0 = NULL;Mat_pleine* Aafin = NULL; Mat_pleine* gamma0 = NULL;Mat_pleine* gammafin = NULL; Mat_pleine* beta0 = NULL;Mat_pleine* betafin = NULL; if (besoin_matrice_chg_base) // dans le cas où on n'est pas en absolue => on sort dans un repère ad hoc donc // il a la dimension locale // sinon on sort dans le repère globale => il a la dimension globale {int dim_effective = dim; // init if (absolue) dim_effective = ParaGlob::Dimension(); Aa0 = new Mat_pleine(dim_effective,dim_effective); Aafin = new Mat_pleine(dim_effective,dim_effective); gamma0 = new Mat_pleine(dim_effective,dim_effective); gammafin = new Mat_pleine(dim_effective,dim_effective); beta0 = new Mat_pleine(dim_effective,dim_effective); betafin = new Mat_pleine(dim_effective,dim_effective); }; switch (cas) { case 1: case 11: // calcul d'une ortho interessance de visualisation des tenseurs // cas de tenseur 3D -> Ia, cas 1D on prend un vecteur norme collineaire // avec g1, dans le cas 2D // la nouvelle base gamma est calculee dans def par projection de "Ipa" sur Galpha // le resultat est une matrice de passage utilisable pour le changement de base // Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a // tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer {if (besoin_matrice_chg_base) {const Met_abstraite::InfoImp& ex = def->RemontImp(absolue,*Aa0,*Aafin); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; giB_0 = ex.giB_0; } else {const Met_abstraite::InfoImp& ex = def->RemontImp(); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; giB_0 = ex.giB_0; }; break; } case -1: case -2: // identique au cas 1 mais avec le fait que la métrique est directement disponible // car le calcul est supposé suivre un calcul implicit au bon pti {if (besoin_matrice_chg_base) {Mat_pleine Aat(dim,dim); // a priori Aat ne sert pas par la suite, mais est nécessaire pour le passage de par const Met_abstraite::Info_et_metrique_0_t_tdt ex = def->Remont_et_metrique_0_t_tdtSansCalMet(absolue,*Aa0,Aat,*Aafin); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; giB_0 = ex.giB_0;giB_t = ex.giB_t; } else {const Met_abstraite::Info_et_metrique_0_t_tdt ex = def->Remont_et_metrique_0_t_tdtSansCalMet(); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; giB_0 = ex.giB_0;giB_t = ex.giB_t; } break; } case 2: case 12: // resultat a t {if (besoin_matrice_chg_base) {const Met_abstraite::InfoExp_tdt& ex= def->RemontExp_tdt(absolue,*Aa0,*Aafin); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; giB_0 = ex.giB_0; } else {const Met_abstraite::InfoExp_tdt& ex= def->RemontExp_tdt(); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; giB_0 = ex.giB_0; } break; }; default: cout << "\n *** cas non prevu cas = "<< cas << "\n ElemMeca::Valeur_multi(..." << endl; Sortie(1); }; if (besoin_matrice_chg_base) {// pour les formules de passage de repère il nous faut : // Ip_a = beta_a^{.j} g_j et Ip^b = gamma^b_{.j} g^j // on a: [beta_a^{.j}] = [Aa^j_{.a}]^T // et [gamma^b_{.j}] = [beta_a^{.j}]^{-1T} = [Aa^j_{.a}]^{-1} (*gamma0) = (Aa0->Inverse()); (*gammafin) = (Aafin->Inverse()); // on détermine également les matrices beta (*beta0) = (Aa0->Transpose()); (*betafin) = (Aafin->Transpose()); }; // définition des tenseurs si nécessaire // ----- maintenant on calcule les grandeurs nécessaires ----- // calcul des tenseurs bool plusZero = true; // s'il faut rajouter des termes, on met des 0 if (besoin_des_contraintes) { if (absolue) {(ptIntegMeca.SigHH())->BaseAbsolue(*sigHH,*giB);}// changement de base finale else {(*sigHH) = *(ptIntegMeca.SigHH());sigHH->ChBase(*gammafin);}; (*sigHB) = *(ptIntegMeca.SigHH()) * (*gijBB); }; if (besoin_des_deformation) {(*epsHB) = (*gijHH) * (*(ptIntegMeca.EpsBB())); if (absolue)// changement de base finale {(ptIntegMeca.DeltaEpsBB())->BaseAbsolue((*DeltaEpsBB),*giH); (ptIntegMeca.EpsBB())->BaseAbsolue((*epsBB),*giH);} else {(*DeltaEpsBB) = *(ptIntegMeca.DeltaEpsBB());DeltaEpsBB->ChBase(*betafin); (*epsBB) = *(ptIntegMeca.EpsBB());epsBB->ChBase(*betafin);}; switch (def->Type_de_deformation()) { case DEFORMATION_STANDART : // c'est à dire almansi { // Green-Lagrange if ( besoin_deformation_greenlagrange) { if (absolue) {(ptIntegMeca.EpsBB())->BaseAbsolue((*eps0BB),*giH_0);}// changement de base finale else {(*eps0BB) = *(ptIntegMeca.EpsBB());eps0BB->ChBase(*beta0);}; }; (*epsAlmBB) = (*epsBB); if (epsAlmTotalBB!=NULL) // cas avec dilatation et demande de def Almansi totale {TenseurBB* epsAlmTotal_local_BB = epsAlmTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsAlmTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsAlmTotal_local_BB); if (absolue)// changement de base finale {epsAlmTotal_local_BB->BaseAbsolue(*epsAlmTotalBB,*giH);} else {epsAlmTotalBB->ChBase(*betafin);}; // ici epsAlmTotal_local_BB == epsAlmTotalBB if (prevoir_change_dim_tenseur) delete epsAlmTotal_local_BB; // car pas utilisé ensuite }; if (epsGLTotalBB!=NULL) // cas avec dilatation et demande de def Green_Lagrange totale {TenseurBB* epsGLTotal_local_BB = epsGLTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsGLTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsGLTotal_local_BB); if (absolue)// changement de base finale {epsGLTotal_local_BB->BaseAbsolue(*epsGLTotalBB,*giH_0);} else {epsGLTotalBB->ChBase(*beta0);}; // ici epsGLTotal_local_BB == epsGLTotalBB if (prevoir_change_dim_tenseur) delete epsGLTotal_local_BB; // car pas utilisé ensuite }; // si l'on veut sortir la déformation logarithmique le plus simple est de la calculer if ((besoin_deformation_logarithmique) || (epsLogTotalBB!=NULL)) {def->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE); if (besoin_deformation_logarithmique) // cas du calcul de la def logarithmique {TenseurBB* epslog_local_BB = epslogBB; // par défaut if (prevoir_change_dim_tenseur) epslog_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epslog_local_BB); if (absolue)// changement de base finale {epslog_local_BB->BaseAbsolue((*epslogBB),*giH);} else {epslogBB->ChBase(*betafin);}; // ici epslog_local_BB == epslogBB if (prevoir_change_dim_tenseur) delete epslog_local_BB; // car pas utilisé ensuite }; if (epsLogTotalBB!=NULL) // cas avec dilatation et demande de def log totale {TenseurBB* epsLogTotal_local_BB = epsLogTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsLogTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsLogTotal_local_BB); if (absolue)// changement de base finale {epsLogTotal_local_BB->BaseAbsolue(*epsLogTotalBB,*giH);} else {epsLogTotalBB->ChBase(*betafin);}; // ici epsLogTotal_local_BB == epsLogTotalBB if (prevoir_change_dim_tenseur) delete epsLogTotal_local_BB; // car pas utilisé ensuite }; def->Change_type_de_deformation(DEFORMATION_STANDART); // on revient au type initial }; break; } case DEFORMATION_LOGARITHMIQUE : { (*epslogBB) = (*epsBB); if (epsLogTotalBB!=NULL) // cas avec dilatation et demande de def log totale {TenseurBB* epsLogTotal_local_BB = epsLogTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsLogTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsLogTotal_local_BB); if (absolue)// changement de base finale {epsLogTotal_local_BB->BaseAbsolue(*epsLogTotalBB,*giH);} else {epsLogTotalBB->ChBase(*betafin);}; // ici epsLogTotal_local_BB == epsLogTotalBB if (prevoir_change_dim_tenseur) delete epsLogTotal_local_BB; // car pas utilisé ensuite }; // si l'on veut sortir la déformation d'Almansi ou de green-lagrange le plus simple est de les calculer if (( besoin_deformation_greenlagrange || besoin_deformation_almansi) || (epsAlmTotalBB!=NULL) || (epsGLTotalBB!=NULL)) {def->Change_type_de_deformation(DEFORMATION_STANDART); if ( ( besoin_deformation_greenlagrange || besoin_deformation_almansi) ) // cas de la def d'almansi { TenseurBB* eps_local_BB = epsAlmBB; // par défaut if (prevoir_change_dim_tenseur) eps_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*eps_local_BB); if (absolue)// changement de base finale {eps_local_BB->BaseAbsolue(*epsAlmBB,*giH); eps_local_BB->BaseAbsolue(*eps0BB,*giH_0);} else {epsAlmBB->ChBase(*betafin); eps0BB->ChBase(*beta0);};// ici eps_local_BB == epsAlmBB }; if (epsAlmTotalBB!=NULL) // cas avec dilatation et demande de def Almansi totale {TenseurBB* epsAlmTotal_local_BB = epsAlmTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsAlmTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsAlmTotal_local_BB); if (absolue)// changement de base finale {epsAlmTotal_local_BB->BaseAbsolue(*epsAlmTotalBB,*giH);} else {epsAlmTotalBB->ChBase(*betafin);}; // ici epsAlmTotal_local_BB == epsAlmTotalBB if (prevoir_change_dim_tenseur) delete epsAlmTotal_local_BB; // car pas utilisé ensuite }; if (epsGLTotalBB!=NULL) // cas avec dilatation et demande de def Green_Lagrange totale {TenseurBB* epsGLTotal_local_BB = epsGLTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsGLTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsGLTotal_local_BB); if (absolue)// changement de base finale {epsGLTotal_local_BB->BaseAbsolue(*epsGLTotalBB,*giH_0);} else {epsGLTotalBB->ChBase(*beta0);}; // ici epsGLTotal_local_BB == epsGLTotalBB if (prevoir_change_dim_tenseur) delete epsGLTotal_local_BB; // car pas utilisé ensuite }; def->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE); // on revient au type initial }; break; } default: cout << "\n cas de deformation non encore implante en sortie de visualisation " << Nom_type_deformation(def->Type_de_deformation()) << " affichage donc errone des valeurs !!!"; }; }; if (besoin_des_vitesses_deformation) { if (absolue) {(ptIntegMeca.DepsBB())->BaseAbsolue(*DepsBB,*giH);}// changement de base finale else {(*DepsBB) = *(ptIntegMeca.DepsBB());DepsBB->ChBase(*betafin);}; (*DepsHB) = (*gijHH) * (*(ptIntegMeca.DepsBB())); }; int caas=0; Coordonnee valPropreSig,valPropreEps,valPropreDeps; // init a dim=0 if (besoin_des_valpropre_sigma) {valPropreSig = sigHB->ValPropre(caas); if (caas == -1) { cout << "\n warning *** erreur dans le calcul des valeurs propres de la contrainte (Valeur_multi)"; if (ParaGlob::NiveauImpression() >= 7) {sigHB->Ecriture(cout); cout << "\nElemMeca::Valeur_multi(...";}; cout << endl; // valPropreSig = sigHB.ValPropre(caas); // !!!!!!! pour débug }; }; if (besoin_des_valpropre_deformation) {valPropreEps = epsHB->ValPropre(caas); if (caas == -1) { cout << "\n warning *** erreur dans le calcul des valeurs propres de la deformation"; if (ParaGlob::NiveauImpression() >= 7) {epsHB->Ecriture(cout); cout << "\nElemMeca::Valeur_multi(...";}; cout << endl;}; }; if (besoin_des_valpropre_vitdef) {valPropreDeps = DepsHB->ValPropre(caas); if (caas == -1) { cout << "\n warning *** erreur dans le calcul des valeurs propres de la vitesse de deformation"; if (ParaGlob::NiveauImpression() >= 7) {DepsHB->Ecriture(cout); cout << "\nElemMeca::Valeur_multi(...";}; cout << endl;}; }; if (besoin_coordonnees) {Mtdt = new Coordonnee(ParaGlob::Dimension()); *Mtdt = def->Position_tdt(); } if (besoin_coordonnees_t ) {*Mt = def->Position_tdt(); }; if (besoin_deplacements || besoin_coordonnees_t0) {if (M0 == NULL) M0 = new Coordonnee(ParaGlob::Dimension()); (*M0) = def->Position_0(); }; if (Vitesse != NULL) {Vitesse = new Coordonnee(ParaGlob::Dimension()); (*Vitesse) = def->VitesseM_tdt(); } if (besoin_des_contraintes_barre) {double Isig = sigHB->Trace(); // trace de la déformation (*sig_barreHB) = (*sigHB) - (Isig/ParaGlob::Dimension()) * (*Id_dim_HB(dim)); }; if (besoin_des_deformation_barre) {double Ieps = epsHB->Trace(); // trace de la déformation (*eps_barreHB) = (*epsHB) - (Ieps/ParaGlob::Dimension()) * (*Id_dim_HB(dim)); }; if (besoin_des_vitesse_deformation_barre) {double IDeps = DepsHB->Trace(); // trace de la déformation (*Deps_barreHB) = (*DepsHB) - (IDeps/ParaGlob::Dimension()) * (*Id_dim_HB(dim)); }; // def éventuelle de la contrainte de Mises double Mises = 0.; Coordonnee& vv = valPropreSig; int dimvec=vv.Dimension();// pour condenser l'écriture switch (dimvec) // dans le cas où dimvec=0 on ne fait rien, cas ou on n'a pas besoin de mises { case 1: Mises = Dabs(vv(1)); break; case 2: Mises = sqrt( ((vv(1)-vv(2))*(vv(1)-vv(2)) + vv(1) * vv(1) + vv(2) * vv(2)) * 0.5); break; case 3: Mises = sqrt( ((vv(1)-vv(2))*(vv(1)-vv(2)) + (vv(1)-vv(3))*(vv(1)-vv(3)) + (vv(3)-vv(2))*(vv(3)-vv(2))) * 0.5); break; }; // def éventuelle du vecteur normal: ceci n'est correct qu'avec une métrique 2D if (N_surf != NULL) { // on vérifie que la métrique est correcte if (giB->NbVecteur() != 2) {if (ParaGlob::NiveauImpression() > 0) {cout << "\n *** attention il ne s'agit pas d'un element 2D," << " le vecteur normal ne sera pas disponible"; if (ParaGlob::NiveauImpression() > 2) cout << "\n element: " << Num_elt() << " pti "<< iteg << "\n ElemMeca::Valeur_multi(... "; }; cout << endl; } else // sinon c'est ok {// la normale vaut le produit vectoriel des 2 premiers vecteurs (*N_surf) = Util::ProdVec_coorBN( (*giB)(1), (*giB)(2)); N_surf->Normer(); // que l'on norme }; // on n'arrête pas l'exécution, car on pourrait vouloir sortir les normales pour un ensemble // d'éléments contenant des volumes, des surfaces, des lignes: bon... il y aura quand même des // pb au niveau des iso par exemple, du au fait que l'on va faire des moyennes sur des éléments // de type différents (à moins de grouper par type du coup on n'aura pas le warning }; // idem à l'instant t if (N_surf_t != NULL) { // on vérifie que la métrique est correcte if (giB_t->NbVecteur() != 2) {if (ParaGlob::NiveauImpression() > 0) {cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D," << " le vecteur normal ne sera pas correctement calcule"; if (ParaGlob::NiveauImpression() > 2) cout << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(... "; }; cout << endl; } else // sinon c'est ok {// la normale vaut le produit vectoriel des 2 premiers vecteurs (*N_surf_t) = Util::ProdVec_coorBN( (*giB_t)(1), (*giB_t)(2)); N_surf_t->Normer(); // que l'on norme }; }; // idem à l'instant t0 if (N_surf_t0 != NULL) { // on vérifie que la métrique est correcte if (giB_0->NbVecteur() != 2) {if (ParaGlob::NiveauImpression() > 0) {cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D," << " le vecteur normal ne sera pas correctement calcule"; if (ParaGlob::NiveauImpression() > 2) cout << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(... "; }; cout << endl; } else // sinon c'est ok {// la normale vaut le produit vectoriel des 2 premiers vecteurs (*N_surf_t0) = Util::ProdVec_coorBN( (*giB_0)(1), (*giB_0)(2)); N_surf_t0->Normer(); // que l'on norme }; }; // def éventuelle de la déformation duale de mises = sqrt(2/3 * epsB:epsB) //--- on utilise directement la grandeur stockée au pt d'integ, mais ici ne sert à rien, puis cette grandeur n'est pas utilisée par la suite ! // // l'expression est la même que celle de mises, ormis un coeff 4/9 qui permet de passer de 3/2 à 2/3 // double defDualMises = 0.; Coordonnee& vdef = valPropreEps; int dimvdef=vdef.Dimension();// pour condenser l'écriture // switch (dimvdef) // dans le cas où dimvdef=0 on ne fait rien, cas ou on n'a pas besoin de defDualMises // { case 1: defDualMises = Dabs(vdef(1)); break; // case 2: defDualMises = sqrt( ((vdef(1)-vdef(2))*(vdef(1)-vdef(2)) + vdef(1) * vdef(1) // + vdef(2) * vdef(2)) * 2./9.); break; // case 3: defDualMises = sqrt( ((vdef(1)-vdef(2))*(vdef(1)-vdef(2)) + (vdef(1)-vdef(3))*(vdef(1)-vdef(3)) // + (vdef(3)-vdef(2))*(vdef(3)-vdef(2))) * 2./9.); break; // }; // donnees propre a la loi mécanique au pt d'integ Loi_comp_abstraite::SaveResul * sDon = tabSaveDon(iteg); // donnees propre a la loi thermo physique au pt d'integ CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques if (loiTP != NULL) {sTP = tabSaveTP(iteg);}; // au pt d'integ si elles existes // donnees propre a la déformation mécanique au pt d'integ Deformation::SaveDefResul * sDefDon = tabSaveDefDon(iteg); // pour la sortie des grandeurs polaires (def et contrainte) double mini_Q = 5.e-5; double * Q_eps=NULL,* Q_sig=NULL,* Q_Deps; //----- fin du calcul des grandeurs nécessaires ----- // on balaie maintenant la liste des grandeurs à sortir int it; // it est l'indice dans le tableau de retour for (it=1,ie=enu.begin(); ie!=iefin;ie++,it++) { // dans le cas où c'est une contrainte, une déformation ou d'une vitesse de déformation // il y a préparation des grandeurs à sortir if ((Meme_famille((*ie).Enum(),SIG11)) || (Meme_famille((*ie).Enum(),EPS11)) || (Meme_famille((*ie).Enum(),DEPS11)) || (Meme_famille((*ie).Enum(),X1)) || (Meme_famille((*ie).Enum(),UX)) ) {// def du numéro de référence du ddl_enum_etendue int posi = (*ie).Position()-NbEnum_ddl(); // récupération des informations en fonction des différents cas // **** 1 >>>>> -- cas des ddl pur, que l'on sort dans le repère global par défaut // cas des contraintes if ((Meme_famille((*ie).Enum(),SIG11)) && ((*ie).Nom_vide())) { // récup de l'ordre Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*sigHH)(ij.i,ij.j); } else if ((Meme_famille((*ie).Enum(),EPS11)) && ((*ie).Nom_vide())) { // récup de l'ordre Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*epsBB)(ij.i,ij.j); } else if ((Meme_famille((*ie).Enum(),DEPS11)) && ((*ie).Nom_vide())) { // récup de l'ordre Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*DepsBB)(ij.i,ij.j); } else if ((Meme_famille((*ie).Enum(),X1)) && ((*ie).Nom_vide())) { tab_ret(it)= (*Mtdt)((*ie).Enum() - X1 +1); } else if ((Meme_famille((*ie).Enum(),UX)) && ((*ie).Nom_vide())) { int i_cor = (*ie).Enum() - UX +1; // l'indice de coordonnée tab_ret(it)= (*Mtdt)(i_cor) - (*M0)(i_cor); } else if ((Meme_famille((*ie).Enum(),V1)) && ((*ie).Nom_vide())) { int i_cor = (*ie).Enum() - V1 +1; // l'indice de coordonnée tab_ret(it)= (*Vitesse)(i_cor); } // --- a complèter ---- else {// **** 2 >>>>> -- cas des grandeurs déduites des ddl pures switch (posi) { case 1: case 2: case 3: case 4: case 5: case 6: /*Green-Lagrange */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*eps0BB)(ij.i,ij.j);break;} case 7: case 8: case 9: case 10: case 11: case 12: /*Almansi */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*epsAlmBB)(ij.i,ij.j);break;} case 49: case 50: case 51: case 52: case 53: case 54: /*logarithmique */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*epslogBB)(ij.i,ij.j);break;} case 55: case 56: case 57: case 58: case 59: case 60: /*Almansi totale */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*epsAlmTotalBB)(ij.i,ij.j);break;} case 61: case 62: case 63: case 64: case 65: case 66: /*Green_Lagrange totale */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*epsGLTotalBB)(ij.i,ij.j);break;} case 67: case 68: case 69: case 70: case 71: case 72: /*Log totale */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*epsLogTotalBB)(ij.i,ij.j);break;} case 13: case 14: case 15: case 16: case 17: case 18: /*Cauchy_local */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*ptIntegMeca.SigHH())(ij.i,ij.j);break;} case 19: case 20: case 21: case 22: case 23: case 24: /*Almansi_local */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*ptIntegMeca.EpsBB())(ij.i,ij.j);break; } case 25: /*Def_principaleI*/ tab_ret(it)=valPropreEps(1);break; case 26: /*Def_principaleII*/ if (valPropreEps.Dimension() > 1) {tab_ret(it)=valPropreEps(2);} else {tab_ret(it)= 0.;};break; case 27: /*Def_principaleIII*/ if (valPropreEps.Dimension() > 2) {tab_ret(it)=valPropreEps(3);} else {tab_ret(it)= 0.;};break; case 28: /*Sigma_principaleI*/ tab_ret(it)=valPropreSig(1);break; case 29: /*Sigma_principaleII*/ if (valPropreSig.Dimension() > 1) {tab_ret(it)=valPropreSig(2);} else {tab_ret(it)= 0.;};break; case 30: /*Sigma_principaleIII*/ if (valPropreSig.Dimension() > 2) {tab_ret(it)=valPropreSig(3);} else {tab_ret(it)= 0.;};break; case 31: /*contrainte_mises*/ tab_ret(it)=Mises;break; // case 77: /*def_duale_mises*/ tab_ret(it)=defDualMises;break; case 77: /*def_duale_mises*/ tab_ret(it)=ptIntegMeca.Deformation_equi_const()(2);break; case 87: /*def_equivalente*/ tab_ret(it)=ptIntegMeca.Deformation_equi_const()(1);break; case 88: /*def_duale_mises_maxi*/ tab_ret(it)=ptIntegMeca.Deformation_equi_const()(3);break; case 89: /*vitesse_def_equivalente*/ tab_ret(it)=ptIntegMeca.Deformation_equi_const()(4) * unSurDeltat;break; case 32: /*contrainte_tresca*/ { switch (dim) {case 1: tab_ret(it)=0.5 * valPropreSig(1);break; case 2: tab_ret(it)=0.5 * (valPropreSig(1)-valPropreSig(2));break; case 3: tab_ret(it)=0.5 * (valPropreSig(1)-valPropreSig(3));break; } break; } case 33: /*def_plastique_cumulee*/ { string nom_comport(loiComp->Nom_comport()); if ((strstr(nom_comport.c_str(),"PRANDTL_REUSS")!=NULL)) { tab_ret(it) = tabSaveDon(iteg)->Deformation_plastique(); } else { cout << "\n erreur, la déformation plastique n'est pas disponible" << "dans l element " << this->Num_elt() << "\n ElemMeca::Valeur_multi(...."; Sortie(1); } break; } case 34: case 35: case 36: case 37: case 38: case 39: /*Vit_def11 et suite */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*DepsBB)(ij.i,ij.j);break; } case 40: /*Vit_principaleI*/ tab_ret(it)=valPropreDeps(1);break; case 41: /*Vit_principaleII*/ if (valPropreDeps.Dimension() > 1) {tab_ret(it)=valPropreDeps(2);} else {tab_ret(it)= 0.;};break; case 42: /*Vit_principaleIII*/ tab_ret(it)=valPropreDeps(3);break; if (valPropreDeps.Dimension() > 2) {tab_ret(it)=valPropreDeps(3);} else {tab_ret(it)= 0.;};break; case 43: case 44: case 45: case 46: case 47: case 48: /*Delta_def11 et suite */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); tab_ret(it)= (*DeltaEpsBB)(ij.i,ij.j);break; } case 73: /*energie_elastique*/ if (temps==TEMPS_tdt) {tab_ret(it)=tab_energ(iteg).EnergieElastique();} else {tab_ret(it)=tab_energ_t(iteg).EnergieElastique();};break; case 74: /*dissipation_plastique*/ if (temps==TEMPS_tdt) {tab_ret(it)=tab_energ(iteg).DissipationPlastique();} else {tab_ret(it)=tab_energ_t(iteg).DissipationPlastique();};break; case 75: /*dissipation_visqueuse*/ if (temps==TEMPS_tdt) {tab_ret(it)=tab_energ(iteg).DissipationVisqueuse();} else {tab_ret(it)=tab_energ_t(iteg).DissipationVisqueuse();};break; case 78: /*Spherique_eps*/ tab_ret(it)=epsHB->Trace()/3.; break; //ParaGlob::Dimension();break; modi du 5/2/2012 case 79: /*Q_eps*/ {Q_eps = new double; tab_ret(it)=*Q_eps = sqrt(eps_barreHB->II());break;} case 80: /*Cos3phi_eps*/ { double Qepsilon = ( (Q_eps!=NULL) ? *Q_eps : sqrt(eps_barreHB->II())); double Qepsilon3 = Qepsilon * Qepsilon * Qepsilon; if (Qepsilon > mini_Q ) { // on peut calculer un cos3phi pas débile double bIIIb = eps_barreHB->III() / 3.; tab_ret(it) = 3. * sqrt(6.) * bIIIb/ Qepsilon3; ////------ debug //cout << "\n debug: ElemMeca::Valeur_multi " // << "\n Qepsilon3= "<< Qepsilon3 << " bIIIb= "<< bIIIb // << " Cos3phi_eps= " << tab_ret(it) << " "; // ////--- fin debug } else tab_ret(it)=0.; // sinon on le met à 0 break; } case 81: /*Spherique_sig*/ tab_ret(it)=sigHB->Trace()/ParaGlob::Dimension();;break; case 82: /*Q_sig*/ {Q_sig = new double; tab_ret(it)=*Q_sig = sqrt(sig_barreHB->II());break;} case 83: /*Cos3phi_sig*/ { double Qsig = ( (Q_sig!=NULL) ? *Q_sig : sqrt(sig_barreHB->II())); double Qsig3 = Qsig * Qsig * Qsig; if (Qsig > mini_Q ) { // on peut calculer un cos3phi pas débile double bIIIb = sig_barreHB->III() / 3.; tab_ret(it) = 3. * sqrt(6.) * bIIIb/ Qsig3; } else tab_ret(it)=0.; // sinon on le met à 0 break; } case 84: /*Spherique_Deps*/ tab_ret(it)=DepsHB->Trace()/ParaGlob::Dimension();break; case 85: /*Q_Deps*/ {Q_Deps = new double; tab_ret(it)=*Q_Deps = sqrt(Deps_barreHB->II());break;} case 86: /*Cos3phi_Deps*/ { double QDepsilon = ( (Q_Deps!=NULL) ? *Q_Deps : sqrt(Deps_barreHB->II())); double QDepsilon3 = QDepsilon * QDepsilon * QDepsilon; if (QDepsilon > mini_Q ) { // on peut calculer un cos3phi pas débile double bIIIb = Deps_barreHB->III() / 3.; tab_ret(it) = 3. * sqrt(6.) * bIIIb/ QDepsilon3; } else tab_ret(it)=0.; // sinon on le met à 0 break; } // le cas 94 est a priori ok, mais il y a un pb de logique: en fait les pti d'efforts externes ne sont pas // a priori les mêmes que ceux de la raideur, donc cela va entrainer des pb // il faudrait une autre méthode spécifique aux efforts externes. // je laisse le code en prévision, mais je commente pour éviter les pb // case 94: /* "pression_ext */ // { if (lesChargeExtSurEle != NULL) // {if (lesChargeExtSurEle->LesPressionsExternes() != NULL) // cas où il existe des pressions sauvegardées // { Tableau >& lesPressionsExternes = *lesChargeExtSurEle->LesPressionsExternes(); // int nb_face = lesPressionsExternes.Taille(); // if (nb_face != 1) // {cout << "\n pas de sortie de pression possible en dehors d'element face "<< endl; // tab_ret(it)= 0.; // } // else // {int n_face=1; // int t_tail = lesPressionsExternes(n_face).Taille(); // if (t_tail != 0) // { Tableau & tab_press_appliquee = (lesPressionsExternes(n_face)); // pour simplifier // if (t_tail != 0) // cas où on a une face chargée // { tab_ret(it)=tab_press_appliquee(it).press;} // else {tab_ret(it)= 0.;}; // sinon rien // } // else {tab_ret(it)= 0.;}; // sinon rien // }; // } // else {tab_ret(it)= 0.;}; // sinon rien // } // else {tab_ret(it)= 0.;}; // sinon rien // break; // } case 114: // le vecteur normal N_surf_1 {tab_ret(it)= (*N_surf)(1);break;} case 115: // le vecteur normal N_surf_2 {tab_ret(it)= (*N_surf)(2);break;} case 116: // le vecteur normal N_surf_3 {tab_ret(it)= (*N_surf)(3);break;} case 117: // le vecteur normal N_surf_1_t {tab_ret(it)= (*N_surf_t)(1);break;} case 118: // le vecteur normal N_surf_2_t {tab_ret(it)= (*N_surf_t)(2);break;} case 119: // le vecteur normal N_surf_3_t {tab_ret(it)= (*N_surf_t)(3);break;} case 120: // le vecteur normal N_surf_1_t0 {tab_ret(it)= (*N_surf_t0)(1);break;} case 121: // le vecteur normal N_surf_2_t0 {tab_ret(it)= (*N_surf_t0)(2);break;} case 122: // le vecteur normal N_surf_3_t0 {tab_ret(it)= (*N_surf_t0)(3);break;} case 123: // la position géométrique Mt {tab_ret(it)= (*Mt)(1);break;} case 124: // la position géométrique Mt {tab_ret(it)= (*Mt)(2);break;} case 125: // la position géométrique Mt {tab_ret(it)= (*Mt)(3);break;} case 126: // la position géométrique M0 {tab_ret(it)= (*M0)(1);break;} case 127: // la position géométrique M0 {tab_ret(it)= (*M0)(2);break;} case 128: // la position géométrique M0 {tab_ret(it)= (*M0)(3);break;} case 137: // l'erreur relative {if (sigErreur_relative != NULL) {tab_ret(it)= *sigErreur_relative;} else {tab_ret(it)= 0.;}; break; } default : {cout << "\n cas de ddl actuellement non traite " << "\n pas de ddl " << (*ie).Nom() << " dans l'element " << "\n ou cas non implante pour l'instant" << "\n ElemMeca::Valeur_multi(...."; tab_ret(it) = 0.; }; } // fin cas **** 2 >>>>> } // " " " } // -- fin du cas ou c'est une grandeur liée aux contraintes ou déformations // cas de l'erreur else if (( (*ie).Enum() == ERREUR) && ((*ie).Nom_vide())) {if (sigErreur != NULL) tab_ret(it)= *sigErreur; else if (ParaGlob::NiveauImpression()>4) {cout << "\n pas encore de ddl erreur dans l'element " << "\n ElemMeca::Valeur_multi(...."; // this->Affiche(); }; } else if (( (*ie).Enum() == TEMP) && ((*ie).Nom_vide())) {// on vérifie que le ddl existe au premier noeud if (tab_noeud(1)->Existe_ici(TEMP)) {tab_ret(it)= def->DonneeInterpoleeScalaire(TEMP,temps);} else if (ParaGlob::NiveauImpression()>3) {cout << "\n pas de ddl temperature disponible au noeud " << "\n ElemMeca::Valeur_multi(...."; // this->Affiche(); }; } else { tab_ret(it) = 0.; cout << "\n cas de ddl actuellement non traite " << "\n pas de ddl " << (*ie).Nom() << " dans l'element " << "\n ou cas non implante pour l'instant, on retourne 0" << "\n ElemMeca::Valeur_multi(...."; }; };// -- fin de la boucle sur la liste de Ddl_enum_etendu // delete des tenseurs if (sigHH != NULL); delete sigHH; if (eps0BB != NULL); delete eps0BB; if (epsBB != NULL); delete epsBB; if (epslogBB != NULL); delete epslogBB; if (epsAlmBB != NULL); delete epsAlmBB; if (sigHB != NULL); delete sigHB; if (epsHB != NULL); delete epsHB; if (DepsBB != NULL); delete DepsBB; if (DepsHB != NULL); delete DepsHB; if (DeltaEpsBB != NULL); delete DeltaEpsBB; if (eps_barreHB != NULL); delete eps_barreHB; if (Deps_barreHB != NULL); delete Deps_barreHB; if (sig_barreHB != NULL); delete sig_barreHB; // cas des pointeurs if (epsAlmTotalBB!=NULL) delete epsAlmTotalBB; // pour la déformation totale d'almansi if (epsGLTotalBB!=NULL) delete epsGLTotalBB; // pour la déformation totale de green_lagrange if (epsLogTotalBB!=NULL) delete epsLogTotalBB; // pour la déformation totale logarithmique if (Q_sig != NULL) delete Q_sig; // grandeurs polaires if (Q_eps != NULL) delete Q_eps; // grandeurs polaires if (Mtdt != NULL) delete Mtdt; // coordonnée du point à t if (Mt != NULL ) delete Mt; // la position à t if (M0 != NULL ) delete M0; // coordonnée du point à 0 if (N_surf != NULL) delete N_surf; // vecteur normal à la surface if (N_surf_t != NULL) delete N_surf_t; // vecteur normal à la surface à t if (N_surf_t0 != NULL) delete N_surf_t0; // vecteur normal à la surface à t0 if (Vitesse != NULL) delete Vitesse; // vitesse // pointeurs de matrice if (Aa0 != NULL) delete Aa0; if (Aafin != NULL) delete Aafin; if (gamma0 != NULL) delete gamma0; if (gammafin != NULL) delete gammafin; if (beta0 != NULL) delete beta0; if (betafin != NULL) delete betafin; // liberation des tenseurs intermediaires LibereTenseur(); def->Retour_pti_precedant(); // on revient au pti précédent return tab_ret; }; // récupération des valeurs Tensorielles (et non scalaire comme avec Valeur_multi) // au numéro d'ordre = iteg pour les grandeur enu // enu contiend les grandeurs de retour // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière // cas = -1 ou -2: cela signifie que la métrique vient juste d'être calculée au pti iteg // sinon il faut recalculer qq éléments de la métrique // NB Importante: il faut faire attention à ce que ces métriques soient identiques à celles qui ont servit // pour le calcul des tenseurs: en particulier si c'est utilisé pour calculer les grandeurs pour le chargement // il faut s'assurer que ce sont les "mêmes pti" qui servent pour la charge et pour la raideur !! void ElemMeca::Valeurs_Tensorielles(bool absolue, Enum_dure temps,List_io& enu ,int iteg,int cas ) // --- la méthode est pratiquement la même que celle dans Loi_comp_abstraite // ici il y a le paramètre cas qui permet de flêcher sur une métrique particulière // Dans loi_comp_abstraite on récupère en paramètre, directement la métrique éventuellement ... // ici on recalcule Mises et Tresca, dans Loi_comp_abstraite on ressort les valeurs déjà calculée // ici on { // ----- def de grandeurs de travail // def de la dimension des tenseurs // il y a deux pb a gérer: 1) le fait que la dimension absolue peut-être différente de la dimension des tenseurs // 2) le fait que l'on veut une sortie dans une base ad hoc ou pas int dim = lesPtIntegMecaInterne->DimTens();int dim_sortie_tenseur = dim; // dans le cas ou l'on veut une sortie en base absolue, il faut que dim_sortie_tenseur = la dimension de la base absolue int dim_espace = ParaGlob::Dimension(); if (absolue) dim_sortie_tenseur = dim_espace; // pour ne faire qu'un seul test ensuite bool prevoir_change_dim_tenseur = false; // initialement on faisait le test suivant, // if ((absolue) && (dim != dim_sortie_tenseur)) if (absolue) // mais en fait, quand on sort en absolu on est obligé d'utiliser un tenseur intermédiaire // car BaseAbsolue(.. modifie tenseur passé en paramètre, donc dans tous les cas de sortie absolue // il faut un tenseur intermédiaire qui a ou non une dimension différente prevoir_change_dim_tenseur = true; PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(iteg); // def de tenseurs pour la sortie TenseurHH* sigHH=NULL;TenseurBB* epsBB=NULL;TenseurBB* epslogBB=NULL; TenseurBB* epsAlmBB=NULL;TenseurBB* eps0BB=NULL;TenseurBB* DepsBB = NULL; TenseurBB* epsAlmTotalBB=NULL;TenseurBB* epsGLTotalBB=NULL;TenseurBB* epsLogTotalBB=NULL; TenseurBB* DeltaEpsBB = NULL; Coordonnee* Mtdt=NULL; // coordonnées éventuelles du point d'intégration considéré Coordonnee* Mt=NULL; // coordonnées à t Coordonnee* M0=NULL; // coordonnées à t0 Coordonnee* N_surf = NULL; // coordonnée d'un vecteur normal si c'est adéquate Coordonnee* N_surf_t = NULL; // coordonnée d'un vecteur normal à t si c'est adéquate Coordonnee* N_surf_t0 = NULL; // coordonnée d'un vecteur normal à t0 si c'est adéquate Coordonnee* Vitesse = NULL; // cas des vitesses Coordonnee* Deplacement = NULL; // cas du déplacement // pour les valeurs propres Coordonnee* valPropreSig=NULL;Coordonnee* valPropreEps=NULL; Coordonnee* valPropreDeps=NULL; // pour les vecteurs propres Tableau base_propre_sigma , base_propre_eps , base_propre_D; // pour les bases Tableau base_ad_hoc , base_giH , base_giB; // grandeurs scalaires double* Mises = NULL; double* defDualMises = NULL; double* Tresca = NULL; double* erreur = NULL;double* erreur_rel = NULL; double* spherique_eps=NULL,* Q_eps=NULL,* cos3phi_eps=NULL; double* spherique_sig=NULL,* Q_sig=NULL,* cos3phi_sig=NULL; double* spherique_Deps=NULL,* Q_Deps=NULL,* cos3phi_Deps=NULL; double* def_equivalente=NULL, * defDualMisesMaxi=NULL; // --- dev d'un ensemble de variable booléenne pour gérer les sorties en une passe ----- // on se réfère au informations définit dans la méthode: Les_type_evolues_internes() bool deformationCourante=false; bool vitesseDeformationCourante=false; bool almansi=false; bool greenLagrange=false; bool logarithmique=false; bool deltaDef=false; bool almansiTotal=false; bool greenLagrangeTotale=false; bool logarithmiqueTotale=false; bool contrainteCourante=false; bool defPrincipales=false; bool sigmaPrincipales=false; bool vitPrincipales=false; bool contrainteMises=false; bool contraintesTresca=false; bool erreurQ=false; bool erreurRel = false; bool defPlastiqueCumulee=false; bool def_duale_mises=false; bool besoin_des_contraintes=false; bool besoin_des_deformation=false; bool besoin_des_contraintes_barre=false; bool besoin_des_deformation_barre=false; bool besoin_des_vitesses_deformation=false; bool besoin_des_vitesses_deformation_barre=false; bool besoin_des_valpropre_sigma=false; bool besoin_des_valpropre_deformation = false; bool besoin_des_valpropre_vitdef = false; bool besoin_coordonnees = false; bool besoin_coordonnees_t = false;bool besoin_coordonnees_t0 = false; bool besoin_dir_princ_sig = false; bool besoin_dir_princ_eps = false; bool besoin_dir_princ_D = false; bool besoin_rep_local_ortho=false; bool besoin_rep_giH = false; bool besoin_rep_giB = false; bool def_SPHERIQUE_EPS = false; bool def_Q_EPS = false; bool def_COS3PHI_EPS = false; bool def_SPHERIQUE_SIG = false; bool def_Q_SIG = false; bool def_COS3PHI_SIG = false; bool def_SPHERIQUE_DEPS = false; bool def_Q_DEPS = false; bool def_COS3PHI_DEPS = false; bool def_def_equivalente = false; bool def_duale_mises_maxi = false; // on initialise ces variables booléennes et les conteneurs List_io::iterator ipq,ipqfin=enu.end(); for (ipq=enu.begin();ipq!=ipqfin;ipq++) {switch ((*ipq).EnuTypeQuelconque().EnumTQ()) { case CONTRAINTE_COURANTE : {contrainteCourante=true; besoin_des_contraintes=true; Grandeur_TenseurHH& gr= *((Grandeur_TenseurHH*) ((*ipq).Grandeur_pointee())); sigHH = gr.ConteneurTenseur(); break;} case DEFORMATION_COURANTE : {deformationCourante=true; besoin_des_deformation=true; Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee())); epsBB = gr.ConteneurTenseur(); break;} case VITESSE_DEFORMATION_COURANTE : {vitesseDeformationCourante=true; Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee())); DepsBB = gr.ConteneurTenseur(); besoin_des_vitesses_deformation=true; break;} case ALMANSI : {almansi=true; besoin_des_deformation=true; Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee())); epsAlmBB = gr.ConteneurTenseur(); break;} case GREEN_LAGRANGE : {greenLagrange=true; besoin_des_deformation=true; Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee())); eps0BB = gr.ConteneurTenseur(); break;} case LOGARITHMIQUE : {logarithmique=true; besoin_des_deformation=true; Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee())); epslogBB = gr.ConteneurTenseur(); break;} case DELTA_DEF : {deltaDef=true; besoin_des_deformation=true; Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee())); DeltaEpsBB = gr.ConteneurTenseur(); break;} case ALMANSI_TOTAL : {almansiTotal=true; besoin_des_deformation=true; Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee())); epsAlmTotalBB = gr.ConteneurTenseur(); break;} case GREEN_LAGRANGE_TOTAL : {greenLagrangeTotale=true; besoin_des_deformation=true; Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee())); epsGLTotalBB = gr.ConteneurTenseur(); break;} case LOGARITHMIQUE_TOTALE : {logarithmiqueTotale=true; besoin_des_deformation=true; Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee())); epsLogTotalBB = gr.ConteneurTenseur(); break;} case DEF_PRINCIPALES : {defPrincipales=true; besoin_des_deformation=true; besoin_des_valpropre_deformation=true; Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); valPropreEps = gr.ConteneurCoordonnee(); break;} case SIGMA_PRINCIPALES : {sigmaPrincipales=true; besoin_des_contraintes=true; besoin_des_valpropre_sigma=true; Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); valPropreSig = gr.ConteneurCoordonnee(); break;} case VIT_PRINCIPALES : {vitPrincipales=true; besoin_des_vitesses_deformation=true; besoin_des_valpropre_vitdef=true; Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); valPropreDeps = gr.ConteneurCoordonnee(); break;} case CONTRAINTE_MISES : {contrainteMises=true; besoin_des_contraintes=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); Mises = gr.ConteneurDouble(); besoin_des_valpropre_sigma=true; break;} case DEF_DUALE_MISES : {def_duale_mises=true; besoin_des_deformation=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); defDualMises = gr.ConteneurDouble(); //besoin_des_valpropre_deformation=true; // non car maintenant on utilise directement ce qui est stocké au pt d'integ break;} case CONTRAINTE_TRESCA : {contraintesTresca=true; besoin_des_contraintes=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); Tresca = gr.ConteneurDouble(); besoin_des_valpropre_sigma=true; break;} case ERREUR_Q : {erreurQ=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); erreur = gr.ConteneurDouble(); break;} case ERREUR_SIG_RELATIVE : {erreurRel=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); erreur_rel = gr.ConteneurDouble(); break;} case POSITION_GEOMETRIQUE : {besoin_coordonnees=true; Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); Mtdt = gr.ConteneurCoordonnee(); break;} case POSITION_GEOMETRIQUE_t : {besoin_coordonnees_t=true; Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); Mt = gr.ConteneurCoordonnee(); break;} case POSITION_GEOMETRIQUE_t0 : {besoin_coordonnees_t0=true; Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); M0 = gr.ConteneurCoordonnee(); break;} case SPHERIQUE_EPS : {def_SPHERIQUE_EPS=true; besoin_des_deformation=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); spherique_eps = gr.ConteneurDouble(); break;} case Q_EPS : {def_Q_EPS=true; besoin_des_deformation=true;besoin_des_deformation_barre=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); Q_eps = gr.ConteneurDouble(); break;} case COS3PHI_EPS : {def_COS3PHI_EPS=true; besoin_des_deformation=true;besoin_des_deformation_barre=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); cos3phi_eps = gr.ConteneurDouble(); besoin_des_valpropre_deformation=true; break;} case SPHERIQUE_SIG : {def_SPHERIQUE_SIG=true; besoin_des_contraintes=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); spherique_sig = gr.ConteneurDouble(); besoin_des_valpropre_sigma=true; break;} case Q_SIG : {def_Q_SIG=true; besoin_des_contraintes=true;besoin_des_contraintes_barre=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); Q_sig = gr.ConteneurDouble(); break;} case COS3PHI_SIG : {def_COS3PHI_SIG=true; besoin_des_contraintes=true;besoin_des_contraintes_barre=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); cos3phi_sig = gr.ConteneurDouble(); break;} case SPHERIQUE_DEPS : {def_SPHERIQUE_DEPS=true; besoin_des_vitesses_deformation=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); spherique_Deps = gr.ConteneurDouble(); break;} case Q_DEPS : {def_Q_DEPS=true; besoin_des_vitesses_deformation=true;besoin_des_vitesses_deformation_barre=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); Q_Deps = gr.ConteneurDouble(); break;} case COS3PHI_DEPS : {def_COS3PHI_DEPS=true; besoin_des_vitesses_deformation=true;besoin_des_vitesses_deformation_barre=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); cos3phi_Deps = gr.ConteneurDouble(); break;} case DEF_EQUIVALENTE : {def_def_equivalente=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); def_equivalente = gr.ConteneurDouble(); break;} case DEF_DUALE_MISES_MAXI : {def_duale_mises_maxi=true; Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); defDualMisesMaxi = gr.ConteneurDouble(); break;} case REPERE_LOCAL_ORTHO : {besoin_rep_local_ortho=true; Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee())); base_ad_hoc.Change_taille(dim_espace); switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init switch (dim_espace) {case 3: base_ad_hoc(3) = &(gr(3)); case 2: base_ad_hoc(2) = &(gr(2)); case 1: base_ad_hoc(1) = &(gr(1)); default:break; }; break;} case REPERE_LOCAL_H : {besoin_rep_giH=true; Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee())); base_giH.Change_taille(dim); switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init switch (dim) {case 3: base_giH(3) = &(gr(3)); case 2: base_giH(2) = &(gr(2)); case 1: base_giH(1) = &(gr(1)); default:break; }; break;} case REPERE_LOCAL_B : {besoin_rep_giB=true; Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee())); base_giB.Change_taille(dim); switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init switch (dim) {case 3: base_giB(3) = &(gr(3)); case 2: base_giB(2) = &(gr(2)); case 1: base_giB(1) = &(gr(1)); default:break; }; break;} case DIRECTIONS_PRINC_SIGMA : {besoin_des_contraintes=true; besoin_des_valpropre_sigma=true;besoin_dir_princ_sig=true; Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee())); switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init base_propre_sigma.Change_taille(dim); switch (dim) {case 3: base_propre_sigma(3) = &(gr(3)); case 2: base_propre_sigma(2) = &(gr(2)); case 1: base_propre_sigma(1) = &(gr(1)); default:break; }; break;} case DIRECTIONS_PRINC_DEF : {besoin_des_deformation=true; besoin_des_valpropre_deformation=true;besoin_dir_princ_eps=true; Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee())); switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init base_propre_eps.Change_taille(dim); switch (dim) {case 3: base_propre_eps(3) = &(gr(3)); case 2: base_propre_eps(2) = &(gr(2)); case 1: base_propre_eps(1) = &(gr(1)); default:break; }; break;} case DIRECTIONS_PRINC_D : {vitesseDeformationCourante=true;besoin_des_vitesses_deformation=true; besoin_des_valpropre_deformation=true;besoin_dir_princ_D=true; Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee())); switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init base_propre_D.Change_taille(dim); switch (dim) {case 3: base_propre_D(3) = &(gr(3)); case 2: base_propre_D(2) = &(gr(2)); case 1: base_propre_D(1) = &(gr(1)); default:break; }; break;} case NN_SURF:{ Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); N_surf = gr.ConteneurCoordonnee(); break;} case NN_SURF_t:{ Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); N_surf_t = gr.ConteneurCoordonnee(); break;} case NN_SURF_t0:{ Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); N_surf_t0 = gr.ConteneurCoordonnee(); break;} case DEPLACEMENT:{besoin_coordonnees=true; Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); Deplacement = gr.ConteneurCoordonnee(); break;} case VITESSE:{ Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); Vitesse = gr.ConteneurCoordonnee(); break;} case ACCELERATION:{ Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); Vitesse = gr.ConteneurCoordonnee(); break;} // dans le cas des numéros, traitement direct ici case NUM_ELEMENT: { *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()))= num_elt; break;} case NUM_MAIL_ELEM: { *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()))= num_maillage; break;} case NUM_PTI: { *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()))= iteg; break;} default : {// on initialise la grandeur pour éviter d'avoir des valeurs aléatoires ((*ipq).Grandeur_pointee())->InitParDefaut(); if (ParaGlob::NiveauImpression() > 0) {cout << "\nWarning : attention cas non traite: " << (*ipq).EnuTypeQuelconque().NomPlein() << "!\n"; // on initialise la grandeur pour éviter d'avoir des valeurs aléatoires if (ParaGlob::NiveauImpression() > 5) cout << "\n ElemMeca::Valeurs_Tensorielles(...."; }; } }; }; // on complète les définitions de tenseurs si besoin // les tenseurs restants en locale TenseurHB* sigHB = NULL ;TenseurHB* sig_barreHB = NULL ; TenseurHB* epsHB = NULL ;TenseurHB* eps_barreHB = NULL ; TenseurHB* DepsHB = NULL ; TenseurHB* Deps_barreHB = NULL ; if ((besoin_des_contraintes && !contrainteCourante) || (besoin_des_valpropre_sigma && !sigmaPrincipales) ) {sigHH = NevezTenseurHH(dim_sortie_tenseur) ; sigHB = (NevezTenseurHB(dim)) ; sig_barreHB = (NevezTenseurHB(dim)) ; } else if (besoin_des_contraintes || (besoin_des_valpropre_sigma && !sigmaPrincipales) ) // là il s'agit uniquement des tenseurs locaux {sigHB = (NevezTenseurHB(dim)) ; sig_barreHB = (NevezTenseurHB(dim)) ; }; if ((besoin_des_deformation && !deformationCourante) || (besoin_des_valpropre_deformation && !defPrincipales) ) {epsBB = NevezTenseurBB(dim_sortie_tenseur) ; epsHB = (NevezTenseurHB(dim)) ; eps_barreHB = (NevezTenseurHB(dim)) ; } else if (besoin_des_deformation || (besoin_des_valpropre_deformation && !defPrincipales) ) // là il s'agit uniquement des tenseurs locaux {epsHB = (NevezTenseurHB(dim)) ; eps_barreHB = (NevezTenseurHB(dim)) ; }; if ((besoin_des_vitesses_deformation && !vitesseDeformationCourante) || (besoin_des_valpropre_vitdef && !vitPrincipales) ) {DepsBB = NevezTenseurBB(dim_sortie_tenseur); DepsHB = (NevezTenseurHB(dim)); Deps_barreHB = (NevezTenseurHB(dim)) ; } else if(besoin_des_vitesses_deformation || (besoin_des_valpropre_vitdef && !vitPrincipales) ) // là il s'agit uniquement des tenseurs locaux {DepsHB = (NevezTenseurHB(dim)); Deps_barreHB = (NevezTenseurHB(dim)) ; }; if (besoin_des_valpropre_deformation && !defPrincipales) valPropreEps = new Coordonnee(dim_sortie_tenseur); if (besoin_des_valpropre_vitdef && !vitPrincipales) valPropreDeps = new Coordonnee(dim_sortie_tenseur); if (besoin_des_valpropre_sigma && !sigmaPrincipales) valPropreSig = new Coordonnee(dim_sortie_tenseur); // on statut sur le fait d'avoir besoin ou non des chgt de base bool besoin_matrice_chg_base = false; if ( besoin_des_contraintes || besoin_des_deformation || besoin_des_vitesses_deformation || besoin_rep_local_ortho ) {besoin_matrice_chg_base = true;}; // definition générale def->ChangeNumInteg(iteg); // on change le numéro de point d'intégration courant if (((cas == 1) || (cas == 2))) { // cas d'une premiere initialisation Tableau tab(5); tab(1) = iM0; tab(2) = iMt; tab(3) = iMtdt; tab(4) = igiH_0; tab(5) = igijHH_0; met->PlusInitVariables(tab) ; }; // éléments de métrique et matrices de passage TenseurHH* gijHH=NULL;TenseurBB* gijBB=NULL;BaseB* giB=NULL; BaseH* giH_0=NULL;BaseH* giH=NULL; BaseB* giB_0=NULL; BaseB* giB_t=NULL; // n'est définit "que" pour certains cas Mat_pleine* Aa0 = NULL;Mat_pleine* Aafin = NULL; Mat_pleine* gamma0 = NULL;Mat_pleine* gammafin = NULL; Mat_pleine* beta0 = NULL;Mat_pleine* betafin = NULL; if (besoin_matrice_chg_base) // dans le cas où on n'est pas en absolue => on sort dans un repère ad hoc donc // il a la dimension locale // sinon on sort dans le repère globale => il a la dimension globale {int dim_effective = dim; // init if (absolue) dim_effective = ParaGlob::Dimension(); Aa0 = new Mat_pleine(dim_effective,dim_effective); Aafin = new Mat_pleine(dim_effective,dim_effective); gamma0 = new Mat_pleine(dim_effective,dim_effective); gammafin = new Mat_pleine(dim_effective,dim_effective); beta0 = new Mat_pleine(dim_effective,dim_effective); betafin = new Mat_pleine(dim_effective,dim_effective); }; switch (cas) { case 1: case 11: // calcul d'une ortho interessance de visualisation des tenseurs // cas de tenseur 3D -> Ia, cas 1D on prend un vecteur norme collineaire // avec g1, dans le cas 2D // la nouvelle base gamma est calculee dans def par projection de "Ipa" sur Galpha // le resultat est une matrice de passage utilisable pour le changement de base // Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a // tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer // resultat a t+dt {if (besoin_matrice_chg_base) {const Met_abstraite::InfoImp& ex = def->RemontImp(absolue,*Aa0,*Aafin); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; giB_0 = ex.giB_0; // a priori ici il s'agit d'un post-traitement donc la valeur à t n'est pas pertinente giB_t = giB; // donc on met la même valeur à t et tdt } else {const Met_abstraite::InfoImp& ex = def->RemontImp(); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; giB_0 = ex.giB_0; // a priori ici il s'agit d'un post-traitement donc la valeur à t n'est pas pertinente giB_t = giB; // donc on met la même valeur à t et tdt } break; } case -1: case -2: // identique au cas 1 mais avec le fait que la métrique est directement disponible // car le calcul est supposé suivre un calcul implicit au bon pti {if (besoin_matrice_chg_base) {Mat_pleine Aat(dim,dim); // a priori Aat ne sert pas par la suite, mais est nécessaire pour le passage de par const Met_abstraite::Info_et_metrique_0_t_tdt ex = def->Remont_et_metrique_0_t_tdtSansCalMet(absolue,*Aa0,Aat,*Aafin); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; giB_0 = ex.giB_0;giB_t = ex.giB_t; } else {const Met_abstraite::Info_et_metrique_0_t_tdt ex = def->Remont_et_metrique_0_t_tdtSansCalMet(); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; giB_0 = ex.giB_0;giB_t = ex.giB_t; } break; } case 2: case 12: // resultat a t {if (besoin_matrice_chg_base) {const Met_abstraite::InfoExp_tdt& ex= def->RemontExp_tdt(absolue,*Aa0,*Aafin); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; } else {const Met_abstraite::InfoExp_tdt& ex= def->RemontExp_tdt(); gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; } break; }; default: cout << "\n *** cas non prevu cas = "<< cas << "\n ElemMeca::Valeurs_Tensorielles(..." << endl; Sortie(1); }; if (besoin_matrice_chg_base) {// pour les formules de passage de repère il nous faut : // Ip_a = beta_a^{.j} g_j et Ip^b = gamma^b_{.j} g^j // on a: [beta_a^{.j}] = [Aa^j_{.a}]^T // et [gamma^b_{.j}] = [beta_a^{.j}]^{-1T} = [Aa^j_{.a}]^{-1} (*gamma0) = (Aa0->Inverse()); (*gammafin) = (Aafin->Inverse()); // on détermine également les matrices beta (*beta0) = (Aa0->Transpose()); (*betafin) = (Aafin->Transpose()); }; // récup des bases si besoin if (besoin_rep_local_ortho) {if ((!absolue)&&(dim_espace==3)&&(dim==2)) // cas d'éléments 2D, pour lesquels on veut un repère local ad hoc // on ramène une base à 3 vecteurs { *(base_ad_hoc(1)) = (*gammafin)(1,1) * giB->Coordo(1) + (*gammafin)(2,1) * giB->Coordo(2); *(base_ad_hoc(2)) = (*gammafin)(1,2) * giB->Coordo(1) + (*gammafin)(2,2) * giB->Coordo(2); *(base_ad_hoc(3)) = Util::ProdVec_coor(*(base_ad_hoc(1)),*(base_ad_hoc(2))); } else if((!absolue)&&(dim_espace>1)&&(dim==1)) // cas d'éléments 1D, dans un espace 2D ou 3D // on ramène un seul vecteur non nul, les autres ne peuvent être calculé sans info supplémentaire { *(base_ad_hoc(1)) = (*gammafin)(1,1) * giB->Coordo(1); (base_ad_hoc(2))->Zero(); (base_ad_hoc(3))->Zero(); // init à 0 par défaut } else // dans tous les autres cas { switch (dim_espace) // on se contente de ramener le repère identité {case 3: (base_ad_hoc(3))->Zero();(*(base_ad_hoc(3)))(3)=1.; case 2: (base_ad_hoc(2))->Zero();(*(base_ad_hoc(2)))(2)=1.; case 1: (base_ad_hoc(1))->Zero();(*(base_ad_hoc(1)))(1)=1.; default:break; }; }; }; if (besoin_rep_giH) {switch (dim) {case 3: *(base_giH(3)) = giH->Coordo(3); case 2: *(base_giH(2)) = giH->Coordo(2); case 1: *(base_giH(1)) = giH->Coordo(1); default:break; }; }; if (besoin_rep_giB) {switch (dim) {case 3: *(base_giB(3)) = giB->Coordo(3); case 2: *(base_giB(2)) = giB->Coordo(2); case 1: *(base_giB(1)) = giB->Coordo(1); default:break; }; }; // ----- calcul des grandeurs à sortir // calcul des tenseurs initiaux bool plusZero = true; // s'il faut rajouter des termes, on met des 0 if (besoin_des_contraintes) {(*sigHB) = *(ptIntegMeca.SigHH()) * (*gijBB); if (contrainteCourante) // on n'intervient ici que si on veut une sortie des contraintes courantes {if (absolue) {(ptIntegMeca.SigHH())->BaseAbsolue(*sigHH,*giB);}// changement de base finale else { *sigHH = *(ptIntegMeca.SigHH());sigHH->ChBase(*gammafin);}; }; }; if (besoin_des_deformation) {// cas de delta_eps if(DeltaEpsBB != NULL) {if (absolue)// changement de base finale {(ptIntegMeca.DeltaEpsBB())->BaseAbsolue(*DeltaEpsBB,*giH);} else {*DeltaEpsBB = *(ptIntegMeca.DeltaEpsBB());DeltaEpsBB->ChBase(*betafin);}; }; // cas de la déformation (*epsHB) = (*gijHH) * (*(ptIntegMeca.EpsBB())); if (deformationCourante)// on n'intervient ici que si on veut une sortie des déformations courantes {if (absolue)// changement de base finale {(ptIntegMeca.EpsBB())->BaseAbsolue(*epsBB,*giH);} else {*epsBB = *(ptIntegMeca.EpsBB());epsBB->ChBase(*betafin);}; }; switch (def->Type_de_deformation()) {case DEFORMATION_STANDART : // c'est à dire almansi {if (greenLagrange) {if (absolue) {(ptIntegMeca.EpsBB())->BaseAbsolue(*eps0BB,*giH_0);}// changement de base finale else {eps0BB->ChBase(*beta0);}; }; if (almansi) *epsAlmBB = *epsBB;// le changement de bas a été fait juste plus haut if (almansiTotal) // cas avec dilatation et demande de def Almansi totale {TenseurBB* epsAlmTotal_local_BB = epsAlmTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsAlmTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsAlmTotal_local_BB); if (absolue)// changement de base finale {epsAlmTotal_local_BB->BaseAbsolue(*epsAlmTotalBB,*giH);} else {epsAlmTotalBB->ChBase(*betafin);}; // ici epsAlmTotal_local_BB == epsAlmTotalBB if (prevoir_change_dim_tenseur) delete epsAlmTotal_local_BB; // car pas utilisé ensuite }; if (greenLagrangeTotale) // cas avec dilatation et demande de def Green_Lagrange totale {TenseurBB* epsGLTotal_local_BB = epsGLTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsGLTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsGLTotal_local_BB); if (absolue)// changement de base finale {epsGLTotal_local_BB->BaseAbsolue(*epsGLTotalBB,*giH_0);} else {epsGLTotalBB->ChBase(*beta0);}; // ici epsGLTotal_local_BB == epsGLTotalBB if (prevoir_change_dim_tenseur) delete epsGLTotal_local_BB; // car pas utilisé ensuite }; // si l'on veut sortir la déformation logarithmique le plus simple est de la calculer if (logarithmiqueTotale || logarithmique) {def->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE); if (logarithmique) // cas du calcul de la def logarithmique {TenseurBB* epslog_local_BB = epslogBB; // par défaut if (prevoir_change_dim_tenseur) epslog_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epslog_local_BB); if (absolue)// changement de base finale {epslog_local_BB->BaseAbsolue(*epslogBB,*giH);} else {epslogBB->ChBase(*betafin);}; // ici epslog_local_BB == epslogBB if (prevoir_change_dim_tenseur) delete epslog_local_BB; // car pas utilisé ensuite }; if (logarithmiqueTotale) // cas avec dilatation et demande de def log totale {TenseurBB* epsLogTotal_local_BB = epsLogTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsLogTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsLogTotal_local_BB); if (absolue)// changement de base finale {epsLogTotal_local_BB->BaseAbsolue(*epsLogTotalBB,*giH);} else {epsLogTotalBB->ChBase(*betafin);}; // ici epsLogTotal_local_BB == epsLogTotalBB if (prevoir_change_dim_tenseur) delete epsLogTotal_local_BB; // car pas utilisé ensuite }; def->Change_type_de_deformation(DEFORMATION_STANDART); // on revient au type initial }; break; } case DEFORMATION_LOGARITHMIQUE : { if (logarithmique) *epslogBB=*epsBB; if (logarithmiqueTotale) // cas avec dilatation et demande de def log totale {TenseurBB* epsLogTotal_local_BB = epsLogTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsLogTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsLogTotal_local_BB); if (absolue)// changement de base finale {epsLogTotal_local_BB->BaseAbsolue(*epsLogTotalBB,*giH);} else {epsLogTotalBB->ChBase(*betafin);}; // ici epsLogTotal_local_BB == epsLogTotalBB if (prevoir_change_dim_tenseur) delete epsLogTotal_local_BB; // car pas utilisé ensuite }; // si l'on veut sortir la déformation d'Almansi ou de green-lagrange le plus simple est de les calculer if (almansi || greenLagrangeTotale || almansiTotal) {def->Change_type_de_deformation(DEFORMATION_STANDART); if (almansi) // cas de la def d'almansi { TenseurBB* eps_local_BB = epsAlmBB; // par défaut if (prevoir_change_dim_tenseur) eps_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*eps_local_BB); if (absolue)// changement de base finale {eps_local_BB->BaseAbsolue(*epsAlmBB,*giH);} else {epsAlmBB->ChBase(*betafin);};// ici eps_local_BB == epsAlmBB if(greenLagrange) {if (absolue)// changement de base finale {eps_local_BB->BaseAbsolue(*eps0BB,*giH_0);} else {eps0BB->ChBase(*beta0);}; // ici eps_local_BB == eps0BB }; if (prevoir_change_dim_tenseur) delete eps_local_BB; // car pas utilisé ensuite }; if (almansiTotal) // cas avec dilatation et demande de def Almansi totale {TenseurBB* epsAlmTotal_local_BB = epsAlmTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsAlmTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsAlmTotal_local_BB); if (absolue)// changement de base finale {epsAlmTotal_local_BB->BaseAbsolue(*epsAlmTotalBB,*giH);} else {epsAlmTotalBB->ChBase(*betafin);}; // ici epsAlmTotal_local_BB == epsAlmTotalBB if (prevoir_change_dim_tenseur) delete epsAlmTotal_local_BB; // car pas utilisé ensuite }; if (greenLagrangeTotale) // cas avec dilatation et demande de def Green_Lagrange totale {TenseurBB* epsGLTotal_local_BB = epsGLTotalBB; // par défaut if (prevoir_change_dim_tenseur) epsGLTotal_local_BB = NevezTenseurBB(dim); def->Cal_deformation (temps,*epsGLTotal_local_BB); if (absolue)// changement de base finale {epsGLTotal_local_BB->BaseAbsolue(*epsGLTotalBB,*giH_0);} else {epsGLTotalBB->ChBase(*beta0);}; // ici epsGLTotal_local_BB == epsGLTotalBB if (prevoir_change_dim_tenseur) delete epsGLTotal_local_BB; // car pas utilisé ensuite }; def->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE); // on revient au type initial }; break; } default: cout << "\n cas de deformation non encore implante en sortie de visualisation " << Nom_type_deformation(def->Type_de_deformation()) << " affichage donc errone des valeurs !!!"; }; }; //-- fin de if (besoin_des_deformation) if (besoin_des_vitesses_deformation) { *DepsHB = (*gijHH) * (*(ptIntegMeca.DepsBB())); if (absolue)// changement de base finale {(ptIntegMeca.DepsBB())->BaseAbsolue(*DepsBB,*giH);} else {*DepsBB = *(ptIntegMeca.DepsBB());DepsBB->ChBase(*betafin);}; }; if (besoin_des_contraintes_barre) {double Isig = sigHB->Trace(); // trace de la déformation *sig_barreHB = (*sigHB) - (Isig/dim_espace) * (*Id_dim_HB(dim)); }; if (besoin_des_deformation_barre) {double Ieps = epsHB->Trace(); // trace de la déformation *eps_barreHB = (*epsHB) - (Ieps/dim_espace) * (*Id_dim_HB(dim)); }; if (besoin_des_vitesses_deformation_barre) {double IDeps = DepsHB->Trace(); // trace de la déformation *Deps_barreHB = (*DepsHB) - (IDeps/dim_espace) * (*Id_dim_HB(dim)); }; // cas des valeurs propres et éventuellement des vecteurs propres {int caas=0; ////----- debug //cout << "\n besoin_des_valpropre_sigma= "<< besoin_des_valpropre_sigma // << " besoin_dir_princ_sig= "<< besoin_dir_princ_sig << endl; ////--- fin debug if (besoin_des_valpropre_sigma && besoin_dir_princ_sig)// on veut les directions et les valeurs propres { Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH Coordonnee inter(sigHB->ValPropre(caas,mat)); inter.Change_dim(dim_sortie_tenseur); // on étant la taille éventuellement *valPropreSig = inter; // sauvegarde des valeurs propres // puis des directions principales, que l'on doit exprimer dans le repère absolue Tableau tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor // on regarde le cas particulier ou les valeurs propres sont toutes nulles // si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles if (inter.Max_val_abs() > ConstMath::petit) // et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite {if (caas != -1) {switch (dim) {case 3: giB->BaseAbsolue( *(base_propre_sigma(3)),tab(3));(base_propre_sigma(3))->Normer(); giB->BaseAbsolue( *(base_propre_sigma(2)),tab(2));(base_propre_sigma(2))->Normer(); case 1: giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer(); break; case 2: // en 2D le premier vecteur en giB et le second en giH !! giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer(); giH->BaseAbsolue( *(base_propre_sigma(2)),mat.CoordonneeB_Base_associee(2)); (base_propre_sigma(2))->Normer(); break; default:break; }; }; }; } else if (besoin_dir_princ_sig) // cas où on ne veut que les directions principales { Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH Coordonnee inter(sigHB->ValPropre(caas,mat)); // sauvegarde des directions principales Tableau tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor // on regarde le cas particulier ou les valeurs propres sont toutes nulles // si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles if (inter.Max_val_abs() > ConstMath::petit) // et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite {if (caas != -1) {switch (dim) {case 3: giB->BaseAbsolue( *(base_propre_sigma(3)),tab(3));(base_propre_sigma(3))->Normer(); giB->BaseAbsolue( *(base_propre_sigma(2)),tab(2));(base_propre_sigma(2))->Normer(); case 1: giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer(); break; case 2: // en 2D le premier vecteur en en giB et le second en giH !! giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer(); giH->BaseAbsolue( *(base_propre_sigma(2)),mat.CoordonneeB_Base_associee(2)); (base_propre_sigma(2))->Normer(); break; default:break; }; }; }; } else if (besoin_des_valpropre_sigma)// on ne veut que les valeurs propres { Coordonnee inter(sigHB->ValPropre(caas)); inter.Change_dim(dim_sortie_tenseur); // on étend la taille éventuellement *valPropreSig = inter; // sauvegarde des valeurs propres }; // gestion d'une erreur éventuelle if (caas == -1) { cout << "\n warning *** erreur dans le calcul des valeurs propres et-ou vecteurs propres de la contrainte "; if (ParaGlob::NiveauImpression() > 5) {sigHB->Ecriture(cout);cout << "\n cas = "<< cas ; cout << "\nElemMeca::Valeurs_Tensorielles(...";}; cout << endl; }; }; // -- fin de l'encapsulation du cas des contraintes {int caas=0; if (besoin_des_valpropre_deformation && besoin_dir_princ_eps)// on veut les directions et les valeurs propres { Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH // 1) calcul des valeurs propres dans un vecteur de travail inter, et des vecteurs propres stockés Coordonnee inter(epsHB->ValPropre(caas,mat)); // dans la matrice mat inter.Change_dim(dim_sortie_tenseur); // on étend la taille, dans le cas ou le nb de valeurs // propre est inférieur à la dimension de l'espace *valPropreEps = inter; // sauvegarde des valeurs propres // 2) puis des directions principales, que l'on doit exprimer dans le repère absolue Tableau tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor // on regarde le cas particulier ou les valeurs propres sont toutes nulles // si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles if (inter.Max_val_abs() > ConstMath::petit) // et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite {if (caas != -1) {switch (dim) {case 3: giB->BaseAbsolue( *(base_propre_eps(3)),tab(3));(base_propre_eps(3))->Normer(); giB->BaseAbsolue( *(base_propre_eps(2)),tab(2));(base_propre_eps(2))->Normer(); case 1: giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer(); break; case 2: // en 2D le premier vecteur en en giB et le second en giH !! giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer(); giH->BaseAbsolue( *(base_propre_eps(2)),mat.CoordonneeB_Base_associee(2)); (base_propre_eps(2))->Normer(); break; default:break; }; }; }; } else if (besoin_dir_princ_eps) // cas où on ne veut que les directions principales { Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH Coordonnee inter(epsHB->ValPropre(caas,mat)); // dans la matrice mat // puis des directions principales, que l'on doit exprimer dans le repère absolue Tableau tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor // on regarde le cas particulier ou les valeurs propres sont toutes nulles // si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles if (inter.Max_val_abs() > ConstMath::petit) // et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite {if (caas != -1) {switch (dim) {case 3: giB->BaseAbsolue( *(base_propre_eps(3)),tab(3));(base_propre_eps(3))->Normer(); giB->BaseAbsolue( *(base_propre_eps(2)),tab(2));(base_propre_eps(2))->Normer(); case 1: giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer(); break; case 2: // en 2D le premier vecteur en en giB et le second en giH !! giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer(); giH->BaseAbsolue( *(base_propre_eps(2)),mat.CoordonneeB_Base_associee(2)); (base_propre_eps(2))->Normer(); break; default:break; }; }; }; } else if (besoin_des_valpropre_deformation)// on ne veut que les valeurs propres { Coordonnee inter(epsHB->ValPropre(caas)); inter.Change_dim(dim_sortie_tenseur); // on étant la taille éventuellement *valPropreEps = inter; // sauvegarde des valeurs propres }; // gestion d'une erreur éventuelle if (caas == -1) { cout << "\n warning *** erreur dans le calcul des valeurs propres et-ou vecteurs propres de la deformation"; if (ParaGlob::NiveauImpression() >= 7) {epsHB->Ecriture(cout);cout << "\n cas = "<< cas ; cout << "\nElemMeca::Valeurs_Tensorielles(...";}; cout << endl; }; }; // -- fin de l'encapsulation du cas des déformations {int caas=0; // il faut l'initialiser sinon il peut prendre la valeur du cas précédant if (besoin_des_valpropre_vitdef && besoin_dir_princ_D)// on veut les directions et les valeurs propres { Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH // 1) calcul des valeurs propres dans un vecteur de travail inter, et des vecteurs propres stockés Coordonnee inter(DepsHB->ValPropre(caas,mat)); // dans la matrice mat inter.Change_dim(dim_sortie_tenseur); // on étend la taille, dans le cas ou le nb de valeurs // propre est inférieur à la dimension de l'espace *valPropreDeps = inter; // sauvegarde des valeurs propres // 2) puis des directions principales, que l'on doit exprimer dans le repère absolue Tableau tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor // on regarde le cas particulier ou les valeurs propres sont toutes nulles // si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles if (inter.Max_val_abs() > ConstMath::petit) // et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite {if (caas != -1) {switch (dim) {case 3: giB->BaseAbsolue( *(base_propre_D(3)),tab(3));(base_propre_D(3))->Normer(); giB->BaseAbsolue( *(base_propre_D(2)),tab(2));(base_propre_D(2))->Normer(); case 1: giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer(); break; case 2: // en 2D le premier vecteur en en giB et le second en giH !! giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer(); giH->BaseAbsolue( *(base_propre_D(2)),mat.CoordonneeB_Base_associee(2)); (base_propre_D(2))->Normer(); break; default:break; }; }; }; } else if (besoin_dir_princ_D) // cas où on ne veut que les directions principales { Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH Coordonnee inter(DepsHB->ValPropre(caas,mat)); // dans la matrice mat // puis des directions principales, que l'on doit exprimer dans le repère absolue Tableau tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor // on regarde le cas particulier ou les valeurs propres sont toutes nulles // si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles if (inter.Max_val_abs() > ConstMath::petit) // et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite {if (caas != -1) {switch (dim) {case 3: giB->BaseAbsolue( *(base_propre_D(3)),tab(3));(base_propre_D(3))->Normer(); giB->BaseAbsolue( *(base_propre_D(2)),tab(2));(base_propre_D(2))->Normer(); case 1: giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer(); break; case 2: // en 2D le premier vecteur en en giB et le second en giH !! giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer(); giH->BaseAbsolue( *(base_propre_D(2)),mat.CoordonneeB_Base_associee(2)); (base_propre_D(2))->Normer(); break; default:break; }; }; }; } else if (besoin_des_valpropre_vitdef)// on ne veut que les valeurs propres { Coordonnee inter(DepsHB->ValPropre(caas)); inter.Change_dim(dim_sortie_tenseur); // on étant la taille *valPropreDeps = inter; // sauvegarde des valeurs propres }; // gestion d'une erreur éventuelle if (caas == -1) { cout << "\n warning *** erreur dans le calcul des valeurs propres et-ou vecteurs propres de la vitesse de deformation"; if (ParaGlob::NiveauImpression() >= 7) {DepsHB->Ecriture(cout);cout << "\n cas = "<< cas ; cout << "\nElemMeca::Valeurs_Tensorielles(...";}; cout << endl; }; }; // -- fin de l'encapsulation du cas des vitesses de déformations if (besoin_coordonnees) (*Mtdt) = def->Position_tdt(); if (besoin_coordonnees_t) (*Mt) = def->Position_t(); if (besoin_coordonnees_t0) (*M0) = def->Position_0(); // def éventuelle du vecteur normal: ceci n'est correct qu'avec une métrique 2D if (N_surf != NULL) { // on vérifie que la métrique est correcte if (giB->NbVecteur() != 2) {if (ParaGlob::NiveauImpression() > 0) {cout << "\n *** attention il ne s'agit pas d'un element 2D," << " le vecteur normal ne sera pas disponible"; if (ParaGlob::NiveauImpression() > 2) cout << "\n element: " << Num_elt() << " pti "<< iteg << "\n ElemMeca::Valeurs_Tensorielles(... "; }; cout << endl; } else // sinon c'est ok {// la normale vaut le produit vectoriel des 2 premiers vecteurs (*N_surf) = Util::ProdVec_coorBN( (*giB)(1), (*giB)(2)); N_surf->Normer(); // que l'on norme }; // on n'arrête pas l'exécution, car on pourrait vouloir sortir les normales pour un ensemble // d'éléments contenant des volumes, des surfaces, des lignes: bon... il y aura quand même des // pb au niveau des iso par exemple, du au fait que l'on va faire des moyennes sur des éléments // de type différents (à moins de grouper par type du coup on n'aura pas le warning }; // idem à l'instant t if (N_surf_t != NULL) { // on vérifie que la métrique est correcte if (giB_t->NbVecteur() != 2) {if (ParaGlob::NiveauImpression() > 0) {cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D," << " le vecteur normal ne sera pas correctement calcule"; if (ParaGlob::NiveauImpression() > 2) cout << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(... "; }; cout << endl; } else // sinon c'est ok {// la normale vaut le produit vectoriel des 2 premiers vecteurs (*N_surf_t) = Util::ProdVec_coorBN( (*giB_t)(1), (*giB_t)(2)); N_surf_t->Normer(); // que l'on norme }; }; // idem à l'instant t0 if (N_surf_t0 != NULL) { // on vérifie que la métrique est correcte if (giB_0->NbVecteur() != 2) {if (ParaGlob::NiveauImpression() > 0) {cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D," << " le vecteur normal ne sera pas correctement calcule"; if (ParaGlob::NiveauImpression() > 2) cout << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(... "; }; cout << endl; } else // sinon c'est ok {// la normale vaut le produit vectoriel des 2 premiers vecteurs (*N_surf_t0) = Util::ProdVec_coorBN( (*giB_0)(1), (*giB_0)(2)); N_surf_t0->Normer(); // que l'on norme }; }; // cas du déplacement et de la vitesse if (Deplacement != NULL) (*Deplacement) = (*Mtdt) - def->Position_0(); if (Vitesse != NULL) (*Vitesse) = def->VitesseM_tdt(); // def éventuelle de la contrainte de Mises if (contrainteMises) { Coordonnee& vv = *valPropreSig; int dimvec=vv.Dimension();// pour condenser l'écriture switch (dimvec) // dans le cas où dimvec=0 on ne fait rien, cas ou on n'a pas besoin de mises { case 1: *Mises = Dabs(vv(1)); break; case 2: *Mises = sqrt( ((vv(1)-vv(2))*(vv(1)-vv(2)) + vv(1) * vv(1) + vv(2) * vv(2)) * 0.5); break; case 3: *Mises = sqrt( ((vv(1)-vv(2))*(vv(1)-vv(2)) + (vv(1)-vv(3))*(vv(1)-vv(3)) + (vv(3)-vv(2))*(vv(3)-vv(2))) * 0.5); break; }; }; // def éventuelle de la déformation duale de mises = sqrt(2/3 * epsB:epsB) //--- on utilise directement la grandeur stockée au pt d'integ, mais ici ne sert à rien, puis cette grandeur n'est pas utilisée par la suite ! // // l'expression est la même que celle de mises, ormis un coeff 4/9 qui permet de passer de 3/2 à 2/3 // if (def_duale_mises) // { Coordonnee& vdef = *valPropreEps; int dimvdef=vdef.Dimension();// pour condenser l'écriture // switch (dimvdef) // dans le cas où dimvdef=0 on ne fait rien, cas ou on n'a pas besoin de defDualMises // { case 1: *defDualMises = Dabs(vdef(1)); break; // case 2: *defDualMises = sqrt( ((vdef(1)-vdef(2))*(vdef(1)-vdef(2)) + vdef(1) * vdef(1) // + vdef(2) * vdef(2)) * 2./9.); break; // case 3: *defDualMises = sqrt( ((vdef(1)-vdef(2))*(vdef(1)-vdef(2)) + (vdef(1)-vdef(3))*(vdef(1)-vdef(3)) // + (vdef(3)-vdef(2))*(vdef(3)-vdef(2))) * 2./9.); break; // }; // }; // donnees propre a la loi mécanique au pt d'integ Loi_comp_abstraite::SaveResul * sDon = tabSaveDon(iteg); // donnees propre a la loi thermo physique au pt d'integ CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques if (loiTP != NULL) {sTP = tabSaveTP(iteg);}; // au pt d'integ si elles existes // donnees propre a la déformation mécanique au pt d'integ Deformation::SaveDefResul * sDefDon = tabSaveDefDon(iteg); //contrainte_tresca if (contraintesTresca) { switch (dim) {case 1: *Tresca=0.5 * (*valPropreSig)(1);break; case 2: *Tresca=0.5 * ((*valPropreSig)(1)-(*valPropreSig)(2));break; case 3: *Tresca=0.5 * ((*valPropreSig)(1)-(*valPropreSig)(3));break; }; }; // cas de l'erreur if (erreur != NULL) {if (sigErreur != NULL) *erreur= *sigErreur; else {if (ParaGlob::NiveauImpression() > 4) cout << "\n pas encore de ddl erreur dans l'element " << "\n ElemMeca::Valeurs_Tensorielles(...."; }; }; // cas de l'erreur relative sur les contraintes if (erreur_rel != NULL) {if (sigErreur_relative != NULL) *erreur_rel= *sigErreur_relative; else {if (ParaGlob::NiveauImpression() > 4) cout << "\n pas encore de ddl erreur relative dans l'element " << "\n ElemMeca::Valeurs_Tensorielles(...."; }; }; // --- cas des grandeurs de la décomposition polaire // cas de la déformation if (def_SPHERIQUE_EPS) {*spherique_eps = epsHB->Trace()/3.;}; //ParaGlob::Dimension();}; modif 5/2/2012 double mini_Q = 5.e-5; if (def_Q_EPS) {*Q_eps = sqrt(eps_barreHB->II());}; if (cos3phi_eps) { double Qepsilon = ( def_Q_EPS ? *Q_eps : sqrt(eps_barreHB->II())); double Qepsilon3 = Qepsilon * Qepsilon * Qepsilon; if (Qepsilon > mini_Q ) { // on peut calculer un cos3phi pas débile double bIIIb = eps_barreHB->III() / 3.; *cos3phi_eps = 3. * sqrt(6.) * bIIIb/ Qepsilon3; } else *cos3phi_eps=0.; // sinon on le met à 0 }; // cas de la contrainte if (def_SPHERIQUE_SIG) {*spherique_sig = sigHB->Trace()/3.;}; //ParaGlob::Dimension();}; modif 5/2/2012 if (def_Q_SIG) {*Q_sig = sqrt(sig_barreHB->II());}; if (cos3phi_sig) { double Qsig = ( def_Q_SIG ? *Q_sig : sqrt(sig_barreHB->II())); double Qsig3 = Qsig * Qsig * Qsig; if (Qsig > mini_Q ) { // on peut calculer un cos3phi pas débile double bIIIb = sig_barreHB->III() / 3.; *cos3phi_sig = 3. * sqrt(6.) * bIIIb/ Qsig3; } else *cos3phi_sig=0.; // sinon on le met à 0 }; // cas de la vitesse de déformation if (def_SPHERIQUE_DEPS) {*spherique_Deps = DepsHB->Trace()/3.;}; //ParaGlob::Dimension();}; modif 5/2/2012 if (def_Q_DEPS) {*Q_Deps = sqrt(Deps_barreHB->II());}; if (cos3phi_Deps) { double QDepsilon = ( def_Q_DEPS ? *Q_Deps : sqrt(Deps_barreHB->II())); double QDepsilon3 = QDepsilon * QDepsilon * QDepsilon; if (QDepsilon > mini_Q ) { // on peut calculer un cos3phi pas débile double bIIIb = Deps_barreHB->III() / 3.; *cos3phi_Deps = 3. * sqrt(6.) * bIIIb/ QDepsilon3; } else *cos3phi_Deps=0.; // sinon on le met à 0 }; // cas de la deformation équivalente cumulée if (def_def_equivalente) {*def_equivalente = ptIntegMeca.Deformation_equi_const()(1);}; // cas de la deformation duale au sens de mises, if (def_duale_mises) {*defDualMises = ptIntegMeca.Deformation_equi_const()(2);}; // cas de la deformation maxi duale au sens de mises, if (def_duale_mises_maxi) {*defDualMisesMaxi = ptIntegMeca.Deformation_equi_const()(3);}; // delete des tenseurs if (sigHB != NULL) delete sigHB; if (epsHB != NULL) delete epsHB; if (DepsHB != NULL) delete DepsHB; if (eps_barreHB != NULL) delete eps_barreHB; if (Deps_barreHB != NULL) delete Deps_barreHB; if (sig_barreHB != NULL) delete sig_barreHB; // effacement conditionnel if (besoin_des_contraintes && !contrainteCourante) delete sigHH; if (besoin_des_deformation && !deformationCourante) delete epsBB; if (vitesseDeformationCourante && !vitesseDeformationCourante) delete DepsBB; if (besoin_des_valpropre_deformation && !defPrincipales) delete valPropreEps; if (besoin_des_valpropre_vitdef && !vitPrincipales) delete valPropreDeps ; if (besoin_des_valpropre_sigma && !sigmaPrincipales) delete valPropreSig; // pointeurs de matrice if (Aa0 != NULL) delete Aa0; if (Aafin != NULL) delete Aafin; if (gamma0 != NULL) delete gamma0; if (gammafin != NULL) delete gammafin; if (beta0 != NULL) delete beta0; if (betafin != NULL) delete betafin; // liberation des tenseurs intermediaires LibereTenseur(); def->Retour_pti_precedant(); // on revient au pti précédant };