// FICHIER : HexaMemb.cc // CLASSE : HexaMemb // 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 using namespace std; //introduces namespace std #include #include "Sortie.h" #include "FrontQuadLine.h" #include "FrontPointF.h" #include "HexaMemb.h" #include "TypeConsTens.h" // -------------- definition de la classe conteneur de donnees communes ------------ HexaMemb::DonnComHexa::DonnComHexa (ElemGeomC0* hexa1,DdlElement& tab,DdlElement& tabErr,DdlElement& tab_Err1Sig, Met_abstraite& met_gene,Tableau & resEr,Mat_pleine& raidEr, ElemGeomC0* hexaeEr,GeomQuadrangle& quadS,GeomSeg& seggS, Vecteur& residu_int,Mat_pleine& raideur_int, Tableau & residus_extN,Tableau & raideurs_extN, Tableau & residus_extA,Tableau & raideurs_extA, Tableau & residus_extS,Tableau & raideurs_extS, Mat_pleine& mat_masse,ElemGeomC0* hexaeMas,int nbi, ElemGeomC0* hexaeHourg ) : tab_ddl(tab),tab_ddlErr(tabErr),tab_Err1Sig11(tab_Err1Sig),metrique(met_gene) ,matGeom(tab.NbDdl(),tab.NbDdl()),matInit(tab.NbDdl(),tab.NbDdl()) ,d_epsBB(tab.NbDdl()),d_sigHH(tab.NbDdl()),d2_epsBB(nbi) ,resErr(resEr),raidErr(raidEr),quadraS(quadS),segS(seggS) ,residu_interne(residu_int),raideur_interne(raideur_int) ,residus_externeN(residus_extN),raideurs_externeN(raideurs_extN) ,residus_externeA(residus_extA),raideurs_externeA(raideurs_extA) ,residus_externeS(residus_extS),raideurs_externeS(raideurs_extS) ,matrice_masse(mat_masse) ,hexaed(hexa1),hexaedEr(hexaeEr),hexaedMas(hexaeMas),hexaedHourg(hexaeHourg) { int nbddl = tab.NbDdl(); for (int ni=1;ni<=nbi;ni++) {d2_epsBB(ni).Change_taille(nbddl); for (int i1=1; i1<= nbddl; i1++) for (int i2=1; i2<= nbddl; i2++) d2_epsBB(ni)(i1,i2) = NevezTenseurBB(3); }; int tailledeps = d_epsBB.Taille(); for (int i=1;i<= tailledeps; i++) d_epsBB(i) = NevezTenseurBB(3); int tailledsig = d_sigHH.Taille(); for (int j=1;j<= tailledsig; j++) d_sigHH(j) = NevezTenseurHH(3); }; // Il ne doit exister qu'un exemplaire de la classe, donc au niveau du constructeur de copie les tableaux // de pointeurs pointes sur les mêmes entités que la valeur passée en paramètre HexaMemb::DonnComHexa::DonnComHexa(DonnComHexa& a) : hexaed(a.hexaed), tab_ddl(a.tab_ddl),tab_ddlErr(a.tab_ddlErr) ,tab_Err1Sig11(a.tab_Err1Sig11), metrique(a.metrique) ,matGeom(a.matGeom),matInit(a.matInit),d2_epsBB(a.d2_epsBB) ,resErr(a.resErr),raidErr(a.raidErr),quadraS(a.quadraS),segS(a.segS) ,hexaedEr(a.hexaedEr) ,d_epsBB(a.d_epsBB),d_sigHH(a.d_sigHH) ,residu_interne(a.residu_interne),raideur_interne(a.raideur_interne) ,residus_externeN(a.residus_externeN),raideurs_externeN(a.raideurs_externeN) ,residus_externeA(a.residus_externeA),raideurs_externeA(a.raideurs_externeA) ,residus_externeS(a.residus_externeS),raideurs_externeS(a.raideurs_externeS) ,matrice_masse(a.matrice_masse),hexaedMas(a.hexaedMas) ,hexaedHourg(a.hexaedHourg) { int nbddl = d_sigHH.Taille(); int nbi=d2_epsBB.Taille(); for (int ni=1;ni<=nbi;ni++) for (int i1=1; i1<= nbddl; i1++) for (int i2=1; i2<= nbddl; i2++) d2_epsBB(ni)(i1,i2) = NevezTenseurBB (*(a.d2_epsBB(ni)(i1,i2))); int tailledeps = d_epsBB.Taille(); for (int i=1;i<= tailledeps; i++) d_epsBB(i) = NevezTenseurBB (3); int tailledsig = d_sigHH.Taille(); for (int j=1;j<= tailledsig; j++) d_sigHH(j) = NevezTenseurHH (3); }; HexaMemb::DonnComHexa::~DonnComHexa() { int nbddl = tab_ddl.NbDdl(); int nbi=d2_epsBB.Taille(); for (int ni=1;ni<=nbi;ni++) for (int i1=1; i1<= nbddl; i1++) for (int i2=1; i2<= nbddl; i2++) delete d2_epsBB(ni)(i1,i2); int tailledeps = d_epsBB.Taille(); for (int i=1;i<= tailledeps; i++) delete d_epsBB(i); int tailledsig = d_sigHH.Taille(); for (int j=1;j<= tailledsig; j++) delete d_sigHH(j); }; // ---------- fin definition de la classe conteneur de donnees communes ------------ // -+-+ definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+ HexaMemb::UneFois::UneFois () : // constructeur par défaut doCoMemb(NULL),CalResPrem_t(0),CalResPrem_tdt(0),CalimpPrem(0),dualSortHexa(0) ,CalSMlin_t(0),CalSMlin_tdt(0),CalSMRlin(0) ,CalSMsurf_t(0),CalSMsurf_tdt(0),CalSMRsurf(0) ,CalSMvol_t(0),CalSMvol_tdt(0),CalSMvol(0) ,CalDynamique(0),CalPt_0_t_tdt(0) ,nbelem_in_Prog(0) {}; HexaMemb::UneFois::~UneFois () { delete doCoMemb; }; // -+-+ fin definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+ // CONSTRUCTEURS : // Constructeur par defaut HexaMemb::HexaMemb () : ElemMeca(),lesPtMecaInt(),unefois(NULL),nombre(NULL) { lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca }; // Constructeur fonction d'un numero // d'identification , d'identificateur d'interpolation et de geometrie HexaMemb::HexaMemb (int num_mail,int num_id,Enum_interpol id_interp_elt,Enum_geom id_geom_elt,string info) : ElemMeca(num_mail,num_id,id_interp_elt,id_geom_elt,info) ,lesPtMecaInt(),unefois(NULL),nombre(NULL) { lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca }; // Constructeur fonction d'un numero de maillage et d'identification, // du tableau de connexite des noeuds, d'identificateur d'interpolation et de geometrie HexaMemb::HexaMemb (int num_mail,int num_id,Enum_interpol id_interp_elt,Enum_geom id_geom_elt, const Tableau& tab,string info) : ElemMeca(num_mail,num_id,tab,id_interp_elt,id_geom_elt,info) ,lesPtMecaInt(),unefois(NULL),nombre(NULL) { lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca }; // Constructeur de copie HexaMemb::HexaMemb (const HexaMemb& hexaMemb) : ElemMeca (hexaMemb) ,lesPtMecaInt(hexaMemb.lesPtMecaInt) ,unefois(hexaMemb.unefois),nombre(hexaMemb.nombre) { lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca // --- cas des différentes puissances et ... --- HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture residu = &(CoHexa->residu_interne); // residu local raideur = &(CoHexa->raideur_interne); // raideur locale // --- cas de la dynamique ----- mat_masse = &(CoHexa->matrice_masse); // --- cas des efforts externes concernant noeuds ------ res_extN = &(CoHexa->residus_externeN); // pour les résidus et second membres raid_extN= &(CoHexa->raideurs_externeN);// pour les raideurs // --- cas des efforts externes concernant les aretes ------ res_extA = &(CoHexa->residus_externeA); // pour les résidus et second membres raid_extA= &(CoHexa->raideurs_externeA);// pour les raideurs // --- cas des efforts externes concernant les faces ------ res_extS= &(CoHexa->residus_externeS); // pour les résidus et second membres raid_extS= &(CoHexa->raideurs_externeS); // pour les raideurs }; // DESTRUCTEUR : HexaMemb::~HexaMemb () { // on va libérer les déformations d'arrête et de surface qui sont un peu particulières /* if (defArete.Taille() != 0) if (defArete(1) != NULL) { delete defArete(1); // on remet à null les pointeurs pour un appel correcte du destructeur d'elemMeca for (int ia=1;ia<= 12; ia++) defArete(ia) = NULL; } if (defSurf.Taille() != 0) if (defSurf(1) != NULL) { delete defSurf(1); // on remet à null les pointeurs pour un appel correcte du destructeur d'elemMeca for (int ia=1;ia<= 6; ia++) defSurf(ia) = NULL; } */ LibereTenseur(); }; // Lecture des donnees de la classe sur fichier void HexaMemb::LectureDonneesParticulieres (UtilLecture * entreePrinc,Tableau * tabMaillageNoeud) { int nb;int nbne=nombre->nbne; tab_noeud.Change_taille(nbne); for (int i=1; i<= nbne; i++) { *(entreePrinc->entree) >> nb; if ((entreePrinc->entree)->rdstate() == 0) // pour mémoire ici on a /* enum io_state { badbit = 1<<0, // -> 1 dans rdstate() eofbit = 1<<1, // -> 2 failbit = 1<<2, // -> 4 goodbit = 0 // -> O };*/ tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale #ifdef ENLINUX else if ((entreePrinc->entree)->fail()) // on a atteind la fin de la ligne et on appelle un nouvel enregistrement { entreePrinc->NouvelleDonnee(); // lecture d'un nouvelle enregistrement *(entreePrinc->entree) >> nb; tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale } #else /* #ifdef SYSTEM_MAC_OS_X_unix else if ((entreePrinc->entree)->fail()) // on a atteind la fin de la ligne et on appelle un nouvel enregistrement { entreePrinc->NouvelleDonnee(); // lecture d'un nouvelle enregistrement *(entreePrinc->entree) >> nb; tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale } #else*/ else if ((entreePrinc->entree)->eof()) // la lecture est bonne mais on a atteind la fin de la ligne { tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture // si ce n'est pas la fin de la lecture on appelle un nouvel enregistrement if (i != nbne) entreePrinc->NouvelleDonnee(); // lecture d'un nouvelle enregistrement } // #endif #endif else // cas d'une erreur de lecture { cout << "\n erreur de lecture inconnue "; entreePrinc->MessageBuffer("** lecture des données particulières **"); cout << "HexaMemb::LecDonPart"; Affiche(); Sortie (1); } } // construction du tableau de ddl des noeuds de HexaMemb ConstTabDdl(); }; // calcul d'un point dans l'élément réel en fonction des coordonnées dans l'élément de référence associé // temps: indique si l'on veut les coordonnées à t = 0, ou t ou tdt // 1) cas où l'on utilise la place passée en argument Coordonnee & HexaMemb::Point_physique(const Coordonnee& c_int,Coordonnee & co,Enum_dure temps) { HexaMemb::DonnComHexa* doCoHexa = unefois->doCoMemb; // pour simplifier l'écriture // a) on commence par définir les bonnes grandeurs dans la métrique if( !(unefois->CalPt_0_t_tdt )) { unefois->CalPt_0_t_tdt += 1; Tableau tab(3); tab(1)=iM0;tab(2)=iMt;tab(3)=iMtdt; doCoHexa->metrique.PlusInitVariables(tab) ; }; // b) calcul de l'interpolation const Vecteur& phi = doCoHexa->hexaed->Phi_point(c_int); // c) calcul du point switch (temps) { case TEMPS_0 : co = met->PointM_0(tab_noeud,phi); break; case TEMPS_t : co = met->PointM_t(tab_noeud,phi); break; case TEMPS_tdt : co = met->PointM_tdt(tab_noeud,phi); break; } // d) retour return co; }; // 3) cas où l'on veut les coordonnées aux 1, 2 ou trois temps selon la taille du tableau t_co void HexaMemb::Point_physique(const Coordonnee& c_int,Tableau & t_co) { HexaMemb::DonnComHexa* doCoHexa = unefois->doCoMemb; // pour simplifier l'écriture // a) on commence par définir les bonnes grandeurs dans la métrique if( !(unefois->CalPt_0_t_tdt )) { unefois->CalPt_0_t_tdt += 1; Tableau tab(3); tab(1)=iM0;tab(2)=iMt;tab(3)=iMtdt; doCoHexa->metrique.PlusInitVariables(tab) ; }; // b) calcul de l'interpolation const Vecteur& phi = doCoHexa->hexaed->Phi_point(c_int); // c) calcul des point switch (t_co.Taille()) { case 3 : t_co(3) = met->PointM_tdt(tab_noeud,phi); case 2 : t_co(2) = met->PointM_t(tab_noeud,phi); case 1 : t_co(1) = met->PointM_0(tab_noeud,phi); } }; // Calcul du residu local à t ou tdt en fonction du booleen atdt Vecteur* HexaMemb::CalculResidu (bool atdt,const ParaAlgoControle & pa) { HexaMemb::DonnComHexa* doCoHexa = unefois->doCoMemb; // pour simplifier l'écriture Tableau & d_epsBB = unefois->doCoMemb->d_epsBB; // dimensionnement de la metrique if (!atdt) {if( !(unefois->CalResPrem_t)) { unefois->CalResPrem_t += 1; Tableau tab(10); 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 ; tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt; doCoHexa->metrique.PlusInitVariables(tab) ; };} else {if( !(unefois->CalResPrem_tdt)) {unefois->CalResPrem_tdt += 1; Tableau tab(11); tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt; tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_tdt;tab(11) = iVtdt; doCoHexa->metrique.PlusInitVariables(tab) ; };}; // initialisation du résidu residu->Zero(); ElemMeca::Cal_explicit ( doCoHexa->tab_ddl,d_epsBB, nombre->nbi,(doCoHexa->hexaed)->TaWi(),pa,atdt); return residu; }; // Calcul du residu local et de la raideur locale, // pour le schema implicite Element::ResRaid HexaMemb::Calcul_implicit (const ParaAlgoControle & pa) { HexaMemb::DonnComHexa* doCoHexa = unefois->doCoMemb; // pour simplifier l'écriture Tableau & d_epsBB = unefois->doCoMemb->d_epsBB; Tableau & d_sigHH = unefois->doCoMemb->d_sigHH; // dans le cas ou cas=true on enchaine les deux cas, sinon uniquement le second cas bool cald_Dvirtuelle = false; if (!(unefois->CalimpPrem)) { unefois->CalimpPrem += 1; Tableau tab(20); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0; tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt; tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt; tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ; tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt; doCoHexa->metrique.PlusInitVariables(tab) ; // on ne calcul la dérivée de la déformation virtuelle qu'une fois // car elle est constante dans le temps et indépendante des coordonnées cald_Dvirtuelle=true; }; // initialisation du résidu residu->Zero(); // initialisation de la raideur raideur->Zero(); ElemMeca::Cal_implicit (doCoHexa->tab_ddl, d_epsBB,(doCoHexa->d2_epsBB),d_sigHH, nombre->nbi,(doCoHexa->hexaed)->TaWi(),pa,cald_Dvirtuelle); Element::ResRaid el; el.res = residu; el.raid = raideur; return el; }; // Calcul de la matrice masse pour l'élément Mat_pleine * HexaMemb::CalculMatriceMasse (Enum_calcul_masse type_calcul_masse) { HexaMemb::DonnComHexa* doCoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionement de la métrique si nécessaire if (!(unefois->CalDynamique )) { unefois->CalDynamique += 1; Tableau tab(5); tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0; tab(4) = igijBB_tdt;tab(5) = igradVmoyBB_t; doCoHexa->metrique.PlusInitVariables(tab) ; // on vérifie le bon dimensionnement de la matrice if (type_calcul_masse == MASSE_CONSISTANTE) // dans le cas où la masse est consistante il faut la redimensionner { int nbddl = doCoHexa->tab_ddl.NbDdl(); (doCoHexa->matrice_masse).Initialise (nbddl,nbddl,0.); } }; // appel de la routine générale ElemMeca::Cal_Mat_masse (doCoHexa->tab_ddl,type_calcul_masse, nombre->nbiMas,(doCoHexa->hexaedMas)->TaPhi(),nombre->nbne ,(doCoHexa->hexaedMas)->TaWi()); return mat_masse; }; //============= lecture écriture dans base info ========== // cas donne le niveau de la récupération // = 1 : on récupère tout // = 2 : on récupère uniquement les données variables (supposées comme telles) void HexaMemb::Lecture_base_info (ifstream& ent,const Tableau * tabMaillageNoeud,const int cas) {// tout d'abord appel de la lecture de la classe elem_meca ElemMeca::Lecture_bas_inf(ent,tabMaillageNoeud,cas); // traitement du cas particulier de l'élément switch (cas) { case 1 : // ------- on récupère tout ------------------------- { // construction du tableau de ddl des noeuds ConstTabDdl(); // récup contraintes et déformation lesPtMecaInt.Lecture_base_info(ent,cas); break; } case 2 : // ----------- lecture uniquement de se qui varie -------------------- { // récup contraintes et déformation lesPtMecaInt.Lecture_base_info(ent,cas); break; } default : { cout << "\nErreur : valeur incorrecte du type de lecture !\n"; cout << "HexaMemb::Lecture_base_info(ifstream& ent,const " << "Tableau * tabMaillageNoeud,const int cas)" << " cas= " << cas << endl; Sortie(1); } } }; // cas donne le niveau de sauvegarde // = 1 : on sauvegarde tout // = 2 : on sauvegarde uniquement les données variables (supposées comme telles) void HexaMemb::Ecriture_base_info(ofstream& sort,const int cas) {// tout d'abord appel de l'écriture de la classe elem_meca ElemMeca::Ecriture_bas_inf(sort,cas); // traitement du cas particulier de l'élément switch (cas) { case 1 : // ------- on sauvegarde tout ------------------------- { // des tenseurs déformation et contrainte, lesPtMecaInt.Ecriture_base_info(sort,cas); break; } case 2 : // ----------- sauvegarde uniquement de se qui varie -------------------- { // des tenseurs déformation et contrainte, lesPtMecaInt.Ecriture_base_info(sort,cas); break; } default : { cout << "\nErreur : valeur incorrecte du type d'écriture !\n"; cout << "HexaMemb::Ecriture_base_info(ofstream& sort,const int cas)" << " cas= " << cas << endl; Sortie(1); } } }; // récupération des valeurs au numéro d'ordre = iteg pour les grandeurs enu // ici il s'agit de grandeurs scalaire, // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière Tableau HexaMemb::Valeur_a_diff_temps(bool absolue,Enum_dure enu_t,const List_io& enu,int iteg) { // appel de la procedure de elem meca int cas; if (!(unefois->dualSortHexa) && (unefois->CalimpPrem)) { cas=1;unefois->dualSortHexa += 1; } else if ((unefois->dualSortHexa) && (unefois->CalimpPrem)) { cas = 11;} else if ( !(unefois->dualSortHexa) && (unefois->CalResPrem_tdt)) { cas=2;unefois->dualSortHexa += 1; } else if ((unefois->dualSortHexa ) && (unefois->CalResPrem_tdt)) { cas = 12;} // sinon problème else { cout << "\n warning: les grandeurs ne sont pas calculees : il faudrait au moins un pas de calcul" << " pour inialiser les conteneurs des tenseurs resultats "; if (ParaGlob::NiveauImpression() >= 4) cout << "\n cas non prévu, unefois->dualSortHexa= " << unefois->dualSortHexa << " unefois->CalimpPrem= " << unefois->CalimpPrem << "\n HexaMemb::Valeur_a_diff_temps(Enum_dure enu_t..."; Sortie(1); }; return ElemMeca::Valeur_multi(absolue,enu_t,enu,iteg,cas); }; // récupération des valeurs au numéro d'ordre = iteg pour les grandeurs enu // ici il s'agit de grandeurs tensorielles, le retour s'effectue dans la liste // de conteneurs quelconque associée void HexaMemb::ValTensorielle_a_diff_temps(bool absolue,Enum_dure enu_t,List_io& enu,int iteg) { // appel de la procedure de elem meca int cas; if (!(unefois->dualSortHexa) && (unefois->CalimpPrem)) { cas=1;unefois->dualSortHexa += 1; } else if ((unefois->dualSortHexa) && (unefois->CalimpPrem)) { cas = 11;} else if (!(unefois->dualSortHexa) && (unefois->CalResPrem_tdt)) { cas=2;unefois->dualSortHexa += 1; } else if ((unefois->dualSortHexa) && (unefois->CalResPrem_tdt)) { cas = 12;} // sinon problème else { cout << "\n warning: les grandeurs ne sont pas calculees : il faudrait au moins un pas de calcul" << " pour inialiser les conteneurs des tenseurs resultats "; if (ParaGlob::NiveauImpression() >= 4) cout << "\n cas non prévu, unefois->dualSortHexa= " << unefois->dualSortHexa << " unefois->CalimpPrem= " << unefois->CalimpPrem << "\n HexaMemb::ValTensorielle_a_diff_temps(Enum_dure enu_t..."; Sortie(1); }; ElemMeca::Valeurs_Tensorielles(absolue,enu_t,enu,iteg,cas); }; // calcul du résidu et de la matrice de raideur pour le calcul d'erreur Element::Er_ResRaid HexaMemb::ContrainteAuNoeud_ResRaid() { HexaMemb::DonnComHexa* doCoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if (!(unefois->CalResPrem_t)) { unefois->CalResPrem_t += 1; Tableau tab(10); 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 ; tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt; doCoHexa->metrique.PlusInitVariables(tab) ; }; // appel du programme général ElemMeca:: SigmaAuNoeud_ResRaid(tab_noeud.Taille(), (doCoHexa->hexaed)->TaPhi(), (doCoHexa->hexaed)->TaWi(), doCoHexa-> resErr,doCoHexa->raidErr, (doCoHexa->hexaedEr)->TaPhi(), (doCoHexa->hexaedEr)->TaWi()); return (Element::Er_ResRaid( &(doCoHexa-> resErr),&(doCoHexa->raidErr))); } ; // calcul du résidu pour le calcul d'erreur ramenée aux noeuds Element::Er_ResRaid HexaMemb::ErreurAuNoeud_ResRaid() { HexaMemb::DonnComHexa* doCoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if ( !(unefois->CalResPrem_t) ) { unefois->CalResPrem_t += 1; Tableau tab(10); 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 ; tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt; doCoHexa->metrique.PlusInitVariables(tab) ; }; // appel du programme général ElemMeca::Cal_ErrAuxNoeuds(tab_noeud.Taille(), (doCoHexa->hexaed)->TaPhi(),(doCoHexa->hexaed)->TaWi() ,doCoHexa-> resErr ); return (Element::Er_ResRaid( &(doCoHexa-> resErr),&(doCoHexa->raidErr))); } ; // actualisation des ddl et des grandeurs actives de t+dt vers t void HexaMemb::TdtversT() { lesPtMecaInt.TdtversT(); // contrainte for (int ni=1;ni<= nombre->nbi; ni++) { if (tabSaveDon(ni) != NULL) tabSaveDon(ni)->TdtversT(); if (tabSaveTP(ni) != NULL) tabSaveTP(ni)->TdtversT(); if (tabSaveDefDon(ni) != NULL) tabSaveDefDon(ni)->TdtversT(); } ElemMeca::TdtversT_(); // appel de la procédure mère }; // actualisation des ddl et des grandeurs actives de t vers tdt void HexaMemb::TversTdt() { lesPtMecaInt.TversTdt(); // contrainte for (int ni=1;ni<= nombre->nbi; ni++) { if (tabSaveDon(ni) != NULL) tabSaveDon(ni)->TversTdt(); if (tabSaveTP(ni) != NULL) tabSaveTP(ni)->TversTdt(); if (tabSaveDefDon(ni) != NULL) tabSaveDefDon(ni)->TversTdt(); } ElemMeca::TversTdt_(); // appel de la procédure mère }; // calcul de l'erreur sur l'élément. void HexaMemb::ErreurElement(int type,double& errElemRelative ,double& numerateur, double& denominateur) { HexaMemb::DonnComHexa* doCoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if ( !(unefois->CalResPrem_t)) { unefois->CalResPrem_t += 1; Tableau tab(10); 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 ; tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt; doCoHexa->metrique.PlusInitVariables(tab) ; }; // on cré une instance de déformation // appel du programme général ElemMeca::Cal_ErrElem(type,errElemRelative,numerateur,denominateur, tab_noeud.Taille(),(doCoHexa->hexaed)->TaPhi(), (doCoHexa->hexaed)->TaWi(), (doCoHexa->hexaedEr)->TaPhi(),(doCoHexa->hexaedEr)->TaWi()); } ; // cas d'un chargement volumique, // force indique la force volumique appliquée // retourne le second membre résultant // ici on considère l'épaisseur de l'élément pour constituer le volume // -> explicite à t ou tdt en fonction du booleen atdt Vecteur HexaMemb::SM_charge_volumique_E (const Coordonnee& force,Fonction_nD* pt_fonct,bool atdt,const ParaAlgoControle & pa,bool sur_volume_finale_) { HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if (!atdt) {if ( !(unefois->CalSMvol_t)) { unefois->CalSMvol_t += 1; Tableau tab(10); 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 ; tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt; CoHexa->metrique.PlusInitVariables(tab) ; };} else {if (!(unefois->CalSMvol_tdt)) { unefois->CalSMvol_tdt += 1; Tableau tab(20); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0; tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt; tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt; tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ; tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt; CoHexa->metrique.PlusInitVariables(tab) ; }; }; // initialisation du résidu residu->Zero(); // appel du programme général d'elemmeca et retour du vecteur second membre // multiplié par l'épaisseur pour avoir le volume return ElemMeca::SM_charge_vol_E (CoHexa->tab_ddl,(CoHexa->hexaed)->TaPhi() ,tab_noeud.Taille(),(CoHexa->hexaed)->TaWi(),force,pt_fonct,pa,sur_volume_finale_,atdt); }; // calcul des seconds membres suivant les chargements // cas d'un chargement volumique, // force indique la force volumique appliquée // retourne le second membre résultant // ici on l'épaisseur de l'élément pour constituer le volume -> implicite, // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur // retourne le second membre et la matrice de raideur correspondant Element::ResRaid HexaMemb::SMR_charge_volumique_I (const Coordonnee& force,Fonction_nD* pt_fonct,const ParaAlgoControle & pa,bool sur_volume_finale_) { HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // initialisation du résidu residu->Zero(); // initialisation de la raideur raideur->Zero(); // -- définition des constantes de la métrique si nécessaire // en fait on fait appel aux même éléments que pour le calcul implicite if (!(unefois->CalSMvol_tdt)) { unefois->CalSMvol_tdt += 1; Tableau tab(20); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0; tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt; tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt; tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ; tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt; CoHexa->metrique.PlusInitVariables(tab) ; }; // appel du programme général d'elemmeca ElemMeca::SMR_charge_vol_I (CoHexa->tab_ddl ,(CoHexa->hexaed)->TaPhi(),tab_noeud.Taille() ,(CoHexa->hexaed)->TaWi(),force,pt_fonct,pa,sur_volume_finale_); Element::ResRaid el; el.res = residu; el.raid = raideur; return el; }; // calcul des seconds membres suivant les chargements // cas d'un chargement surfacique, sur les frontières des éléments // force indique la force surfacique appliquée // numface indique le numéro de la face chargée // retourne le second membre résultant // explicite à t et à tdt en fonction de la variable booléeenne atdt Vecteur HexaMemb::SM_charge_surfacique_E (const Coordonnee& force,Fonction_nD* pt_fonct,int numface,bool atdt,const ParaAlgoControle & pa) { HexaMemb::DonnComHexa* doCoHexa = unefois->doCoMemb; // pour simplifier l'écriture // initialisation du vecteur résidu ((*res_extS)(numface))->Zero(); // on récupère ou on crée la frontière surfacique Frontiere_surfacique(numface,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= tabb(numface)->Metrique(); // définition des constantes de la métrique si nécessaire if (!atdt) {if ( !(unefois->CalSMsurf_t )) { unefois->CalSMsurf_t += 1; Tableau tab(11); 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 ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_t; tab(11) = iVt; meta->PlusInitVariables(tab) ; };} else {if ( !(unefois->CalSMsurf_tdt)) { unefois->CalSMsurf_tdt += 1; Tableau tab(11); tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt; tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_tdt;tab(11) = iVtdt; meta->PlusInitVariables(tab) ; };}; // définition d'une déformation a doc if (defSurf(numface) == NULL) defSurf(numface) = new Deformation(*meta,tabb(numface)->TabNoeud(), (doCoHexa->quadraS).TaDphi(),(doCoHexa->quadraS).TaPhi()); // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SM_charge_surf_E (tabb(numface)->DdlElem(),numface ,(doCoHexa->quadraS).TaPhi(),(tabb(numface)->TabNoeud()).Taille() ,(doCoHexa->quadraS).TaWi(),force,pt_fonct,pa,atdt); }; // cas d'un chargement surfacique, sur les frontières des éléments // force indique la force surfacique appliquée // numface indique le numéro de la face chargée // retourne le second membre résultant // -> implicite, // pa : permet de déterminer si oui ou non on calcul la contribution à la raideur // retourne le second membre et la matrice de raideur correspondant Element::ResRaid HexaMemb::SMR_charge_surfacique_I (const Coordonnee& force,Fonction_nD* pt_fonct,int numface,const ParaAlgoControle & pa) { // initialisation du vecteur résidu et de la raideur ((*res_extS)(numface))->Zero(); ((*raid_extS)(numface))->Zero(); // on récupère ou on crée la frontière surfacique Frontiere_surfacique(numface,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= tabb(numface)->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique // définition des constantes de la métrique si nécessaire if ( !(unefois->CalSMRsurf)) { unefois->CalSMRsurf += 1; Tableau tab(20); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0; tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt; tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt; tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ; tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt; meta->PlusInitVariables(tab) ; }; // on définit la déformation a doc si elle n'existe pas déjà if (defSurf(numface) == NULL) defSurf(numface) = new Deformation(*meta,tabb(numface)->TabNoeud(), (CoHexa->quadraS).TaDphi(),(CoHexa->quadraS).TaPhi()); // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SMR_charge_surf_I (tabb(numface)->DdlElem(),numface ,(CoHexa->quadraS).TaPhi(),(tabb(numface)->TabNoeud()).Taille() ,(CoHexa->quadraS).TaWi(),force,pt_fonct,pa); }; // cas d'un chargement de type pression, sur les frontières des éléments // pression indique la pression appliquée // numface indique le numéro de la face chargée // retourne le second membre résultant // NB: il y a une définition par défaut pour les éléments qui n'ont pas de // surface externe -> message d'erreur d'où le virtuel et non virtuel pur // -> explicite à t ou à tdt en fonction de la variable booléenne atdt Vecteur HexaMemb::SM_charge_pression_E (double pression,Fonction_nD* pt_fonct,int numface,bool atdt,const ParaAlgoControle & pa) { // initialisation du vecteur résidu ((*res_extS)(numface))->Zero(); // on récupère ou on crée la frontière surfacique, si elle n'existe pas il faut la creer d'ou le true Frontiere_surfacique(numface,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= tabb(numface)->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // définition des constantes de la métrique si nécessaire if (!atdt) {if ( !(unefois->CalSMsurf_t)) { unefois->CalSMsurf_t += 1; Tableau tab(11); 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 ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_t; tab(11) = iVt; meta->PlusInitVariables(tab) ; };} else {if ( !(unefois->CalSMsurf_tdt)) { unefois->CalSMsurf_tdt += 1; Tableau tab(11); tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt; tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_tdt;tab(11) = iVtdt; meta->PlusInitVariables(tab) ; };}; // on définit la déformation a doc si elle n'existe pas déjà if (defSurf(numface) == NULL) defSurf(numface) = new Deformation (*meta,tabb(numface)->TabNoeud(),(CoHexa->quadraS).TaDphi(),(CoHexa->quadraS).TaPhi()); // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SM_charge_pres_E (tabb(numface)->DdlElem(),numface ,(CoHexa->quadraS).TaPhi(),(tabb(numface)->TabNoeud()).Taille() ,(CoHexa->quadraS).TaWi(),pression,pt_fonct,pa,atdt); }; // cas d'un chargement de type pression, sur les frontières des éléments // pression indique la pression appliquée // numface indique le numéro de la face chargée // retourne le second membre résultant // NB: il y a une définition par défaut pour les éléments qui n'ont pas de // surface externe -> message d'erreur d'où le virtuel et non virtuel pur // -> implicite, // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur // retourne le second membre et la matrice de raideur correspondant Element::ResRaid HexaMemb::SMR_charge_pression_I (double pression,Fonction_nD* pt_fonct,int numface,const ParaAlgoControle & pa) { // initialisation du vecteur résidu et de la raideur ((*res_extS)(numface))->Zero(); ((*raid_extS)(numface))->Zero(); // on récupère ou on crée la frontière surfacique, si elle n'existe pas il faut la creer d'ou le true Frontiere_surfacique(numface,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= tabb(numface)->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique // définition des constantes de la métrique si nécessaire if ( !(unefois->CalSMRsurf)) { unefois->CalSMRsurf += 1; Tableau tab(20); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0; tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt; tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt; tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ; tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt; meta->PlusInitVariables(tab) ; }; // on définit la déformation a doc si elle n'existe pas déjà if (defSurf(numface) == NULL) defSurf(numface) = new Deformation(*meta,tabb(numface)->TabNoeud(), (CoHexa->quadraS).TaDphi(),(CoHexa->quadraS).TaPhi()); // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SMR_charge_pres_I (tabb(numface)->DdlElem(),numface ,(CoHexa->quadraS).TaPhi(),(tabb(numface)->TabNoeud()).Taille() ,(CoHexa->quadraS).TaWi(),pression,pt_fonct,pa); }; // cas d'un chargement de type pression unidirectionnelle, sur les frontières des éléments // presUniDir indique le vecteur appliquée // numface indique le numéro de la face chargée // retourne le second membre résultant // -> explicite à t ou tdt en fonction de la variable booleenne atdt Vecteur HexaMemb::SM_charge_presUniDir_E (const Coordonnee& presUniDir,Fonction_nD* pt_fonct ,int numface,bool atdt,const ParaAlgoControle & pa) { // initialisation du vecteur résidu ((*res_extS)(numface))->Zero(); // on récupère ou on crée la frontière surfacique Frontiere_surfacique(numface,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= tabb(numface)->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // définition des constantes de la métrique si nécessaire if (!atdt) {if ( !(unefois->CalSMsurf_t )) { unefois->CalSMsurf_t += 1; Tableau tab(11); 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 ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_t; tab(11) = iVt; meta->PlusInitVariables(tab) ; };} else {if ( !(unefois->CalSMsurf_tdt )) { unefois->CalSMsurf_tdt += 1; Tableau tab(11); tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt; tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_tdt;tab(11) = iVtdt; meta->PlusInitVariables(tab) ; };}; // on définit la déformation a doc si elle n'existe pas déjà if (defSurf(numface) == NULL) defSurf(numface) = new Deformation (*meta,tabb(numface)->TabNoeud(),(CoHexa->quadraS).TaDphi(),(CoHexa->quadraS).TaPhi()); // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SM_charge_surf_Suiv_E (tabb(numface)->DdlElem(),numface ,(CoHexa->quadraS).TaPhi(),(tabb(numface)->TabNoeud()).Taille() ,(CoHexa->quadraS).TaWi(),presUniDir,pt_fonct,pa,atdt); }; // -> implicite, // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur // retourne le second membre et la matrice de raideur correspondant Element::ResRaid HexaMemb::SMR_charge_presUniDir_I (const Coordonnee& presUniDir,Fonction_nD* pt_fonct,int numface,const ParaAlgoControle & pa) { // initialisation du vecteur résidu et de la raideur ((*res_extS)(numface))->Zero(); ((*raid_extS)(numface))->Zero(); // on récupère ou on crée la frontière surfacique Frontiere_surfacique(numface,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= tabb(numface)->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique // définition des constantes de la métrique si nécessaire if ( !(unefois->CalSMRsurf)) { unefois->CalSMRsurf += 1; Tableau tab(20); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0; tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt; tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt; tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ; tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt; meta->PlusInitVariables(tab) ; }; // on définit la déformation a doc si elle n'existe pas déjà if (defSurf(numface) == NULL) defSurf(numface) = new Deformation(*meta,tabb(numface)->TabNoeud(), (CoHexa->quadraS).TaDphi(),(CoHexa->quadraS).TaPhi()); // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SMR_charge_surf_Suiv_I (tabb(numface)->DdlElem(),numface ,(CoHexa->quadraS).TaPhi(),(tabb(numface)->TabNoeud()).Taille() ,(CoHexa->quadraS).TaWi(),presUniDir,pt_fonct,pa); }; // cas d'un chargement lineique, sur les arêtes frontières des éléments // force indique la force lineique appliquée // numArete indique le numéro de l'arête chargée // retourne le second membre résultant // -> explicite à t ou tdt en fonction de la variable booléenne atdt Vecteur HexaMemb::SM_charge_lineique_E (const Coordonnee& force,Fonction_nD* pt_fonct ,int numArete,bool atdt,const ParaAlgoControle & pa) { // initialisation du vecteur résidu ((*res_extA)(numArete))->Zero(); // on récupère ou on crée la frontière arrête ElFrontiere* elf = Frontiere_lineique(numArete,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= elf->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if (!atdt) {if( !(unefois->CalSMlin_t)) { unefois->CalSMlin_t += 1; Tableau tab(8); 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 ; tab(8) = igradVBB_t; meta->PlusInitVariables(tab) ; };} else {if( !(unefois->CalSMlin_tdt)) { unefois->CalSMlin_tdt += 1; Tableau tab(8); tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt; tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ; tab(8) = igradVBB_tdt; meta->PlusInitVariables(tab) ; };}; // on définit la déformation a doc si elle n'existe pas déjà if (defArete(numArete) == NULL) defArete(numArete) = new Deformation(*meta,elf->TabNoeud(), (CoHexa->segS).TaDphi(),(CoHexa->segS).TaPhi()); // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SM_charge_line_E (elf->DdlElem(),numArete ,(CoHexa->segS).TaPhi(),(elf->TabNoeud()).Taille() ,(CoHexa->segS).TaWi(),force,pt_fonct,pa,atdt); }; // cas d'un chargement lineique, sur les aretes frontières des éléments // force indique la force lineique appliquée // numarete indique le numéro de l'arete chargée // retourne le second membre résultant // NB: il y a une définition par défaut pour les éléments qui n'ont pas // d'arete externe -> message d'erreur d'où le virtuel et non virtuel pur // -> implicite, // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur // retourne le second membre et la matrice de raideur correspondant Element::ResRaid HexaMemb::SMR_charge_lineique_I (const Coordonnee& force,Fonction_nD* pt_fonct,int numArete,const ParaAlgoControle & pa) { // initialisation du vecteur résidu et de la raideur ((*res_extA)(numArete))->Zero(); ((*raid_extA)(numArete))->Zero(); // on récupère ou on crée la frontière arrête ElFrontiere* elf = Frontiere_lineique(numArete,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les triangles // du même type Met_abstraite * meta= elf->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if( !(unefois->CalSMRlin)) { unefois->CalSMRlin += 1; Tableau tab(15); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0; tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt; tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt; tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; tab(15) = igradVBB_tdt; meta->PlusInitVariables(tab) ; }; // on définit la déformation a doc si elle n'existe pas déjà if (defArete(numArete) == NULL) defArete(numArete) = new Deformation(*meta,elf->TabNoeud(), (CoHexa->segS).TaDphi(),(CoHexa->segS).TaPhi()); // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SMR_charge_line_I (elf->DdlElem(),numArete ,(CoHexa->segS).TaPhi(),(elf->TabNoeud()).Taille() ,(CoHexa->segS).TaWi(),force,pt_fonct,pa); }; // cas d'un chargement de type pression hydrostatique, sur les frontières des éléments // la charge dépend de la hauteur à la surface libre du liquide déterminée par un point // et une direction normale à la surface libre: // nSurf : le numéro de la surface externe // poidvol: indique le poids volumique du liquide // M_liquide : un point de la surface libre // dir_normal_liquide : direction normale à la surface libre // retourne le second membre résultant // calcul à l'instant tdt ou t en fonction de la variable atdt // retourne le second membre résultant // NB: il y a une définition par défaut pour les éléments qui n'ont pas de // surface externe -> message d'erreur d'où le virtuel et non virtuel pur // -> explicite à t ou à tdt en fonction de la variable booléenne atdt Vecteur HexaMemb::SM_charge_hydrostatique_E(const Coordonnee& dir_normal_liquide,const double& poidvol ,int numface,const Coordonnee& M_liquide,bool atdt ,const ParaAlgoControle & pa ,bool sans_limitation) { // initialisation du vecteur résidu ((*res_extS)(numface))->Zero(); // on récupère ou on crée la frontière surfacique Frontiere_surfacique(numface,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= tabb(numface)->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // définition des constantes de la métrique si nécessaire if (!atdt) {if ( !(unefois->CalSMsurf_t)) { unefois->CalSMsurf_t += 1; Tableau tab(11); 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 ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_t; tab(11) = iVt; meta->PlusInitVariables(tab) ; };} else {if ( !(unefois->CalSMsurf_tdt )) { unefois->CalSMsurf_tdt += 1; Tableau tab(11); tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt; tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_tdt;tab(11) = iVtdt; meta->PlusInitVariables(tab) ; };}; // on définit la déformation a doc si elle n'existe pas déjà if (defSurf(numface) == NULL) defSurf(numface) = new Deformation (*meta,tabb(numface)->TabNoeud(),(CoHexa->quadraS).TaDphi(),(CoHexa->quadraS).TaPhi()); // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SM_charge_hydro_E (tabb(numface)->DdlElem(),numface ,(CoHexa->quadraS).TaPhi(),(tabb(numface)->TabNoeud()).Taille() ,(CoHexa->quadraS).TaWi(),dir_normal_liquide,poidvol,M_liquide,sans_limitation,pa,atdt); }; // cas d'un chargement de type pression hydrostatique, sur les frontières des éléments // la charge dépend de la hauteur à la surface libre du liquide déterminée par un point // et une direction normale à la surface libre: // nSurf : le numéro de la surface externe // poidvol: indique le poids volumique du liquide // M_liquide : un point de la surface libre // dir_normal_liquide : direction normale à la surface libre // retourne le second membre résultant // calcul à l'instant tdt ou t en fonction de la variable atdt // retourne le second membre résultant // NB: il y a une définition par défaut pour les éléments qui n'ont pas de // surface externe -> message d'erreur d'où le virtuel et non virtuel pur // -> implicite, // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur // retourne le second membre et la matrice de raideur correspondant Element::ResRaid HexaMemb::SMR_charge_hydrostatique_I (const Coordonnee& dir_normal_liquide,const double& poidvol ,int numface,const Coordonnee& M_liquide ,const ParaAlgoControle & pa ,bool sans_limitation) { // initialisation du vecteur résidu et de la raideur ((*res_extS)(numface))->Zero(); ((*raid_extS)(numface))->Zero(); // on récupère ou on crée la frontière surfacique Frontiere_surfacique(numface,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= tabb(numface)->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique // définition des constantes de la métrique si nécessaire if ( !(unefois->CalSMRsurf)) { unefois->CalSMRsurf += 1; Tableau tab(20); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0; tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt; tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt; tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ; tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt; meta->PlusInitVariables(tab) ; }; // on définit la déformation a doc si elle n'existe pas déjà if (defSurf(numface) == NULL) defSurf(numface) = new Deformation(*meta,tabb(numface)->TabNoeud(), (CoHexa->quadraS).TaDphi(),(CoHexa->quadraS).TaPhi()); // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SMR_charge_hydro_I (tabb(numface)->DdlElem(),numface ,(CoHexa->quadraS).TaPhi(),(tabb(numface)->TabNoeud()).Taille() ,(CoHexa->quadraS).TaWi(),dir_normal_liquide,poidvol,M_liquide,sans_limitation,pa); }; // cas d'un chargement surfacique hydro-dynamique, // Il y a trois forces: une suivant la direction de la vitesse: de type traînée aerodynamique // Fn = poids_volu * fn(V) * S * (normale*u) * u, u étant le vecteur directeur de V (donc unitaire) // une suivant la direction normale à la vitesse de type portance // Ft = poids_volu * ft(V) * S * (normale*u) * w, w unitaire, normal à V, et dans le plan n et V // une suivant la vitesse tangente de type frottement visqueux // T = to(Vt) * S * ut, Vt étant la vitesse tangentielle et ut étant le vecteur directeur de Vt // coef_mul: est un coefficient multiplicateur global (de tout) // retourne le second membre résultant // // -> explicite à t ou à tdt en fonction de la variable booléenne atdt Vecteur HexaMemb::SM_charge_hydrodynamique_E( Courbe1D* frot_fluid,const double& poidvol , Courbe1D* coef_aero_n,int numface,const double& coef_mul , Courbe1D* coef_aero_t,bool atdt ,const ParaAlgoControle & pa) { // initialisation du vecteur résidu ((*res_extS)(numface))->Zero(); // on récupère ou on crée la frontière surfacique Frontiere_surfacique(numface,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= tabb(numface)->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // définition des constantes de la métrique si nécessaire if (!atdt) {if ( !(unefois->CalSMsurf_t )) { unefois->CalSMsurf_t += 1; Tableau tab(11); 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 ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_t; tab(11) = iVt; meta->PlusInitVariables(tab) ; };} else {if ( !(unefois->CalSMsurf_tdt )) { unefois->CalSMsurf_tdt += 1; Tableau tab(11); tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt; tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ; tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ; tab(10) = igradVBB_tdt;tab(11) = iVtdt; meta->PlusInitVariables(tab) ; };}; // on définit la déformation a doc si elle n'existe pas déjà if (defSurf(numface) == NULL) defSurf(numface) = new Deformation (*meta,tabb(numface)->TabNoeud(),(CoHexa->quadraS).TaDphi(),(CoHexa->quadraS).TaPhi()); // appel du programme général d'ElemMeca et retour du vecteur second membre return ElemMeca::SM_charge_hydrodyn_E (poidvol,(CoHexa->quadraS).TaPhi(),(tabb(numface)->TabNoeud()).Taille() ,frot_fluid,(CoHexa->quadraS).TaWi() ,coef_aero_n,numface,coef_mul,coef_aero_t,pa,atdt); }; // -> implicite, // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur // retourne le second membre et la matrice de raideur correspondant Element::ResRaid HexaMemb::SMR_charge_hydrodynamique_I( Courbe1D* frot_fluid,const double& poidvol , Courbe1D* coef_aero_n,int numface,const double& coef_mul , Courbe1D* coef_aero_t ,const ParaAlgoControle & pa) { // initialisation du vecteur résidu et de la raideur ((*res_extS)(numface))->Zero(); ((*raid_extS)(numface))->Zero(); // on récupère ou on crée la frontière surfacique Frontiere_surfacique(numface,true); // on utilise la métrique des éléments de frontière // avec l'instance déformation dédiée pour // récupération de la métrique associée à l'élément qui est commune a tous les quadrangles // du même type Met_abstraite * meta= tabb(numface)->Metrique(); HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique // définition des constantes de la métrique si nécessaire if ( !(unefois->CalSMRsurf)) { unefois->CalSMRsurf += 1; Tableau tab(20); tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0; tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt; tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt; tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ; tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt; meta->PlusInitVariables(tab) ; }; // on définit la déformation a doc si elle n'existe pas déjà if (defSurf(numface) == NULL) defSurf(numface) = new Deformation(*meta,tabb(numface)->TabNoeud(), (CoHexa->quadraS).TaDphi(),(CoHexa->quadraS).TaPhi()); // appel du programme général d'ElemMeca et retour du vecteur second membre return ElemMeca::SM_charge_hydrodyn_I (poidvol,(CoHexa->quadraS).TaPhi(),(tabb(numface)->TabNoeud()).Taille() ,frot_fluid,(CoHexa->quadraS).TaWi(),tabb(numface)->DdlElem() ,coef_aero_n,numface,coef_mul ,coef_aero_t,pa); }; // Calcul des frontieres surfaces de l'element // creation des elements frontieres et retour du tableau de ces elements // la création n'a lieu qu'au premier appel // ou lorsque l'on force le paramètre force a true // dans ce dernier cas seul les frontière effacées sont recréée Tableau const & HexaMemb::Frontiere(bool force) { int cas = 1; // ici on ne crée que les frontières faces return Frontiere_elemeca(cas,force); // // le calcul et la création ne sont effectués qu'au premier appel // // ou lorsque l'on veut forcer une recréation // if (((ind_front_lin == 0) && (ind_front_surf == 0) && (ind_front_point == 0)) // || force ) //// if ( (ind_front_surf==0)|| force ||(ind_front_surf==2) ) // { HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // ElemGeomC0 & el = *(CoHexa->hexaed); // DdlElement & tdd = CoHexa->tab_ddl; // int tail_s = el.Nonf().Taille(); // nombre de faces // // def de la taille // if (((ind_front_point > 0) || (ind_front_lin > 0)) && (ind_front_surf == 0)) // // cas où les frontières points ou linéiques existent mais pas les frontières faces // // on redimentionne et on transfert // { int tail_af = nombre->nbne+el.NonS().Taille(); // nombre d'arête + noeud // int taille_f = tail_af + tail_s; // tabb.Change_taille(taille_f); // for (int i=1;i<= tail_af;i++) // { tabb(i+tail_s) = tabb(i); // tabb(i) = NULL; // } // posi_tab_front_lin += tail_s; // posi_tab_front_point += tail_s; // } // else if (ind_front_surf == 0) // // cas où il n'y a ni points ni lignes de frontières, redimentionnement éventuelle // tabb.Change_taille(tail_s); // // création éventuelle des faces // for (int nface=1;nface<=tail_s;nface++) // { int nbnoe = el.Nonf()(nface).Taille(); // nb noeud de la face // Tableau tab(nbnoe); // les noeuds des faces frontieres // DdlElement ddelem(nombre->nbneS); // les ddlelements des noeuds frontieres // for (int i=1;i<= nbnoe;i++) // { int nib=el.Nonf()(nface)(i); // numéro dans l'élément du numéro i de la face // tab(i) = tab_noeud(nib); // DdlNoeudElement einter = tdd(nib); // ddelem.Change_un_ddlNoeudElement(i,einter); // // ddelem(i) = tdd(el.Nonf()(nface)(i)); // } // if (tabb(nface) == NULL) // tabb(nface) = new_frontiere_surf(tab,ddelem); // } // // mise à jour de l'indicateur de création de faces // ind_front_surf = 1; // } // return (Tableau &)tabb; }; //// ramène la frontière point //// éventuellement création des frontieres points de l'element et stockage dans l'element //// si c'est la première fois sinon il y a seulement retour de l'elements //// a moins que le paramètre force est mis a true //// dans ce dernier cas la frontière effacéee est recréée //// num indique le numéro du point à créer (numérotation EF) //ElFrontiere* const HexaMemb::Frontiere_points(int num,bool force) // { // le calcul et la création ne sont effectués qu'au premier appel // // ou lorsque l'on veut forcer une recréation // #ifdef MISE_AU_POINT // if ((num > nombre->nbne)||(num <=0)) // { cout << "\n *** erreur, le noeud demande pour frontiere: " << num << " esten dehors de la numeration de l'element ! " // << "\n Frontiere_points(int num,bool force)"; // Sortie(1); // } // #endif // // if ((ind_front_point == 0) || force || (ind_front_point == 2)) // {Tableau tab(1); // les noeuds des points frontieres // DdlElement ddelem(1); // les ddlelements des points frontieres // // on regarde si les frontières points existent sinon on les crée // if (ind_front_point == 1) // return (ElFrontiere*)tabb(posi_tab_front_point+num); // else if ( ind_front_point == 2) // // cas où certaines frontières existent // if (tabb(posi_tab_front_point+num) != NULL) // return (ElFrontiere*)tabb(posi_tab_front_point+num); // // dans tous les autres cas on construit la frontière point // // on commence par dimensionner le tableau de frontière, comme les frontières points sont // // les dernières, il suffit de les ajouter, d'où on redimentionne le tableau mais on ne créra // // que la frontière adoc // // def de la taille // int taille_actuelle = tabb.Taille(); // if ((ind_front_point == 0) && ((ind_front_lin > 0) || (ind_front_surf > 0))) // // cas où les frontières lignes ou surfaces existent, mais pas les points // { int tail_p = nombre->nbne; // le nombre de noeuds // int taille_f = taille_actuelle + tail_p; // tabb.Change_taille(taille_f); // for (int i=1;i<= tail_p;i++) // { tabb(i+taille_actuelle) = tabb(i);tabb(i) = NULL;}; // posi_tab_front_point +=taille_actuelle; // if (ind_front_lin > 0) posi_tab_front_lin += taille_actuelle; // } // else if (ind_front_point == 0) // // cas où aucune frontière n'existe // {tabb.Change_taille(nombre->nbne);}; // // dans les autres cas, les frontières points exitent donc pas à dimensionner // // on définit tous les points par simplicité // for (int i=1;i<=nombre->nbne;i++) // {tab(1) = tab_noeud(i);ddelem.Change_un_ddlNoeudElement(1,unefois->doCoMemb->tab_ddl(i)); // if (tabb(i+posi_tab_front_point) == NULL) // tabb(i+posi_tab_front_point) = new FrontPointF (tab,ddelem); // }; // }; // return (ElFrontiere*)tabb(num+posi_tab_front_point); // }; //// ramène la frontière linéique //// éventuellement création des frontieres linéique de l'element et stockage dans l'element //// si c'est la première fois et en 3D sinon il y a seulement retour de l'elements //// a moins que le paramètre force est mis a true //// dans ce dernier cas la frontière effacéee est recréée //// num indique le numéro de l'arête à créer (numérotation EF) //ElFrontiere* const HexaMemb::Frontiere_lineique(int num,bool force) // { // le calcul et la création ne sont effectués qu'au premier appel // // ou lorsque l'on veut forcer une recréation // if ((ind_front_lin == 0) || force || (ind_front_lin == 2)) // { HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture // ElemGeomC0 & el = *(CoHexa->hexaed); // DdlElement & tdd = CoHexa->tab_ddl; // int taille = tabb.Taille(); // la taille initiales des frontières // int tail_s = el.Nonf().Taille(); // nombre de faces // int tail_a = el.NonS().Taille(); // nombre d'arête // posi_tab_front_lin = 0; // init indice de début d'arrête dans tabb // // dimensionnement du tableau de frontières ligne si nécessaire // if (ind_front_lin == 0) // {if (ind_front_surf > 0) // // cas où il y a déjà des surfaces et pas de lignes // { tabb.Change_taille(tail_a + tail_s); // posi_tab_front_lin = tail_s; // } // else // cas pas de lignes et de surfaces // { tabb.Change_taille(tail_a); // // et c'est ok pour posi_tab_front_lin // } // } // // création éventuelle des lignes // for (int nlign=1;nlign<=tail_a;nlign++) // if (tabb(posi_tab_front_lin+nlign) == NULL) // { int nbnoe = el.NonS()(nlign).Taille(); // nb noeud de l'arête // Tableau tab(nbnoe); // les noeuds de l'arête frontiere // DdlElement ddelem(nombre->nbneA); // les ddlelements des noeuds frontieres // for (int i=1;i<= nbnoe;i++) // { tab(i) = tab_noeud(el.NonS()(nlign)(i)); // ddelem.Change_un_ddlNoeudElement(i,tdd(el.NonS()(nlign)(i))); // // ddelem(i) = tdd(el.NonS()(nlign)(i)); // } // tabb(posi_tab_front_lin+nlign) = new_frontiere_lin(tab,ddelem); // } // // mise à jour de l'indicateur // ind_front_lin = 1; // } // // return (ElFrontiere*)tabb(num + posi_tab_front_lin); // }; //// ramène la frontière surfacique //// éventuellement création d'une frontieres surfacique de l'element et stockage dans l'element //// si c'est la première fois sinon il y a seulement retour de l'elements //// a moins que le paramètre force est mis a true //// dans ce dernier cas la frontière effacéee est recréée //// num indique le numéro de la surface à créer (numérotation EF) //ElFrontiere* const HexaMemb::Frontiere_surfacique(int num,bool force) // { // le calcul et la création ne sont effectués qu'au premier appel // // ou lorsque l'on veut forcer une recréation // // on regarde si les frontières surfaciques existent sinon on les crée // if (ind_front_surf == 1) // return (ElFrontiere*)tabb(num); // else if ( ind_front_surf == 2) // // cas où certaines frontières existent // {if (tabb(num) != NULL) // return (ElFrontiere*)tabb(num); // }; // // dans tous les autres cas on reconstruit les frontières // Frontiere(force); // return (ElFrontiere*)tabb(num); // }; // liberation de la place pointee void HexaMemb::Libere () {Element::Libere (); // liberation de residu et raideur LibereTenseur() ; // liberation des tenseur intermediaires }; // acquisition ou modification d'une loi de comportement void HexaMemb::DefLoi (LoiAbstraiteGeneral * NouvelleLoi) { // verification du type de loi if (NouvelleLoi == NULL) { cout << "\n Erreur, la loi de comportement n'est pas definit " << " on ne peut pas continuer !! \n" ; Element::Affiche(1); Sortie(1); }; if ((NouvelleLoi->Dimension_loi() != 3) && (NouvelleLoi->Dimension_loi() != 4)) { cout << "\n Erreur, la loi de comportement a utiliser avec des HexaMembs"; cout << " doit etre de type 3D, \n ici elle est de type = " << (NouvelleLoi->Dimension_loi()) << "D !!! " << endl; Sortie(1); }; // cas d'une loi mécanique if (GroupeMecanique(NouvelleLoi->Id_categorie())) {loiComp = (Loi_comp_abstraite *) NouvelleLoi; // initialisation du stockage particulier, ici en fonction du nb de pt d'integ int imax = tabSaveDon.Taille(); for (int i=1;i<= imax;i++) tabSaveDon(i) = loiComp->New_et_Initialise(); // idem pour le type de déformation mécanique associé int iDefmax = tabSaveDefDon.Taille(); for (int i=1;i<= iDefmax;i++) tabSaveDefDon(i) = def->New_et_Initialise(); // définition du type de déformation associé à la loi loiComp->Def_type_deformation(*def); // on active les données particulières nécessaires au fonctionnement de la loi de comp loiComp->Activation_donnees(tab_noeud,dilatation,lesPtMecaInt); }; // cas d'une loi thermo physique if (GroupeThermique(NouvelleLoi->Id_categorie())) {loiTP = (CompThermoPhysiqueAbstraite *) NouvelleLoi; // initialisation du stockage particulier thermo physique, int imax = tabSaveTP.Taille(); for (int i=1;i<= imax;i++) tabSaveTP(i) = loiTP->New_et_Initialise(); // on active les données particulières nécessaires au fonctionnement de la loi de comp loiTP->Activation_donnees(tab_noeud); }; // cas d'une loi de frottement if (GroupeFrottement(NouvelleLoi->Id_categorie())) loiFrot = (CompFrotAbstraite *) NouvelleLoi; }; // test si l'element est complet int HexaMemb::TestComplet() { int res = ElemMeca::TestComplet(); // test dans la fonction mere if ( tab_noeud(1) == NULL) { cout << "\n les noeuds de l\'element hexaedrique ne sont pas defini \n"; res = 0; } else { int testi =1; int posi = Id_nom_ddl("X1") -1; int dim = ParaGlob::Dimension(); int jmax = tab_noeud.Taille(); for (int i =1; i<= dim; i++) for (int j=1;j<=jmax;j++) if(!(tab_noeud(j)->Existe_ici(Enum_ddl(posi+i)))) testi = 0; if(testi == 0) { cout << "\n les ddls Xi des noeuds ne sont pas defini \n"; cout << " \n utilisez HexaMemb::ConstTabDdl() pour completer " ; res = 0; }; }; return res; }; // procedure permettant de completer l'element apres // sa creation avec les donnees du bloc transmis // ici pour l'instant rien de spécifique Element* HexaMemb::Complete(BlocGen & bloc,LesFonctions_nD* lesFonctionsnD) { // complétion avec bloc ElemMeca::Complete_ElemMeca(bloc,lesFonctionsnD); return this; }; // Compléter pour la mise en place de la gestion de l'hourglass Element* HexaMemb::Complet_Hourglass(LoiAbstraiteGeneral * loiHourglass, const BlocGen & bloc) { // on initialise le traitement de l'hourglass string str_precision; // string vide indique que l'on veut utiliser un élément normal // dans le cas où il s'agit d'élément Quadratique il faut utiliser 27 pti par défaut if ((id_interpol == QUADRATIQUE) || (id_interpol == QUADRACOMPL)) str_precision = "_cm27pti"; ElemMeca::Init_hourglass_comp(*(unefois->doCoMemb->hexaedHourg),str_precision,loiHourglass,bloc); // dans le cas où l'hourglass a été activé mais que l'élément n'a pas // de traitement particulier associé, alors on désactive l'hourglass if ( ((type_stabHourglass == STABHOURGLASS_PAR_COMPORTEMENT) || (type_stabHourglass == STABHOURGLASS_PAR_COMPORTEMENT_REDUIT)) &&(unefois->doCoMemb->hexaedHourg == NULL)) type_stabHourglass = STABHOURGLASS_NON_DEFINIE; return this; }; // ajout du tableau de ddl des noeuds void HexaMemb::ConstTabDdl() { Tableau ta(ParaGlob::Dimension()); int posi = Id_nom_ddl("X1") -1; int dim = ParaGlob::Dimension(); for (int i =1; i<= dim; i++) {Ddl inter((Enum_ddl(i+posi)),0.,LIBRE); ta(i) = inter; }; // attribution des ddls aux noeuds for (int ne=1; ne<= nombre->nbne; ne++) tab_noeud(ne)->PlusTabDdl(ta); }; // =====>>>> methodes privées appelees par les classes dérivees <<<<===== // fonction d'initialisation servant dans les classes derivant // au niveau du constructeur // les pointeurs d'éléments sont non nul uniquement lorsque doCoMemb est null // c'est-à-dire pour l'initialisation HexaMemb::DonnComHexa* HexaMemb::Init (ElemGeomC0* hexa,ElemGeomC0* hexaEr,ElemGeomC0* hexaMas ,ElemGeomC0* hexaeHourg,bool sans_init_noeud) { id_geom=HEXAEDRE; // // le fait de mettre les pointeurs a null permet // de savoir que l'element n'est pas complet // dans le cas d'un constructeur avec tableau de noeud, il ne faut pas mettre // les pointeurs à nuls d'où le test if (!sans_init_noeud) for (int i =1;i<= nombre->nbne;i++) tab_noeud(i) = NULL; // definition des donnees communes aux HexaMembxxx // a la premiere definition d'une instance if (unefois->doCoMemb == NULL) unefois->doCoMemb = HexaMemb::Def_DonneeCommune(hexa,hexaEr,hexaMas,hexaeHourg); HexaMemb::DonnComHexa* doCoHexa = unefois->doCoMemb; // pour simplifier l'écriture met = &(doCoHexa->metrique); // met est defini dans elemeca // def pointe sur la deformation specifique a l'element def = new Deformation(*met,tab_noeud,(doCoHexa->hexaed)->TaDphi(),(doCoHexa->hexaed)->TaPhi()); // idem pour la remonte aux contraintes et le calcul d'erreur defEr = new Deformation(*met,tab_noeud,(doCoHexa->hexaedEr)->TaDphi(),(doCoHexa->hexaedEr)->TaPhi()); // idem pour le calcul de la masse defMas = new Deformation(*met,tab_noeud,(doCoHexa->hexaedMas)->TaDphi(),(doCoHexa->hexaedMas)->TaPhi()); // idem pour le calcul de second membre int nbFace = (doCoHexa->hexaed)->NbFe(); defSurf.Change_taille(nbFace); // une seule surface pour le second membre for (int ii=1;ii<=nbFace;ii++) defSurf(ii) = NULL; // la déformation sera construite si nécessaire au moment du calcul de // second membre // idem pour le calcul de second membre int nbArete = (doCoHexa->hexaed)->NbSe(); defArete.Change_taille(nbArete); // 12 arrêtes utilisées pour le second membre // la déformation sera construite si nécessaire au moment du calcul de second membre for (int ia=1;ia<=nbArete;ia++) defArete(ia) = NULL; //dimensionnement des deformations et contraintes etc.. int dimtens = 3; lesPtMecaInt.Change_taille_PtIntegMeca(nombre->nbi,dimtens); // attribution des numéros de référencement dans le conteneur for (int ni = 1; ni<= nombre->nbi; ni++) {lesPtMecaInt(ni).Change_Nb_mail(this->num_maillage); lesPtMecaInt(ni).Change_Nb_ele(this->num_elt); lesPtMecaInt(ni).Change_Nb_pti(ni); }; // stockage des donnees particulieres de la loi de comportement mécanique au point d'integ tabSaveDon.Change_taille(nombre->nbi); // stockage des donnees particulieres de la loi de comportement thermo physique au point d'integ tabSaveTP.Change_taille(nombre->nbi); // stockage des donnees particulieres de la déformation mécanique au point d'integ tabSaveDefDon.Change_taille(nombre->nbi); tab_energ.Change_taille(nombre->nbi); tab_energ_t.Change_taille(nombre->nbi); // initialisation des pointeurs définis dans la classe Element concernant les résidus et // raideur // --- cas de la puissance interne --- residu = &(doCoHexa->residu_interne); // residu local raideur = &(doCoHexa->raideur_interne); // raideur locale // --- cas de la dynamique ----- mat_masse = &(doCoHexa->matrice_masse); // --- cas des efforts externes concernant les points ------ res_extN = &(doCoHexa->residus_externeN); // pour les résidus et second membres raid_extN= &(doCoHexa->raideurs_externeN);// pour les raideurs // --- cas des efforts externes concernant les aretes ------ res_extA = &(doCoHexa->residus_externeA); // pour les résidus et second membres raid_extA= &(doCoHexa->raideurs_externeA);// pour les raideurs // --- cas des efforts externes concernant les faces ------ res_extS= &(doCoHexa->residus_externeS); // pour les résidus et second membres raid_extS= &(doCoHexa->raideurs_externeS); // pour les raideurs return doCoHexa; }; // fonction permettant le calcul de doCoHexa HexaMemb::DonnComHexa* HexaMemb::Def_DonneeCommune (ElemGeomC0* hexa,ElemGeomC0* hexaEr,ElemGeomC0* hexaMas, ElemGeomC0* hexaeHourg) { // degre de liberte int dim = ParaGlob::Dimension(); DdlElement tab_ddl(nombre->nbne,dim); int posi = Id_nom_ddl("X1") -1; for (int i =1; i<= dim; i++) for (int j=1; j<= nombre->nbne; j++) // tab_ddl (j,i) = Enum_ddl(i+posi); tab_ddl.Change_Enum(j,i,Enum_ddl(i+posi)); // cas des ddl éléments secondaires pour le calcul d'erreur // def du nombre de composantes du tenseur de contrainte en absolu // en fait 6 ici int nbcomposante = 6; DdlElement tab_ddlErr(nombre->nbne,nbcomposante); posi = Id_nom_ddl("SIG11") -1; // uniquement un cas est considéré 6 composantes for (int j=1; j<= nombre->nbne; j++) { // on definit le nombre de composante de sigma en absolu // les composantes sont a suivre dans l'enumération for (int i= 1;i<= nbcomposante; i++) tab_ddlErr.Change_Enum(j,i,Enum_ddl(i+posi)); }; // egalement pour tab_Err1Sig11, def d'un tableau de un ddl : enum SIG11 // par noeud DdlElement tab_Err1Sig11(nombre->nbne,DdlNoeudElement(SIG11)); // pour le calcul des seconds membres, def de l'élément géométrique de facette GeomQuadrangle quadraS(nombre->nbiS,nombre->nbneS); GeomSeg seggS(nombre->nbiA,nombre->nbneA); // def metrique // on definit les variables a priori toujours utiles Tableau tab(15); tab(1) = iM0; tab(2) = iMt; tab(3) = iMtdt ; tab(4)=igiB_0;tab(5)=igiB_t;tab(6)=igiB_tdt; tab(7)=igiH_0;tab(8)=igiH_t;tab(9)=igiH_tdt ; tab(10)=igijBB_0;tab(11)=igijBB_t;tab(12)=igijBB_tdt; tab(13)=igijHH_0;tab(14)=igijHH_t;tab(15)=igijHH_tdt ; // dim du pb , nb de vecteur de la base = 3 ici, tableau de ddl et la def de variables Met_abstraite metri(dim,3,tab_ddl,tab,nombre->nbne) ; // ---- cas du calcul d'erreur sur sigma ou epsilon Tableau resEr(nbcomposante); for (int i=1; i<= nbcomposante; i++) resEr(i)=new Vecteur (nombre->nbne); Mat_pleine raidEr(nombre->nbne,nombre->nbne); // la raideur pour l'erreur // dimensionnement des différents résidus et raideurs pour le calcul mécanique int nbddl = tab_ddl.NbDdl(); int nbA = hexa->NbSe(); // nombre d'arêtes int nbS = hexa->NbFe(); // nombre de faces Vecteur residu_int(nbddl); Mat_pleine raideur_int(nbddl,nbddl); // cas de la dynamique Mat_pleine matmasse(1,nbddl); // a priori on dimensionne en diagonale // cas des noeuds Tableau residus_extN(nombre->nbne); residus_extN(1) = new Vecteur(dim); Tableau raideurs_extN(nombre->nbne); raideurs_extN(1) = new Mat_pleine(dim,dim); for (int i = 2;i<= nombre->nbne;i++) { residus_extN(i) = residus_extN(1); raideurs_extN(i) = raideurs_extN(1); }; // cas des arêtes int nbddlA = nombre->nbneA * dim; Tableau residus_extA(nbA); residus_extA(1) = new Vecteur(nbddlA); Tableau raideurs_extA(nbA); raideurs_extA(1) = new Mat_pleine(nbddlA,nbddlA); for (int i = 2;i<= nbA;i++) { residus_extA(i) = residus_extA(1); raideurs_extA(i) = raideurs_extA(1); }; // cas des faces int nbddlS = nombre->nbneS * dim; // ddl pour les faces Tableau residus_extS(nbS); residus_extS(1) = new Vecteur(nbddlS); Tableau raideurs_extS(nbS); raideurs_extS(1) = new Mat_pleine(nbddlS,nbddlS); for (int i = 2;i<= nbS;i++) { residus_extS(i) = residus_extS(1); raideurs_extS(i) = raideurs_extS(1); }; // definition de la classe static contenant toute les variables communes aux HexaMemb HexaMemb::DonnComHexa* dodo; dodo = new DonnComHexa(hexa,tab_ddl,tab_ddlErr,tab_Err1Sig11,metri,resEr,raidEr,hexaEr,quadraS,seggS ,residu_int,raideur_int,residus_extN,raideurs_extN,residus_extA,raideurs_extA ,residus_extS,raideurs_extS,matmasse,hexaMas,nombre->nbi,hexaeHourg); return dodo; }; // destructions de certaines grandeurs pointées, créées au niveau de l'initialisation void HexaMemb::Destruction() { // tout d'abord l'idée est de détruire certaines grandeurs pointées que pour le dernier élément if ((unefois->nbelem_in_Prog == 0)&& (unefois->doCoMemb != NULL)) // cas de la destruction du dernier élément { HexaMemb::DonnComHexa* CoHexa = unefois->doCoMemb; // pour simplifier l'écriture int resErrTaille = CoHexa->resErr.Taille(); for (int i=1;i<= resErrTaille;i++) delete CoHexa->resErr(i); delete CoHexa->residus_externeN(1); delete CoHexa->raideurs_externeN(1); delete CoHexa->residus_externeA(1); delete CoHexa->raideurs_externeA(1); delete CoHexa->residus_externeS(1); delete CoHexa->raideurs_externeS(1); // suppression des éléments géométriques de référence delete CoHexa->hexaed; delete CoHexa->hexaedEr ; delete CoHexa->hexaedMas ; if (CoHexa->hexaedHourg != NULL) delete CoHexa->hexaedHourg; }; };