// 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 "CharUtil.h" #include "TypeQuelconqueParticulier.h" // actualisation des ddl et des grandeurs actives de t+dt vers t // appelé par les classes dérivées void ElemMeca::TdtversT_() {// on met à jour l'indicateur de premier calcul // s'il y a des sauvegardes de grandeur aux déformations // on ne regarde que le premier élément de tableau, a priori // il y a toujours un pt d'integ et l'organisation est identique pour tous les pt d'integ if (tabSaveDefDon(1) != NULL) premier_calcul_meca_impli_expli=false; // cas des énergies int nbi= tab_energ.Taille(); for (int ni=1;ni<= nbi; ni++) tab_energ_t(ni) = tab_energ(ni); E_elem_bulk_t = E_elem_bulk_tdt; // énergie due au bulk viscosity energie_totale_t = energie_totale; // les énergies globalisées sur l'élément // cas des intégrales de volumes et temps éventuelles if (integ_vol_typeQuel != NULL) {int taille = integ_vol_typeQuel->Taille(); Tableau & tab_integ = *integ_vol_typeQuel; // pour simplifier Tableau & tab_integ_t = *integ_vol_typeQuel_t; // "" for (int il =1;il <= taille ;il++) {tab_integ_t(il) = tab_integ(il); // tab_integ(il).Grandeur_pointee()->InitParDefaut(); // on initialise pour le prochain calcul }; }; if (integ_vol_t_typeQuel != NULL) {int taille = integ_vol_t_typeQuel->Taille(); Tableau & tab_integ = *integ_vol_t_typeQuel; // pour simplifier Tableau & tab_integ_t = *integ_vol_t_typeQuel_t; // "" for (int il =1;il <= taille ;il++) {tab_integ_t(il) = tab_integ(il); }; }; // force de stabilisation éventuelle if (pt_StabMembBiel != NULL) pt_StabMembBiel->TdtversT(); }; // actualisation des ddl et des grandeurs actives de t vers tdt // appelé par les classes dérivées void ElemMeca::TversTdt_() {// on met à jour l'indicateur de premier calcul // on considère que si l'on revient en arrière, il vaut mieux re-initialiser les // grandeurs correspondantes au premier_calcul // s'il y a des sauvegardes de grandeur aux déformations if (tabSaveDefDon(1) != NULL) premier_calcul_meca_impli_expli=true; // cas des énergies int nbi= tab_energ.Taille(); for (int ni=1;ni<= nbi; ni++) tab_energ(ni) = tab_energ_t(ni); E_elem_bulk_tdt = E_elem_bulk_t; // énergie due au bulk viscosity energie_totale = energie_totale_t; // les énergies globalisées sur l'élément // cas des intégrales de volumes et temps éventuelles // il n'y a rien à faire // if (integ_vol_typeQuel != NULL) // {int taille = integ_vol_typeQuel->Taille(); // Tableau & tab_integ = *integ_vol_typeQuel; // pour simplifier // Tableau & tab_integ_t = *integ_vol_typeQuel_t; // "" // for (int il =1;il <= taille ;il++) // tab_integ(il) = tab_integ_t(il); //// tab_integ(il).Grandeur_pointee()->InitParDefaut(); // on initialise pour le prochain calcul // }; if (integ_vol_t_typeQuel != NULL) {int taille = integ_vol_t_typeQuel->Taille(); Tableau & tab_integ = *integ_vol_t_typeQuel; // pour simplifier Tableau & tab_integ_t = *integ_vol_t_typeQuel_t; // "" for (int il =1;il <= taille ;il++) {tab_integ(il) = tab_integ_t(il); // on met à jour le cumul }; }; // force de stabilisation éventuelle if (pt_StabMembBiel != NULL) pt_StabMembBiel->TdtversT(); }; // calcul du résidu et de la matrice de raideur pour le calcul d'erreur // cas d'une intégration suivant une seule liste void ElemMeca::SigmaAuNoeud_ResRaid(const int nbne,const Tableau & taphi ,const Vecteur& poids,Tableau & resErr,Mat_pleine& raidErr ,const Tableau & taphiEr,const Vecteur& poidsEr) { int dimAbsolue = ParaGlob::Dimension(); // dimension absolue // dimension des tenseurs // int dim = lesPtIntegMecaInterne->DimTens(); // inialisation du second membre et de la raideur int nbSM = resErr.Taille(); // nombre de second membre for (int j =1; j<= nbSM; j++) // boucle sur les seconds membres (*resErr(j)).Zero(); raidErr.Zero(); // Il faut déterminer l'ordre dans lequel on parcours les contraintes qui doit // être compatible avec l'ordre des ddl Tableau2 ordre = OrdreContrainte(nbSM); // création d'un tenseur au dimension absolu pour le calcul des contraintes // dans la base absolue TenseurHH& sigHH = *(NevezTenseurHH(dimAbsolue)) ; // ====== calcul du second membre ======= int ni; // compteur globale de point d'integration bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde // donc il faut calculer tous les éléments de la métrique for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++) {PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni); // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul); // passage dans le repère absolu du tenseur contrainte final sigHH = (*(ptIntegMeca.SigHH_t())).BaseAbsolue(sigHH,*ex.giB_t); for (int ne =1; ne<= nbne; ne++) // 1ere boucle sur les noeuds { // les résidus : mais il faut suivre l'ordre de l'enregistrement des ddl for (int itot = 1; itot<= nbSM; itot++) { int ix = (int) (ordre(itot,1)); int iy = (int) (ordre(itot,2)); (*resErr(itot))(ne) += taphi(ni)(ne)*sigHH(ix,iy) * (poids(ni) * (*ex.jacobien_t)); }; }; }; // ====== calcul de la raideur c'est à dire du hessien ======== // boucle sur les pt d'integ spécifiques à l'erreur for (defEr->PremierPtInteg(), ni = 1;defEr->DernierPtInteg();defEr->NevezPtInteg(),ni++) { // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique const Met_abstraite::Expli& ex = defEr->Cal_explicit_t(premier_calcul); for (int ne =1; ne<= nbne; ne++) // 1ere boucle sur les noeuds for (int me =1; me<= nbne; me++) // 2ere boucle sur les noeuds raidErr(ne,me) += taphiEr(ni)(ne) * taphiEr(ni)(me) * poidsEr(ni) * (*ex.jacobien_t); }; // liberation des tenseurs intermediaires TenseurHH * ptsigHH = &sigHH; delete ptsigHH; LibereTenseur(); // calcul de l'erreur relative }; // calcul de l'erreur sur l'élément. Ce calcul n'est disponible // qu'une fois la remontée aux contraintes effectuées sinon aucune // action. En retour la valeur de l'erreur sur l'élément // = 1 : erreur = (int (delta sigma):(delta sigma) dv)/(int sigma:sigma dv) // le numerateur et le denominateur sont tel que : // errElemRelative = numerateur / denominateur , si denominateur different de 0 // sinon denominateur = numerateur si numerateur est different de 0, sinon // tous sont nuls mais on n'effectue pas la division , les autres variables sont spécifiques // a l'element. void ElemMeca::Cal_ErrElem(int type,double& errElemRelative,double& numerateur , double& denominateur,const int nbne,const Tableau & taphi ,const Vecteur& poids,const Tableau & taphiEr,const Vecteur& poidsEr) { int dimAbsolue = ParaGlob::Dimension(); // dimension absolue // dimension des tenseurs // int dim = lesPtIntegMecaInterne->DimTens(); bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde // donc il faut calculer tous les éléments de la métrique double domaine_integration=0.; // volume ou surface on longueur switch (type) {case 1 : // cas du calcul aux moindres carrés {// création d'un tenseur au dimension absolu pour le calcul des contraintes // dans la base absolue, on le choisit HB pour le double produit contracté // mais en absolu la variance n'a pas d'importance TenseurHB& sigHB = *(NevezTenseurHB(dimAbsolue)) ; // idem au point d'intégration et un tenseur nul pour l'initialisation TenseurHB& signiHB = *(NevezTenseurHB(dimAbsolue)) ; TenseurHB& sigiHB = *(NevezTenseurHB(dimAbsolue)) ; // tenseur de travail TenseurHB& nulHB = *(NevezTenseurHB(dimAbsolue)) ; // ====== calcul des termes de l'erreur ======= // ---- tout d'abord on parcourt les points d'intégration de la mécanique numerateur = 0.; // initialisation double numerateur_bis = 0.; denominateur = 0.;// intégrale double de la contrainte élément fini int ni; // compteur globale de point d'integration for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++) {PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni); // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul); // passage dans le repère absolu du tenseur contrainte initiale // car les contraintes aux noeuds on été calculées et stockées en absolue sigHB = (*(ptIntegMeca.SigHH_t())).BaseAbsolue(sigHB,*ex.giB_t); ////---debug //cout << "\n debug ElemMeca::Cal_ErrElem: sigHB= "; // sigHB.Ecriture(cout); ////---fin debug // calcul du denominateur double integ_sigHB_carre_du_pti = (sigHB && sigHB) * (poids(ni) * (*ex.jacobien_t)); denominateur += integ_sigHB_carre_du_pti; // calcul de la premiere partie du numerateur, celle qui dépend des points d'intégration // mécanique. // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds, signiHB = nulHB; // initialisation for (int ne =1; ne<= nbne; ne++) // boucle sur les noeuds signiHB += taphi(ni)(ne) * ((tab_noeud(ne))->Contrainte(sigiHB)); // 2) intégrale de la partie dépendante de ni numerateur += integ_sigHB_carre_du_pti - ((signiHB && sigHB)+(sigHB && signiHB)) * (poids(ni) * (*ex.jacobien_t)) ; numerateur_bis += (signiHB-sigHB)&&(signiHB-sigHB)* (poids(ni) * (*ex.jacobien_t)); }; // ---- on parcourt maintenant les points d'intégration pour le calcul d'erreur // ---- ce qui permet de calculer la deuxième partie du numérateur // boucle sur les pt d'integ spécifiques à l'erreur double denominateur_continu = 0.;// intégrale double de la contrainte continue domaine_integration=0.; // init for (defEr->PremierPtInteg(), ni = 1;defEr->DernierPtInteg();defEr->NevezPtInteg(),ni++) { // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique const Met_abstraite::Expli& ex = defEr->Cal_explicit_t(premier_calcul); // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds, signiHB = nulHB; // initialisation for (int ne =1; ne<= nbne; ne++) // boucle sur les noeuds signiHB += taphiEr(ni)(ne)*(tab_noeud(ne))->Contrainte(sigiHB); // 2) intégrale de la partie dépendant de ni ////---debug //cout << "\n debug ElemMeca::Cal_ErrElem: signiHB= "; // signiHB.Ecriture(cout); //cout << "\n debug ElemMeca::Cal_ErrElem: (tab_noeud(ne))->Contrainte(sigiHB)= "; //for (int ne =1; ne<= nbne; ne++) // boucle sur les noeuds // {cout <<"\n";(tab_noeud(ne))->Contrainte(sigiHB).Ecriture(cout);} // ////---fin debug double integ_signiHB_carre_du_pti = (signiHB && signiHB) * poidsEr(ni) * (*ex.jacobien_t); denominateur_continu += integ_signiHB_carre_du_pti; numerateur += integ_signiHB_carre_du_pti; // 3) intégrale du domaine d'intégration domaine_integration += poidsEr(ni) * (*ex.jacobien_t); }; // liberation des tenseurs intermediaires TenseurHB * ptsigHB = &sigHB; delete ptsigHB; ptsigHB = &signiHB; delete ptsigHB;ptsigHB = &sigiHB; delete ptsigHB; ptsigHB = &nulHB; delete ptsigHB; LibereTenseur(); // enregistrement de l'erreur pour l'élément if (sigErreur == NULL) {sigErreur = new double;sigErreur_relative = new double;}; // comme le calcul précédent correspond à l'intégrale d'un carré, // donc c'est analogue à "une grandeur moyenne au carré " * le volume //on prend la racine carrée de: // (" grandeur moyenne au carré " * le volume) / le volume c-a-d |grandeur moyenne| // qui est donc pondèré par le domaine ce qui nous donne une grandeur analogue à // un delta contrainte // on sauvegarde avant double numerateur_initale = numerateur; double denominateur_initiale = denominateur; *sigErreur = sqrt(Dabs(numerateur/domaine_integration)); // pour calculer une erreur relative on commence par calculer une grandeur moyenne // qui sera ici analogue à une contrainte denominateur = sqrt(Dabs(denominateur_initiale)/domaine_integration); // on essaie de gérer le cas des contraintes petites if (denominateur < ConstMath::unpeupetit) // on passe en erreur absolue si le dénominateur est trop petit { *sigErreur_relative = *sigErreur; } else {*sigErreur_relative = *sigErreur/(denominateur);}; numerateur = *sigErreur; //////------ debug //if (numerateur > 90) // {cout << "\n debug ElemMeca::Cal_ErrElem: *sigErreur= "<< *sigErreur // << " *sigErreur_relative = " << *sigErreur_relative // << "\n num élément: " << num_elt // << flush; // // création d'un tenseur au dimension absolu pour le calcul des contraintes // // dans la base absolue, on le choisit HB pour le double produit contracté // // mais en absolu la variance n'a pas d'importance // TenseurHB& sigHB = *(NevezTenseurHB(dimAbsolue)) ; // TenseurHH& sigHH = *(NevezTenseurHH(dimAbsolue)) ; // pour voir si identique à sigHB // // idem au point d'intégration et un tenseur nul pour l'initialisation // TenseurHB& signiHB = *(NevezTenseurHB(dimAbsolue)) ; // TenseurHB& sigiHB = *(NevezTenseurHB(dimAbsolue)) ; // tenseur de travail // TenseurHB& nulHB = *(NevezTenseurHB(dimAbsolue)) ; // // ====== calcul des termes de l'erreur ======= // // ---- tout d'abord on parcourt les points d'intégration de la mécanique // double numerateur = 0.; // initialisation // double numerateur_bis = 0.; // double denominateur = 0.;// intégrale double de la contrainte élément fini // int ni; // compteur globale de point d'integration // for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++) // {PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni); // // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul // // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique // const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul); // // passage dans le repère absolu du tenseur contrainte initiale // // car les contraintes aux noeuds on été calculées et stockées en absolue // sigHB = (*(ptIntegMeca.SigHH_t())).BaseAbsolue(sigHB,*ex.giB_t); // (ptIntegMeca.SigHH_t())->BaseAbsolue(sigHH,*ex.giB_t); ////---debug //cout << "\n debug ElemMeca::Cal_ErrElem: num_pti: " << ni // << "\n contrainte initiale sigHB= "; // sigHB.Ecriture(cout); //// cout << "\n et la version en HH : sigHH= "; ////sigHH.Ecriture(cout); ////---fin debug // // calcul du denominateur // double integ_sigHB_carre_du_pti = (sigHB && sigHB) * (poids(ni) * (*ex.jacobien_t)); // denominateur += integ_sigHB_carre_du_pti; // // // calcul de la premiere partie du numerateur, celle qui dépend des points d'intégration // // mécanique. // // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds, // signiHB = nulHB; // initialisation // for (int ne =1; ne<= nbne; ne++) // boucle sur les noeuds // signiHB += taphi(ni)(ne) * ((tab_noeud(ne))->Contrainte(sigiHB)); // // 2) intégrale de la partie dépendante de ni // numerateur += integ_sigHB_carre_du_pti - ((signiHB && sigHB)+(sigHB && signiHB)) * (poids(ni) * (*ex.jacobien_t)) ; // numerateur_bis += (signiHB-sigHB)&&(signiHB-sigHB)* (poids(ni) * (*ex.jacobien_t)); // }; // // ---- on parcourt maintenant les points d'intégration pour le calcul d'erreur // // ---- ce qui permet de calculer la deuxième partie du numérateur // // boucle sur les pt d'integ spécifiques à l'erreur // double denominateur_continu = 0.;// intégrale double de la contrainte continue // domaine_integration=0.; // init // for (defEr->PremierPtInteg(), ni = 1;defEr->DernierPtInteg();defEr->NevezPtInteg(),ni++) // { // // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul // // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique // const Met_abstraite::Expli& ex = defEr->Cal_explicit_t(premier_calcul); // // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds, // signiHB = nulHB; // initialisation // for (int ne =1; ne<= nbne; ne++) // boucle sur les noeuds // signiHB += taphiEr(ni)(ne)*(tab_noeud(ne))->Contrainte(sigiHB); // // 2) intégrale de la partie dépendant de ni ////---debug //cout << "\n\n debug ElemMeca::Cal_ErrElem: contraintes aux noeuds ramené au pti: "<Contrainte(sigiHB)= "; //for (int ne =1; ne<= nbne; ne++) // boucle sur les noeuds // {cout <<"\n";(tab_noeud(ne))->Contrainte(sigiHB).Ecriture(cout);} // ////---fin debug // double integ_signiHB_carre_du_pti = (signiHB && signiHB) * poidsEr(ni) * (*ex.jacobien_t); // denominateur_continu += integ_signiHB_carre_du_pti; // numerateur += integ_signiHB_carre_du_pti; // // 3) intégrale du domaine d'intégration // domaine_integration += poidsEr(ni) * (*ex.jacobien_t); // }; //cout << "\n debug ElemMeca::Cal_ErrElem: domaine_integration= " // << domaine_integration << flush; // // liberation des tenseurs intermediaires // TenseurHB * ptsigHB = &sigHB; delete ptsigHB; // ptsigHB = &signiHB; delete ptsigHB;ptsigHB = &sigiHB; delete ptsigHB; // ptsigHB = &nulHB; delete ptsigHB; // TenseurHH * ptsigHH = &sigHH; delete ptsigHH; // LibereTenseur(); // } ////------ fin debug // *sigErreur = Dabs(numerateur); // // numerateur = sqrt(Dabs(numerateur/domaine_integration)); // *sigErreur_relative = numerateur; //// numerateur = sqrt(Dabs(numerateur)/domaine_integration)*domaine_integration; // *sigErreur = numerateur; ////---essai sur des grandeurs relatives, mais en fait ce n'est pas une bonne idée, car ici c'est // la différence totale qui nous intéresse, et en relativ cela donne n'importe quoi // *sigErreur = numerateur; {// on commence par calculer le maxi des contraintes moyennes au carré, soient via le calcul élément fini initial // soit via la contrainte continu, et on évite d'avoir 0 //// double sig_moy_2 = DabsMaX(denominateur_continu,denominateur,ConstMath::trespetit); // double sig_moy_2 = DabsMaX(denominateur,ConstMath::trespetit); //////---debug //cout << "\n debug ElemMeca::Cal_ErrElem: " // << "\n Int_sig2_continu=" << denominateur_continu <<", Int_sig2_EF=" << denominateur // <<", rapport="< erreur relative = 1 // errElemRelative = 1./domaine_integration; // } // else // // cas du dénominateur non nul // { errElemRelative = numerateur / denominateur;}//(denominateur*domaine_integration);}; errElemRelative = *sigErreur_relative; }; // calcul de l'erreur aux noeuds. Contrairement au cas des contraintes // seul le résidu est calculé. Cas d'une intégration suivant une liste void ElemMeca::Cal_ErrAuxNoeuds(const int nbne, const Tableau & taphi, const Vecteur& poids,Tableau & resErr ) { //int dimAbsolue = ParaGlob::Dimension(); // dimension absolue // inialisation du second membre, on utilise le premier vecteur uniquement (*resErr(1)).Zero(); // ====== calcul du second membre ======= int ni; // compteur globale de point d'integration double volume = 0.; // le volume de l'élément bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde // donc il faut calculer tous les éléments de la métrique for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++) { // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul); for (int ne =1; ne<= nbne; ne++) // 1ere boucle sur les noeuds // (*resErr(1))(ne) += taphi(ni)(ne)*(*sigErreur) * (poids(ni) * (*ex.jacobien_t)); (*resErr(1))(ne) += taphi(ni)(ne)*(*sigErreur_relative) * (poids(ni) * (*ex.jacobien_t)); // calcul du volume volume += (poids(ni) * (*ex.jacobien_t)); } // on relativise par rapport au volume, car initialement sigErreur représente l'erreur // totale sur tout l'élément. Donc on divise par le volume pour retrouver après résolution // une valeur au noeud qui représente une valeur ponctuelle et non une valeur qui // qui est relative à un volume *resErr(1) /= volume; // ensuite on prend la racine carrée pour que ce soit homogène à une contrainte int tail = (*resErr(1)).Taille(); Vecteur& resErr_1 = *resErr(1); // pour simplifier for (int i=1;i<= tail;i++) {if (resErr_1(1) > 0.) resErr_1(1) = sqrt(resErr_1(1)); else { if (ParaGlob::NiveauImpression() > 3) cout << "\n *** warning : bizarre la composante "< 5) cout << "\n ElemMeca::Cal_ErrAuxNoeuds "; resErr_1(1) = sqrt(-resErr_1(1)); }; }; ////---debug // cout << "\n debug ElemMeca::Cal_ErrAuxNoeuds "; // resErr_1.Affiche(); // cout << endl; ////---fin debug // *resErr(1) = -sigErreur; LibereTenseur(); // calcul de l'erreur relative }; // ajout des ddl relatif aux contraintes pour les noeuds de l'élément void ElemMeca::Ad_ddl_Sigma(const DdlElement& tab_ddlErr) { int nbne = tab_noeud.Taille(); for (int ne=1; ne<= nbne; ne++) // pour chaque noeud { // création du tableau de ddl int tab_ddlErr_Taille = tab_ddlErr(ne).tb.Taille(); // nb de ddl du noeud meca Tableau ta(tab_ddlErr_Taille); for (int i =1; i<= tab_ddlErr_Taille; i++) ta(i) = Ddl (tab_ddlErr(ne).tb(i),0.,LIBRE); // ajout du tableau dans le noeud tab_noeud(ne)->PlusTabDdl(ta); } }; // inactive les ddl du problème primaire de mécanique void ElemMeca::Inact_ddl_primaire(DdlElement& tab_ddl) { int nbne = tab_noeud.Taille(); for (int ne=1; ne<= nbne; ne++) tab_noeud(ne)->Met_hors_service(tab_ddl(ne).tb); }; // active les ddl du problème primaire de mécanique void ElemMeca::Act_ddl_primaire(DdlElement& tab_ddl) { int nbne = tab_noeud.Taille(); for (int ne=1; ne<= nbne; ne++) tab_noeud(ne)->Met_en_service(tab_ddl(ne).tb); }; // inactive les ddl du problème de recherche d'erreur : les contraintes void ElemMeca::Inact_ddl_Sigma(DdlElement& tab_ddlErr) { int nbne = tab_noeud.Taille(); for (int ne=1; ne<= nbne; ne++) tab_noeud(ne)->Met_hors_service(tab_ddlErr(ne).tb); }; // active les ddl du problème de recherche d'erreur : les contraintes void ElemMeca::Act_ddl_Sigma(DdlElement& tab_ddlErr) { int nbne = tab_noeud.Taille(); for (int ne=1; ne<= nbne; ne++) tab_noeud(ne)->Met_en_service(tab_ddlErr(ne).tb); }; // active le premier ddl du problème de recherche d'erreur : SIGMA11 void ElemMeca::Act_premier_ddl_Sigma() { int nbne = tab_noeud.Taille(); for (int ne=1; ne<= nbne; ne++) tab_noeud(ne)->Met_en_service(SIG11); }; // retourne un tableau de ddl element, correspondant à la // composante de sigma -> SIG11, pour chaque noeud qui contiend // des ddl de contrainte // -> utilisé pour l'assemblage de la raideur d'erreur DdlElement& ElemMeca::Tableau_de_Sig1() const { cout << "\n erreur, fonction non defini pour cette element " << "\n ElemMeca::Tableau_de_Sig1()" << endl; Sortie(1); DdlElement * toto = new DdlElement(); return *toto; }; // lecture des contraintes sur le flot d'entrée void ElemMeca::LectureDesContraintes(bool cas,UtilLecture * entreePrinc,Tableau & tabSigHH) { // dimensionnement de la metrique identique au cas d'un calcul explicite if( cas) { Tableau tab(7); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t; tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ; met->PlusInitVariables(tab) ; }; int dimAbsolue = ParaGlob::Dimension(); // dimension absolue // dimension des tenseurs //int dim = (*(tabSigHH(1))).Dimension(); // récup du nombre de composantes du tenseur int nbcomposante= ParaGlob::NbCompTens(); // Il faut déterminer l'ordre dans lequel on parcours les contraintes qui doit // être compatible avec l'ordre des ddl Tableau2 ordre = OrdreContrainte(nbcomposante); // création d'un tenseur au dimension absolu pour récupérer les contraintes // dans la base absolue TenseurHH& sigHH = *(NevezTenseurHH(dimAbsolue)) ; int ni; // compteur globale de point d'integration bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde // donc il faut calculer tous les éléments de la métrique for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++) { // calcul des éléments de la métrique, on utilise le même calcul // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul); for (int i=1;i<= nbcomposante;i++) // récupération des coordonnées du tenseur en absolu *(entreePrinc->entree) >> sigHH.Coor(ordre(i,1),ordre(i,2)); // passage dans le repère local du tenseur contrainte final (*(tabSigHH(ni))) = sigHH.Baselocale((*(tabSigHH(ni))),*ex.giH_t); }; }; // retour des contraintes en absolu retour true si elle existe sinon false void ElemMeca::ContraintesEnAbsolues (bool cas,Tableau & tabSigHH,Tableau & tabSigAbs) { // dimensionnement de la metrique identique au cas d'un calcul explicite if( cas) { Tableau tab(7); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t; tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ; met->PlusInitVariables(tab) ; }; int dimAbsolue = ParaGlob::Dimension(); // dimension absolue // dimension des tenseurs //int dim = (*(tabSigHH(1))).Dimension(); // récup du nombre de composantes du tenseur int nbcomposante= ParaGlob::NbCompTens(); // redimensionnement éventuel du tableau de sortie int nbi = tabSigHH.Taille(); if (tabSigAbs.Taille() != nbi) tabSigAbs.Change_taille(nbi); for (int ni=1;ni<= nbi;ni++) if ( tabSigAbs(ni).Taille() != nbcomposante) tabSigAbs(ni).Change_taille(nbcomposante); // Il faut déterminer l'ordre dans lequel on parcours les contraintes qui doit // être compatible avec l'ordre des ddl Tableau2 ordre = OrdreContrainte(nbcomposante); // création d'un tenseur au dimension absolu pour récupérer les contraintes // dans la base absolue TenseurHH& sigHH = *(NevezTenseurHH(dimAbsolue)) ; int ni; // compteur globale de point d'integration bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde // donc il faut calculer tous les éléments de la métrique for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++) { // calcul des éléments de la métrique, on utilise le même calcul // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul); // passage dans le repère global du tenseur contrainte local sigHH = (*(tabSigHH(ni))).BaseAbsolue(sigHH,*ex.giB_t); // remplissage du vecteur pour le retour d'info for (int i=1;i<= nbcomposante;i++) tabSigAbs(ni)(i) = sigHH(ordre(i,1),ordre(i,2)); }; }; // ---------------- lecture écriture dans base info ---------------- // programmes utilisés par les classes dérivées // cas donne le niveau de la récupération // = 1 : on récupère tout // = 2 : on récupère uniquement les données variables (supposées comme telles) void ElemMeca::Lecture_bas_inf (ifstream& ent,const Tableau * tabMaillageNoeud,const int cas) { // appel de la routine d'élément Element::Lect_bas_inf_element(ent,tabMaillageNoeud,cas); switch (cas) { case 1 : // ------- on récupère tout ------------------------- { string toto,nom; // récup de la masse volumique ent >> toto >> masse_volumique ; // données thermique ent >> toto >> dilatation; // blocage éventuelle d'hourglass ent >> toto >> nom; type_stabHourglass=Id_Nom_StabHourglass(nom.c_str()); // blocage éventuel transversal ent >> nom; if (nom == "stabilisation_transversale:") { pt_StabMembBiel = new StabMembBiel; pt_StabMembBiel->Lecture_bas_inf(ent,cas); // string nom1; string nom_fonction; // ent >> nom1; // double valeur = 0.; // if (nom1=="par_valeur") {ent >> valeur;} // else {ent >> nom_fonction;}; // // maintenant il faut mettre à jour le conteneur // if ((pt_StabMembBiel == NULL)&&(nom1=="par_valeur")) // {pt_StabMembBiel=new StabMembBiel(valeur,NULL);} // else // {pt_StabMembBiel=new StabMembBiel(valeur,&nom_fonction);}; // // ensuite il faudra définir la fonction de stabilisation // // au moment de compléter l'élément } else // cas sans stabilisation {if (pt_StabMembBiel != NULL) {delete pt_StabMembBiel; pt_StabMembBiel=NULL;} }; break; } case 2 : // ----------- lecture uniquement de se qui varie -------------------- {break; } default : { cout << "\nErreur : valeur incorrecte du type de lecture !\n"; cout << "ElemMeca::Lecture_bas_inf (ifstream& ent,const int cas)" << " cas= " << cas << endl; Sortie(1); }; }; // ------ lecture dans tous les cas ------- // résultat d'erreur string toto; ent >> toto; if (toto == "erreur_de_contrainte") { if (sigErreur == NULL) sigErreur = new double; ent >> (*sigErreur) ; }; // données particulière pour les lois de comportement mécanique int tabSaveDonTaille = tabSaveDon.Taille(); if ((tabSaveDonTaille != 0) && (tabSaveDon(1) != NULL)) ent >> toto; int num ; // numéro du pt d'integ, n'est pas vraiment utilisé mais c'est mieux que de lire un string for (int i=1; i<= tabSaveDonTaille; i++) if (tabSaveDon(i) != NULL) { ent >> toto >> num; tabSaveDon(i)->Lecture_base_info(ent,cas); }; // données particulière pour les lois de comportement thermo physique int tabSaveTPTaille = tabSaveTP.Taille(); if ((tabSaveTPTaille != 0) && (tabSaveTP(1) != NULL)) ent >> toto; for (int i=1; i<= tabSaveTPTaille; i++) if (tabSaveTP(i) != NULL) { ent >> toto >> num; tabSaveTP(i)->Lecture_base_info(ent,cas); }; // données particulière pour la déformation mécanique int tabSaveDefDonTaille = tabSaveDefDon.Taille(); if ((tabSaveDefDonTaille != 0) && (tabSaveDefDon(1) != NULL)) ent >> toto; for (int i=1; i<= tabSaveDefDonTaille; i++) if (tabSaveDefDon(i) != NULL) { ent >> toto >> num; tabSaveDefDon(i)->Lecture_base_info(ent,cas); }; // ---- lecture des énergies // énergie totale ent >> toto >> energie_totale_t ; // par point d'intégration int energ_taille = tab_energ_t.Taille(); ent >> toto; for (int i=1;i<= energ_taille;i++) ent >> tab_energ_t(i) ; // énergie et puissance éventuelle de la partie bulk viscosity ent >> toto >> E_elem_bulk_t >> toto >> P_elem_bulk; E_elem_bulk_tdt = E_elem_bulk_t; // blocage éventuel transversal {int taille = 0; ent >> toto >> taille; if (taille != 0) {// cas avec stabilisation enregistrée // on lit systématiquement le cas 2 pt_StabMembBiel->Lecture_bas_inf(ent,2) ; // // on considère qu'il faut la récupérer // pt_StabMembBiel->Change_nb_pti(taille); // for (int i=1;i<= taille;i++) // {ent >> pt_StabMembBiel->FF_StabMembBiel(i); // pt_StabMembBiel->FF_StabMembBiel_t(i)=pt_StabMembBiel->FF_StabMembBiel(i); // }; }; // dans le cas où la taille est nulle, cela signifie qu'il n'y avait pas de stabilisation // auparavant, on laisse en l'état, car l'utilisateur a peut-être déjà définit une stabilisation }; // if (pt_StabMembBiel == NULL) // {// pas de stabilisation active, on passe la grandeur // ent >> toto >> toto; // } // else // {// on a une stabilisation enregistrée, on lit la stabilisation // int taille=0; // init // ent >> toto >> taille; // pt_StabMembBiel->Change_nb_pti(taille); // for (int i=1;i<= taille;i++) // {ent >> pt_StabMembBiel->FF_StabMembBiel(i); // pt_StabMembBiel->FF_StabMembBiel_t(i)=pt_StabMembBiel->FF_StabMembBiel(i); // }; // }; // les forces externes éventuelles int indic=0; ent >> toto >> indic; if (indic == 0) {if (lesChargeExtSurEle != NULL) delete lesChargeExtSurEle; } else {if (lesChargeExtSurEle == NULL) {lesChargeExtSurEle = new LesChargeExtSurElement;}; lesChargeExtSurEle->Lecture_base_info(ent,cas); }; }; // cas donne le niveau de sauvegarde // = 1 : on sauvegarde tout // = 2 : on sauvegarde uniquement les données variables (supposées comme telles) void ElemMeca::Ecriture_bas_inf(ofstream& sort,const int cas) { // appel de la routine d'élément Element::Ecri_bas_inf_element(sort,cas); // en fait ici on sauvegarde la même chose dans tous les cas, par contre la sortie // totale est documentée. switch (cas) { case 1 : // ------- on sauvegarde tout ------------------------- { // écriture de la masse volumique, sort << "masse_volumique " << masse_volumique <<" " << "\n"; // données thermique sort << "dilatation_thermique " << dilatation << " "; // blocage éventuel d'hourglass sort << "\n hourglass: " << Nom_StabHourglass(type_stabHourglass) << " "; // blocage éventuel transversal if (pt_StabMembBiel == NULL) {sort << "\n pas_de_stabilisation_transversale ";} else {sort << "\n stabilisation_transversale: "; pt_StabMembBiel->Ecriture_bas_inf(sort,cas); }; break; } case 2 : // ----------- sauvegarde uniquement de se qui varie -------------------- { break; } default : { cout << "\nErreur : valeur incorrecte du type d'écriture !\n"; cout << "ElemMeca::Ecriture_bas_inf(ofstream& sort,const int cas)" << " cas= " << cas << endl; Sortie(1); }; }; // --- informations à sortir dans tous les cas ---- // résultat d'erreur if (sigErreur != NULL) sort << "erreur_de_contrainte " << (*sigErreur) <<" " << "\n"; else sort << "pas_d'erreur_contrainte \n"; // données particulière pour les lois de comportement mécaniques int tabSaveDonTaille = tabSaveDon.Taille(); if ((tabSaveDonTaille != 0) && (tabSaveDon(1) != NULL)) sort << "\n\n +-+-+-data_spec_loi_comp_meca-+-+-+ "; for (int i=1; i<= tabSaveDonTaille; i++) if (tabSaveDon(i) != NULL) { if (i==1) {sort << "\n...pt_int= " <Ecriture_base_info(sort,cas);} // données particulière pour les lois de comportement thermo physiques int tabSaveTPTaille = tabSaveTP.Taille(); if ((tabSaveTPTaille != 0) && (tabSaveTP(1) != NULL)) sort << "\n\n +-+-+-data_spec_loi_comp_TP-+-+-+ "; for (int i=1; i<= tabSaveTPTaille; i++) if (tabSaveTP(i) != NULL) { if (i==1) {sort << "\n...pt_int= " <Ecriture_base_info(sort,cas);} // données particulière pour la déformation mécanique int tabSaveDefDonTaille = tabSaveDefDon.Taille(); if ((tabSaveDefDonTaille != 0) && (tabSaveDefDon(1) != NULL)) sort << "\n\n +-+-+-data_spec_def-+-+-+ "; for (int i=1; i<= tabSaveDefDonTaille; i++) if (tabSaveDefDon(i) != NULL) { if (i==1) {sort << "\n...pt_int= " <Ecriture_base_info(sort,cas); }; // --- sortie des énergies // énergie totale sort << "\n\n E_elem_totale: " << energie_totale << " "; // par point d'intégration int energ_taille = tab_energ_t.Taille(); sort << "\n par_pti: "; for (int i=1;i<= energ_taille;i++) sort << tab_energ_t(i) << " " ; // énergie et puissance éventuelle de la partie bulk viscosity sort << "\n E_el_bulk= " << E_elem_bulk_t << " P_el_bulk= " << P_elem_bulk; // blocage éventuel transversal if (pt_StabMembBiel == NULL) {sort << "\n StabTrans 0 ";} else {sort << "\n StabTrans 1 "; // on sort systématiquement le cas 2 pt_StabMembBiel->Ecriture_bas_inf(sort,2) ; }; // les forces externes éventuelles if (lesChargeExtSurEle != NULL) {sort << "\n lesChargesExt 1 "; lesChargeExtSurEle->Ecriture_base_info(sort,cas);} else {sort << "\n lesChargesExt 0 ";}; }; // calcul de la longueur d'arrête de l'élément minimal // divisé par la célérité la plus rapide dans le matériau // appelé par les classes dérivées // nb_noeud : =0 indique que l'on utilise tous les noeuds du tableau de noeuds // = un nombre > 0, indique le nombre de noeuds à utiliser au début du tableau double ElemMeca::Interne_Long_arrete_mini_sur_c(Enum_dure temps,int nb_noeud) { int nbne = tab_noeud.Taille(); // récup du nombre de noeud if (nb_noeud != 0) nbne = nb_noeud; // tout d'abord l'objectif est de déterminer la distance minimum // entre les différents noeuds // initialisation de la distance entre les deux noeuds double dist = ConstMath::tresgrand; for (int i=1;i<= nbne;i++) // on itère sur les noeuds restants switch (temps) { case TEMPS_0: for (int j=i+1;j<= nbne;j++) { double dist_new = (tab_noeud(i)->Coord0() - tab_noeud(j)->Coord0()).Norme(); if (dist_new < dist) dist = dist_new; } break; case TEMPS_t: for (int j=i+1;j<= nbne;j++) { double dist_new = (tab_noeud(i)->Coord1() - tab_noeud(j)->Coord1()).Norme(); if (dist_new < dist) dist = dist_new; } break; case TEMPS_tdt: for (int j=i+1;j<= nbne;j++) { double dist_new = (tab_noeud(i)->Coord2() - tab_noeud(j)->Coord2()).Norme(); if (dist_new < dist) dist = dist_new; } break; default : cout << "\n cas du temps non implante temps= " << Nom_dure(temps) << "\n ElemMeca::Interne_Long_arrete_mini_sur_c(Enum_dure temps)"; Sortie(1); } // traitement d'une erreur éventuelle if (dist <= ConstMath::petit) { cout << "\n **** ERREUR une longueur d'arrete de l'element est nulle" << "\n ElemMeca::Interne_Long_arrete_mini_sur_c(..." << "\n element: "; this->Affiche(1); #ifdef MISE_AU_POINT cout << "\n *** version mise au point: on continue neanmoins avec une longueur " << " arbitrairement tres petite (" <NbPti(); for (int i= 1; i<= taille; i++) (*lesPtIntegMecaInterne)(i).Change_statut_Invariants_contrainte(true); }; // idem pour la déformation void ElemMeca::ActivCalculInvariantsDeformation() { int taille = lesPtIntegMecaInterne->NbPti(); for (int i= 1; i<= taille; i++) (*lesPtIntegMecaInterne)(i).Change_statut_Invariants_deformation(true); }; // idem pour la vitesse de déformation void ElemMeca::ActivCalculInvariantsVitesseDeformation() { int taille = lesPtIntegMecaInterne->NbPti(); for (int i= 1; i<= taille; i++) (*lesPtIntegMecaInterne)(i).Change_statut_Invariants_vitesseDeformation(true); }; // initialisation pour le calcul de la matrice masse dans le cas de l'algorithme // de relaxation dynamique avec optimisation en continu de la matrice masse // casMass_relax: permet de choisir entre différentes méthodes de calcul de la masse void ElemMeca::InitCalculMatriceMassePourRelaxationDynamique(int ) { // on active le calcul des invariants de contraintes int taille = lesPtIntegMecaInterne->NbPti(); for (int i= 1; i<= taille; i++) (*lesPtIntegMecaInterne)(i).Change_statut_Invariants_contrainte(true); // on définit le ddl étendu correspondant, s'il n'existe pas // on définit le Ddl_enum_etendu correspondant à la masse au noeud pour la méthode if (masse_relax_dyn.Nom_vide() ) masse_relax_dyn = (Ddl_enum_etendu("masse_relax_dyn")); }; // phase de calcul de la matrice masse dans le cas de l'algo de relaxation dynamique // mi= lambda * delta_t**2 / 2 * (alpha*K+beta*mu+gamma*Isig/3+theta/2*Sig_mises) // avec delta_t=1 par défaut a priori, donc n'intervient pas ici // ep: epaisseur, K module de compressibilite, mu: module de cisaillement, Isig trace de sigma, // Sig_mises la contrainte de mises // casMass_relax: permet de choisir entre différentes méthodes de calcul de la masse void ElemMeca::CalculMatriceMassePourRelaxationDynamique (const double& alph, const double& beta, const double & lambda ,const double & gamma,const double & theta, int casMass_relax) { // 1) // on calcule la moyenne des grandeurs qui servent pour le calcul de la pseudo-masse que l'on distribue // également sur les noeuds double trace_moyenne= 0.; double mises_moyen= 0.; double compress_moy=0.; double cisaille_moy = 0.; int taille = lesPtIntegMecaInterne->NbPti(); // on boucle sur les points d'intégrations for (int i= 1; i<= taille; i++) {PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(i); const Vecteur & invariant = ptIntegMeca.SigInvar(); trace_moyenne += invariant(1) ; // le premier invariant c'est la trace if (invariant.Taille() > 1) mises_moyen += sqrt((3.*invariant(2)-invariant(1)*invariant(1))/2.); compress_moy += ptIntegMeca.ModuleCompressibilite(); cisaille_moy += ptIntegMeca.ModuleCisaillement(); }; // --- tout d'abord la partie contrainte // double ensemble = (val_propre_sup + val_propre_inf + val_cisaillement) / taille; // on divise par taille, ce qui conduit à la moyenne sur les points d'intégration double ensemble = (gamma*trace_moyenne /3.+ 0.5*theta*mises_moyen ) / taille; // --- maintenant la partie raideur ensemble += (alph * compress_moy + beta * cisaille_moy) / taille; // --- pour avoir l'intégration totale on utilise le volume // ensemble *= volume; // --- et maintenant on distribue sur tous les noeuds de manière identique // ne fonctionne que si tous les noeuds sont actifs !!! (ce qui est un cas particulier courant) // la masse élémentaire pour chaque noeud // double masse_elementaire = ensemble * lambda * 0.5 / tab_noeud.Taille(); // pour des membrannes on a: // mi= lambda * delta_t**2 / 2 * epaisseur/4 * (alpha*K+beta*mu+gamma*Isig/3+theta/2*Sig_mises) // avec delta_t=1 par défaut a priori, donc n'intervient pas ici // d'où le facteur 0.125, que l'on garde également pour le volume double masse_elementaire = ensemble * lambda * 0.125 ; // prise en compte du type d'élément: double ParNoeud = 1.; // n'est utilisé que par la méthode de Barnes classique int nbn=tab_noeud.Taille(); switch (Type_geom_generique(ElementGeometrique_const().TypeGeometrie())) { case VOLUME : // par rapport à la surface, on utilise la racine cubique du volume { double long_caracteristique = pow(Volume(),1./3.); masse_elementaire *= long_caracteristique; // ParNoeud = volume / nombre de noeuds, donc c'est la partie du volume // attribuée à chaque noeud ParNoeud = Volume() / nbn; break; } case SURFACE : // l'algorithme historique { double epaisseur_moyenne = EpaisseurMoyenne(TEMPS_tdt); masse_elementaire *= epaisseur_moyenne; // ParNoeud = section / nombre de noeuds, donc c'est la partie de la section // attribuée à chaque noeud ParNoeud = Volume() / epaisseur_moyenne/ nbn; break; } case LIGNE : { double section_moyenne = SectionMoyenne(TEMPS_tdt); double long_caracteristique = Volume()/section_moyenne; masse_elementaire *= long_caracteristique; // ParNoeud = longueur / nombre de noeuds, donc c'est la partie de la longueur // attribuée à chaque noeud ParNoeud = Volume() / section_moyenne/ nbn; break; } default : cout << "\nErreur : pour l'instant les types autres que volume, surface, ligne ne sont pas pris en compte dans la relaxation " << " dynamique selon barnes !\n"; cout << "\n ElemMeca::CalculMatriceMassePourRelaxationDynamique(.. \n"; Sortie(1); }; // on parcours les noeuds de l'élément et on remplace le cas échéent, la valeur de la masse au noeud. switch (casMass_relax) { case 1: case 3: // cas 1 : on cumule aux noeuds, cas 3 on fait la moyenne (pas effectué ici, on ne fait // que préparer le travail pour le cas 3) {for (int inoe=1;inoe<=nbn;inoe++) {Noeud& noe = *tab_noeud(inoe); // on cumule noe.ModifDdl_etendu(masse_relax_dyn).Valeur() += masse_elementaire; }; break; } case 2: // cas 2 : on prend la masse maxi {for (int inoe=1;inoe<=nbn;inoe++) {Noeud& noe = *tab_noeud(inoe); const Ddl_etendu& masse_actuel =noe.DdlEtendue(masse_relax_dyn); // dans le cas où la masse actuelle est plus petite on la remplace if (masse_actuel.ConstValeur() <= masse_elementaire) noe.ModifDdl_etendu(masse_relax_dyn).Valeur() = masse_elementaire; }; break; } case 4: case 5: // on cumule (puis on moyennera pour le cas 4, autre part) et on divise par la section comme dans la formule de barnes {for (int inoe=1;inoe<=nbn;inoe++) {Noeud& noe = *tab_noeud(inoe); const Ddl_etendu& masse_actuel =noe.DdlEtendue(masse_relax_dyn); // on cumule if (masse_actuel.ConstValeur() <= masse_elementaire) // !!!! je pense que ce cas est très bizarre il faudra revoir ce que cela signifie.... !! noe.ModifDdl_etendu(masse_relax_dyn).Valeur() += (masse_elementaire/ParNoeud); }; break; } default : cout << "\nErreur : pour l'instant le cas casMass_relax= " << casMass_relax << " n'est pas pris en compte !\n"; cout << "\n ElemMeca::CalculMatriceMassePourRelaxationDynamique(.. \n"; Sortie(1); }; //-- fin du switch }; // METHODES VIRTUELLES: // recuperation des coordonnées du point de numéro d'ordre = iteg pour // la grandeur enu // temps: dit si c'est à 0 ou t ou tdt // si erreur retourne erreur à true // utilisé par les classes dérivées Coordonnee ElemMeca::CoordPtInt(Enum_dure temps,Enum_ddl enu,int iteg,bool& err) { Coordonnee ptret;err = false; // récupération de l'élément géométrique correspondant à Enu ElemGeomC0& ele = this->ElementGeometrie(enu); // récupération du tableau des fonctions d'interpolations const Tableau & tabphi = ele.TaPhi(); if ((iteg < 0) || (iteg > tabphi.Taille())) { err = true;} else { switch (temps) { case TEMPS_0 : ptret = met->PointM_0(tab_noeud,tabphi(iteg)); break; case TEMPS_t : ptret = met->PointM_t(tab_noeud,tabphi(iteg)); break; case TEMPS_tdt : ptret = met->PointM_tdt(tab_noeud,tabphi(iteg)); break; }; }; return ptret; // retour }; // recuperation des coordonnées du point d'intégration numéro = iteg pour // la face : face // temps: dit si c'est à 0 ou t ou tdt // si erreur retourne erreur à true Coordonnee ElemMeca::CoordPtIntFace(int face, Enum_dure temps,int iteg,bool& err) { Coordonnee ptret(ParaGlob::Dimension());err = false; Deformation* definter=NULL; // variable de travail transitoire if (!(ParaGlob::AxiSymetrie())) {if (defSurf.Taille()) {definter = defSurf(face);} // le cas normal est non axisymétrique else {err = true;}; } else // en axisymétrie, c'est une def d'arête { if (defArete.Taille()) {definter = defArete(face);} else {err = true;}; }; if (definter == NULL) err = true; if (!err) { Deformation & defS = *definter; // pour simplifier l'écriture defS.ChangeNumInteg(iteg); switch (temps) { case TEMPS_0 : ptret = defS.Position_0(); break; case TEMPS_t : ptret = defS.Position_t(); break; case TEMPS_tdt : ptret = defS.Position_tdt(); break; }; }; return ptret; // retour }; // recuperation des coordonnées du point d'intégration numéro = iteg pour // la face : face // temps: dit si c'est à 0 ou t ou tdt // si erreur retourne erreur à true Coordonnee ElemMeca::CoordPtIntArete(int arete, Enum_dure temps,int iteg,bool& err) { Coordonnee ptret(ParaGlob::Dimension());err = false; Deformation* definter=NULL; // variable de travail transitoire if (defArete.Taille()) {definter = defArete(arete);} else {err = true;}; if (definter == NULL) err = true; if (!err) { Deformation & defA = *definter; // pour simplifier l'écriture defA.ChangeNumInteg(iteg); switch (temps) { case TEMPS_0 : ptret = defA.Position_0(); break; case TEMPS_t : ptret = defA.Position_t(); break; case TEMPS_tdt : ptret = defA.Position_tdt(); break; }; }; return ptret; // retour }; // retourne le numero du pt d'ing le plus près ou est exprimé la grandeur enum // temps: dit si c'est à 0 ou t ou tdt // utilisé par les classes dérivées int ElemMeca::PtLePlusPres(Enum_dure temps,Enum_ddl enu, const Coordonnee& M) { int iret; // récupération de l'élément géométrique correspondant à Enu ElemGeomC0& ele = this->ElementGeometrie(enu); // récupération du tableau des fonctions d'interpolations const Tableau & tabphi = ele.TaPhi(); // on boucle sur les pt d'integ pour avoir le point le plus près int tabphitaille = tabphi.Taille(); Coordonnee P; iret=1; double dist= ConstMath::tresgrand; for (int ipt = 1;ipt<=tabphitaille;ipt++) { switch (temps) { case TEMPS_0 : P = met->PointM_0(tab_noeud,tabphi(ipt)); break; case TEMPS_t : P = met->PointM_t(tab_noeud,tabphi(ipt)); break; case TEMPS_tdt : P = met->PointM_tdt(tab_noeud,tabphi(ipt)); break; }; double di=(M-P).Norme(); if (di < dist) { dist = di; iret = ipt;}; }; return iret; // retour }; // ==== >>>> methodes virtuelles redéfini éventuellement dans les classes dérivées ============ // ramene l'element geometrique correspondant au ddl passé en paramètre ElemGeomC0& ElemMeca::ElementGeometrie(Enum_ddl ddl) const { Enum_ddl enu = PremierDdlFamille(ddl); switch (enu) { case X1 : return ElementGeometrique(); break; case SIG11 : return ElementGeometrique(); break; case EPS11 : return ElementGeometrique(); break; case DEPS11 : return ElementGeometrique(); break; case ERREUR : return ElementGeometrique(); break; case TEMP : return ElementGeometrique(); break; case V1 : return ElementGeometrique(); break; case GAMMA1 : return ElementGeometrique(); break; case DELTA_TEMP : return ElementGeometrique(); break; case R_TEMP : return ElementGeometrique(); break; case R_X1 : return ElementGeometrique(); break; case R_V1 : return ElementGeometrique(); break; case R_GAMMA1 : return ElementGeometrique(); break; default : {cout << "\n cas non prevu ou non encore implante: ddl= " << Nom_ddl(ddl) << "\n ElemMeca::ElementGeometrie(Enum_ddl ddl) " ; Sortie(1);} } return ElementGeometrique(); // pour taire le compilo }; // ramène le nombre de grandeurs génératrices pour un pt d'integ, correspondant à un type enuméré // peut-être surchargé pour des éléments particuliers int ElemMeca::NbGrandeurGene(Enum_ddl ddl) const { Enum_ddl enu = PremierDdlFamille(ddl); int nbGG = 0.; // par défaut switch (enu) { case X1 : case SIG11 : case EPS11 : case DEPS11 : { // cas d'un calcul de mécanique classique switch (Type_geom_generique(id_geom)) { case LIGNE : case POINT_G : nbGG = 1; break; // SIG11 case SURFACE : nbGG = 3; break; // SIG11, SIG22, SIG12 case VOLUME : nbGG = 6; break; // les 6 composantes default : cout << "\nErreur : cas non traite, id_geom= :" << Nom_type_geom(Type_geom_generique(id_geom)) << "ElemMeca::NbGrandeurGene(.. \n"; Sortie(1); }; break; } default : {cout << "\n cas non prevu ou non encore implante: ddl= " << Nom_ddl(ddl) << "\n ElemMeca::NbGrandeurGene(Enum_ddl ddl) " ; Sortie(1);} }; return nbGG; }; // modification de l'orientation de l'élément en fonction de cas_orientation // =0: inversion simple (sans condition) de l'orientation // si cas_orientation est diff de 0: on calcul le jacobien aux différents points d'intégration // 1. si tous les jacobiens sont négatifs on change d'orientation // 2. si tous les jacobiens sont positifs on ne fait rien // 3. si certains jacobiens sont positifs et d'autres négatifs message // d'erreur et on ne fait rien // ramène true: s'il y a eu changement effectif, sinon false bool ElemMeca::Modif_orient_elem(int cas_orientation) { // retour: bool retour=false; // par défaut pas de changement if (cas_orientation == 0) // cas où on inverse l'orientation sans condition particulière { // on change l'orientation de l'élément retour = true; int nbnoe = tab_noeud.Taille(); Tableau tab_inter(tab_noeud); // on crée un tableau intermédiaire // on récupère la numérotation locale inverse const Tableau tabi = ElementGeometrique().InvConnec(); // on met à jour le tableau actuel for (int n=1;n<=nbnoe;n++) tab_noeud(n)=tab_inter(tabi(n)); } else { // si cas_orientation est diff de 0: on calcul le jacobien aux différents points d'intégration // 1. si tous les jacobiens sont négatifs on change d'orientation // 2. si tous les jacobiens sont positifs on ne fait rien // 3. si certains jacobiens sont positifs et d'autres négatifs message // d'erreur et on ne fait rien int cas=1; // a priori tout est ok // boucle sur les pt d'integ for (def->PremierPtInteg();def->DernierPtInteg();def->NevezPtInteg()) { double jacobien_0 = def->JacobienInitial(); if (jacobien_0 < 0) {if (cas ==1) // on a trouvé un jacobien négatif {cas =0; } }// si c'était positif --> négatif else // cas positif {if (cas == 0) // on a déjà changé { cas = 2; break;} // on sort de la boucle }; }; // gestion de pb if (cas == 2) { cout << "\n **** Attention **** element nb= "<< this->Num_elt() << " peut-etre trop distordu ?" << " pt d'integ meca positif et negatif !! "; if (ParaGlob::NiveauImpression() > 7) {// on va essayer de sortir plus d'info int pti=1;cout << "\n les jacobiens initiaux \n"; // on sort les valeurs des jacobiens for (def->PremierPtInteg();def->DernierPtInteg();def->NevezPtInteg(),pti++) { cout << " pti= "<Coord0(); }; }; } else if (cas == 0) { switch (cas_orientation) { case 1: // on change l'orientation de l'élément {retour = true; int nbnoe = tab_noeud.Taille(); Tableau tab_inter(tab_noeud); // on crée un tableau intermédiaire // on récupère la numérotation locale inverse const Tableau tabi = ElementGeometrique().InvConnec(); for (int n=1;n<=nbnoe;n++) tab_noeud(n)=tab_inter(tabi(n)); break; } case -1: // on sort une information à l'écran { cout << "\n element nb= "<< this->Num_elt() << "jacobien negatif " ; break; } default: cout << "\n erreur le cas : " << cas_orientation << " n'est pas actuellement pris en compte" << "\n ElemMeca::Modif_orient_elem(..."; Sortie(1); }; }; }; // retour return retour; }; // récupération de valeurs interpolées pour les grandeur enu ou directement calculées // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière // une seule des 3 métriques doit-être renseigné, les autres doivent être un pointeur nul Tableau ElemMeca::Valeur_multi_interpoler_ou_calculer (bool absolue, Enum_dure temps,const List_io& enu ,Deformation & defor ,const Met_abstraite::Impli* ex_impli ,const Met_abstraite::Expli_t_tdt* ex_expli_tdt ,const Met_abstraite::Expli* ex_expli ) { // 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) // 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; // 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 Coordonnee* Mtdt = NULL; // coordonnées finales éventuelles du point d'intégration considéré Coordonnee* Mt=NULL; // coordonnées à t Coordonnee* M0 = NULL; // coordonnées initiales éventuelles du point d'intégration considéré Coordonnee* N_surf = NULL; // coordonnée d'un vecteur normal actuel 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 Tableau tab_ret (enu.size()); // éléments de métrique et matrices de passage TenseurHH* gijHH;TenseurBB* gijBB;BaseB* giB; BaseH* giH_0;BaseH* giH; BaseB* giB_0;BaseB* giB_t; if (ex_impli != NULL) {gijHH = ex_impli->gijHH_tdt;gijBB = ex_impli->gijBB_tdt; giB = ex_impli->giB_tdt; giH_0 = ex_impli->giH_0;giH = ex_impli->giH_tdt; giB_0 = ex_impli->giB_0;giB_t = ex_impli->giB_t; } else if (ex_expli_tdt != NULL) {gijHH = ex_expli_tdt->gijHH_tdt;gijBB = ex_expli_tdt->gijBB_tdt; giB = ex_expli_tdt->giB_tdt; giH_0 = ex_expli_tdt->giH_0;giH = ex_expli_tdt->giH_tdt; giB_0 = ex_expli_tdt->giB_0; giB_t = ex_expli_tdt->giB_t; } else if (ex_expli != NULL) {gijHH = ex_expli->gijHH_t;gijBB = ex_expli->gijBB_t; giB = giB_t = ex_expli->giB_t; giH_0 = ex_expli->giH_0;giH = ex_expli->giH_t; giB_0 = ex_expli->giB_0; } else { cout << "\n *** cas non prevu : aucune metrique transmise " << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(..." << endl; Sortie(1); }; // 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_coordonnees = false; bool besoin_deplacements = false; bool besoin_coordonnees_t = false;bool besoin_coordonnees_t0 = false; for (ie=enu.begin(); ie!=iefin;ie++) { int posi = (*ie).Position()-NbEnum_ddl(); switch (posi) {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; }; }; // définition des tenseurs si nécessaire // ----- maintenant on calcule les grandeurs nécessaires ----- if (besoin_coordonnees) {Mtdt = new Coordonnee(ParaGlob::Dimension()); *Mtdt = defor.Position_tdt(); }; if (besoin_coordonnees_t ) {*Mt = defor.Position_tdt(); }; if (besoin_deplacements || besoin_coordonnees_t0) {if (M0 == NULL) M0 = new Coordonnee(ParaGlob::Dimension()); (*M0) = defor.Position_0(); }; if (Vitesse != NULL) {Vitesse = new Coordonnee(ParaGlob::Dimension()); (*Vitesse) = defor.VitesseM_tdt(); }; // 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'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) = Util::ProdVec_coorBN( (*giB)(1), (*giB)(2)); N_surf->Normer(); // que l'on norme }; }; // 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 }; // 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 }; //----- 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++) { 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)) ) { // on recherche en générale une interpolation en fonction des noeuds: il faut donc que la grandeur soit // présente aux noeuds !! // 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())) { tab_ret(it)= defor.DonneeInterpoleeScalaire((*ie).Enum(),temps); } else if ((Meme_famille((*ie).Enum(),EPS11)) && ((*ie).Nom_vide())) { tab_ret(it)= defor.DonneeInterpoleeScalaire((*ie).Enum(),temps); } else if ((Meme_famille((*ie).Enum(),DEPS11)) && ((*ie).Nom_vide())) { tab_ret(it)= defor.DonneeInterpoleeScalaire((*ie).Enum(),temps); } 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 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;} 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_interpoler_ou_calculer(..."; tab_ret(it) = 0.; }; } // fin cas **** 2 >>>>> } // " " " } // -- fin 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_interpoler_ou_calculer(..."; // 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_interpoler_ou_calculer(..."; }; };// -- fin de la boucle sur la liste de Ddl_enum_etendu // delet e des tenseurs if (Mtdt != NULL) delete Mtdt; // coordonnée du point à tdt 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 à t à la surface if (N_surf_t0 != NULL) delete N_surf_t0; // vecteur normal à la surface if (Vitesse != NULL) delete Vitesse; // vitesse // liberation des tenseurs intermediaires LibereTenseur(); return tab_ret; }; // récupération de valeurs interpolées pour les grandeur enu // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière // une seule des 3 métriques doit-être renseigné, les autres doivent être un pointeur nul void ElemMeca::Valeurs_Tensorielles_interpoler_ou_calculer (bool absolue, Enum_dure temps,List_io& enu ,Deformation & defor ,const Met_abstraite::Impli* ex_impli ,const Met_abstraite::Expli_t_tdt* ex_expli_tdt ,const Met_abstraite::Expli* ex_expli ) { // ----- 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; // éléments de métrique et matrices de passage TenseurHH* gijHH;TenseurBB* gijBB;BaseB* giB; BaseH* giH_0;BaseH* giH; BaseB* giB_0;BaseB* giB_t; if (ex_impli != NULL) {gijHH = ex_impli->gijHH_tdt;gijBB = ex_impli->gijBB_tdt; giB = ex_impli->giB_tdt; giH_0 = ex_impli->giH_0;giH = ex_impli->giH_tdt; giB_t = ex_impli->giB_t; giB_0 = ex_impli->giB_0; } else if (ex_expli_tdt != NULL) {gijHH = ex_expli_tdt->gijHH_tdt;gijBB = ex_expli_tdt->gijBB_tdt; giB = ex_expli_tdt->giB_tdt; giH_0 = ex_expli_tdt->giH_0;giH = ex_expli_tdt->giH_tdt; giB_t = ex_expli_tdt->giB_t; giB_0 = ex_expli_tdt->giB_0; } else if (ex_expli != NULL) {gijHH = ex_expli->gijHH_t;gijBB = ex_expli->gijBB_t; giB = giB_t = ex_expli->giB_t; giH_0 = ex_expli->giH_0;giH = ex_expli->giH_t; giB_0 = ex_expli->giB_0; } else { cout << "\n *** cas non prevu : aucune metrique transmise " << "\n ElemMeca::Valeurs_Tensorielles_interpoler_ou_calculer(..." << endl; Sortie(1); }; // def de tenseurs pour la sortie // les tenseurs restants en locale 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 // pour les vecteurs propres // pour les bases Tableau base_ad_hoc , base_giH , base_giB; // --- 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 besoin_coordonnees = false; bool besoin_coordonnees_t = false;bool besoin_coordonnees_t0 = false; bool besoin_rep_local_ortho=false; bool besoin_rep_giH = false; bool besoin_rep_giB = 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 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 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 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;} 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_interpoler_ou_calculer(...."; }; } }; }; // récup des bases si besoin if (besoin_rep_local_ortho) {// on calcule éventuellement la matrice de passage Mat_pleine Aa0(dim,dim),Aafin(dim,dim); if ( ((!absolue)&&(dim_espace==3)&&(dim==2)) || ((!absolue)&&(dim_espace>1)&&(dim==1)) ) { def->BasePassage(absolue,*giB_0,*giB,*giH_0,*giH,Aa0,Aafin); // g^i = Aa^i_{.a} * Ip^a // et on a : Ip_a = beta_a^{.j} g_j = [A^j_{.a}]^T g_j }; 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)) = Aafin(1,1) * giB->Coordo(1) + Aafin(2,1) * giB->Coordo(2); *(base_ad_hoc(2)) = Aafin(1,2) * giB->Coordo(1) + Aafin(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)) = Aafin(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 if (besoin_coordonnees) (*Mtdt) = defor.Position_tdt(); if (besoin_coordonnees_t) (*Mt) = defor.Position_t(); if (besoin_coordonnees_t0) (*M0) = defor.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 ElemMeca::Valeurs_Tensorielles_interpoler_ou_calculer(... "; }; 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) - defor.Position_0(); if (Vitesse != NULL) (*Vitesse) = defor.VitesseM_tdt(); // --- cas des grandeurs de la décomposition polaire // delete des tenseurs // liberation des tenseurs intermediaires LibereTenseur(); }; // mise à jour éventuel de repère d'anisotropie void ElemMeca::Mise_a_jour_repere_anisotropie (BlocGen & bloc,LesFonctions_nD* lesFonctionsnD) { // cas de la définition d'un repère d'anisotropie // pour l'instant le repère d'anisotropie est éventuellement // utilisable par la loi de comportement // 1) tout d'abord on doit définir le repère au niveau de chaque point d'intégration // définition d'un repère d'anisotropie à un point d'intégration int nb_pti = lesPtIntegMecaInterne->NbPti(); for (int ni=1;ni<= nb_pti;ni++) {BaseH repH = DefRepereAnisotropie(ni,lesFonctionsnD,bloc); // on crée les grandeurs de passages pour renseigner les grandeurs particulières // de la loi de comportement int dim = ParaGlob::Dimension(); Tableau tab_coor(dim); for (int i=1;i<=dim;i++) tab_coor(i)=repH.Coordo(i); // passage des infos: on utilise la méthode Complete_SaveResul // qui au premier passage via ElemMeca::Complete_ElemMeca( // a créé les conteneurs et au passage courant met à jour seulement if (tabSaveDon(ni) != NULL) tabSaveDon(ni)->Complete_SaveResul(bloc,tab_coor,loiComp); }; };