// FICHIER : BielletteC1.cc // CLASSE : BielletteC1 // This file is part of the Herezh++ application. // // The finite element software Herezh++ is dedicated to the field // of mechanics for large transformations of solid structures. // It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600) // INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) . // // Herezh++ is distributed under GPL 3 license ou ultérieure. // // Copyright (C) 1997-2022 Université Bretagne Sud (France) // AUTHOR : Gérard Rio // E-MAIL : gerardrio56@free.fr // // This program is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, // or (at your option) any later version. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty // of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. // See the GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program. If not, see . // // For more information, please consult: . //#include "Debug.h" # include using namespace std; //introduces namespace std #include #include "Sortie.h" #include "BielletteC1.h" #include "FrontPointF.h" #include "FrontSegLine.h" #include "TypeConsTens.h" //---------------------------------------------------------------- // def des donnees commune a tous les elements // la taille n'est pas defini ici car elle depend de la lecture //---------------------------------------------------------------- BielletteC1::DonneeCommune * BielletteC1::doCo = NULL; BielletteC1::UneFois BielletteC1::unefois; BielletteC1::NombresConstruire BielletteC1::nombre_V; BielletteC1::ConstrucElementbiel BielletteC1::construcElementbiel; // fonction définissant les nombres (de noeud, de point d'integ ..) // utilisé dans la construction des éléments //**** à modifier **** BielletteC1::NombresConstruire::NombresConstruire() { nbne = 2; // le nombre de noeud de l'élément nbneA = 2;// le nombre de noeud des aretes nbi = 1; // le nombre de point d'intégration pour le calcul mécanique nbiEr = 2;// le nombre de point d'intégration pour le calcul d'erreur nbiA = 1; // le nombre de point d'intégration pour le calcul de second membre linéique nbiMas = 2; // le nombre de point d'intégration pour le calcul de la matrice masse nbiHour = 0; // le nombre de point d'intégration un blocage d'hourglass }; // ---- definition du constructeur de la classe conteneur de donnees communes ------------ BielletteC1::DonneeCommune::DonneeCommune (GeomSeg& seg,DdlElement& tab,DdlElement& tabErr,DdlElement& tab_Err1Sig, Met_BielletteC1& met_bie, Tableau & resEr,Mat_pleine& raidEr, GeomSeg& seEr,Vecteur& residu_int,Mat_pleine& raideur_int, Tableau & residus_extN,Tableau & raideurs_extN, Tableau & residus_extA,Tableau & raideurs_extA, Mat_pleine& mat_masse,GeomSeg& seMa,int nbi,GeomSeg* segHourg) : segment(seg),tab_ddl(tab),tab_ddlErr(tabErr),tab_Err1Sig11(tab_Err1Sig) ,met_BielletteC1(met_bie) ,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) ,segmentEr(seEr) ,residu_interne(residu_int),raideur_interne(raideur_int) ,residus_externeN(residus_extN),raideurs_externeN(raideurs_extN) ,residus_externeA(residus_extA),raideurs_externeA(raideurs_extA) ,matrice_masse(mat_masse) ,segmentMas(seMa),segmentHourg(segHourg) { 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 (1); }; int tailledeps = d_epsBB.Taille(); for (int i=1;i<= tailledeps; i++) d_epsBB(i) = NevezTenseurBB (1); int tailledsig = d_sigHH.Taille(); for (int j=1;j<= tailledsig; j++) d_sigHH(j) = NevezTenseurHH (1); }; BielletteC1::DonneeCommune::DonneeCommune(DonneeCommune& a) : segment(a.segment),tab_ddl(a.tab_ddl),tab_ddlErr(a.tab_ddlErr),tab_Err1Sig11(a.tab_Err1Sig11) ,met_BielletteC1(a.met_BielletteC1),matGeom(a.matGeom),matInit(a.matInit) ,d2_epsBB(a.d2_epsBB),resErr(a.resErr),raidErr(a.raidErr),segmentEr(a.segmentEr) ,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) ,matrice_masse(a.matrice_masse),segmentMas(a.segmentMas),segmentHourg(a.segmentHourg) { 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 (*(a.d_epsBB(i))); int tailledsig = d_sigHH.Taille(); for (int j=1;j<= tailledsig; j++) d_sigHH(j) = NevezTenseurHH (*(d_sigHH(j))); }; BielletteC1::DonneeCommune::~DonneeCommune() { 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 -+-+-+ BielletteC1::UneFois::UneFois () : // constructeur par défaut doCoMemb(NULL),CalResPrem_t(0),CalResPrem_tdt(0),CalimpPrem(0),dualSortbiel(0) ,CalSMlin_t(0),CalSMlin_tdt(0),CalSMRlin(0) ,CalSMvol_t(0),CalSMvol_tdt(0),CalSMvol(0) ,CalDynamique(0),CalPt_0_t_tdt(0) ,nbelem_in_Prog(0) {}; BielletteC1::UneFois::~UneFois () { delete doCoMemb; }; // -+-+ fin definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+ // ---------- fin definition de la classe conteneur de donnees communes ------------ BielletteC1::BielletteC1 () : // Constructeur par defaut ElemMeca(),lesPtMecaInt(),donnee_specif() { lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca id_interpol=HERMITE; // donnees de la classe mere id_geom=POUT; // tab_noeud.Change_taille(nombre_V.nbne); // initialisation par défaut Init (); unefois.nbelem_in_Prog++; }; BielletteC1::BielletteC1 (double sect,int num_maill,int num_id): // Constructeur utile si la section de l'element et // le numero de l'element sont connus ElemMeca(num_maill,num_id,HERMITE,POUT),lesPtMecaInt(),donnee_specif(sect) { lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca tab_noeud.Change_taille(nombre_V.nbne); // initialisation Init (donnee_specif); unefois.nbelem_in_Prog++; }; // Constructeur fonction d'un numero de maillage et d'identification BielletteC1::BielletteC1 (int num_maill,int num_id) : ElemMeca(num_maill,num_id,HERMITE,POUT),lesPtMecaInt(),donnee_specif() { lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca tab_noeud.Change_taille(nombre_V.nbne); // initialisation par défaut Init (); unefois.nbelem_in_Prog++; }; BielletteC1::BielletteC1 (double sect,int num_maill,int num_id,const Tableau& tab): // Constructeur utile si la section de l'element, le numero de l'element et // le tableau des noeuds sont connus ElemMeca(num_maill,num_id,tab,HERMITE,POUT),donnee_specif(sect) { lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca if (tab_noeud.Taille() != 2) { cout << "\n erreur de dimensionnement du tableau de noeud \n"; cout << " BielletteC1::BielletteC1 (double sect,int num_mail,int num_id,const Tableau& tab)\n"; Sortie (1); } // construction du tableau de ddl spécifique à l'élément pour ses ConstTabDdl(); // initialisation bool sans_init_noeud = true; Init (donnee_specif,sans_init_noeud); unefois.nbelem_in_Prog++; }; BielletteC1::BielletteC1 (const BielletteC1& biel) : ElemMeca (biel),lesPtMecaInt(biel.lesPtMecaInt),donnee_specif(biel.donnee_specif) // Constructeur de copie { lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca // initialisation Init (donnee_specif); unefois.nbelem_in_Prog++; }; BielletteC1::~BielletteC1 () // Destructeur { LibereTenseur(); unefois.nbelem_in_Prog--; // defArete pointe sur la même grandeur que def donc pour éviter // une destruction lors du destructeur générale dans elemméca on le met à null defArete(1) = NULL; Destruction(); }; // Lecture des donnees de la classe sur fichier void BielletteC1::LectureDonneesParticulieres (UtilLecture * entreePrinc,Tableau * tabMaillageNoeud) { int nb; tab_noeud.Change_taille(2); for (int i=1; i<= 2; 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 != 2) 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 << "BielletteC1::LectureDonneesParticulieres"; Affiche(); Sortie (1); } } // construction du tableau de ddl des noeuds de BielletteC1 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 & BielletteC1::Point_physique(const Coordonnee& c_int,Coordonnee & co,Enum_dure temps) { BielletteC1::DonneeCommune* doCo = 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 == 0) { unefois.CalPt_0_t_tdt = 1; Tableau tab(3); tab(1)=iM0;tab(2)=iMt;tab(3)=iMtdt; doCo->met_BielletteC1.PlusInitVariables(tab) ; }; // b) calcul de l'interpolation const Vecteur& phi = doCo->segment.Phi_point(c_int); // c) calcul du point switch (temps) { case TEMPS_0 : co = doCo->met_BielletteC1.PointM_0(tab_noeud,phi); break; case TEMPS_t : co = doCo->met_BielletteC1.PointM_t(tab_noeud,phi); break; case TEMPS_tdt : co = doCo->met_BielletteC1.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 BielletteC1::Point_physique(const Coordonnee& c_int,Tableau & t_co) { BielletteC1::DonneeCommune* doCo = 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 == 0) { unefois.CalPt_0_t_tdt = 1; Tableau tab(3); tab(1)=iM0;tab(2)=iMt;tab(3)=iMtdt; doCo->met_BielletteC1.PlusInitVariables(tab) ; }; // b) calcul de l'interpolation const Vecteur& phi = doCo->segment.Phi_point(c_int); // c) calcul des point switch (t_co.Taille()) { case 3 : t_co(3) = doCo->met_BielletteC1.PointM_tdt(tab_noeud,phi); case 2 : t_co(2) = doCo->met_BielletteC1.PointM_t(tab_noeud,phi); case 1 : t_co(1) = doCo->met_BielletteC1.PointM_0(tab_noeud,phi); } }; // Calcul du residu local à t ou tdt en fonction du booleen atdt Vecteur* BielletteC1::CalculResidu (bool atdt,const ParaAlgoControle & pa) { BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture Tableau & d_epsBB = (doCo->d_epsBB);// " // dimensionnement de la metrique if (!atdt) {if( unefois.CalResPrem_t == 0) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; };} else {if( unefois.CalResPrem_tdt == 0) {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; doCo->met_BielletteC1.PlusInitVariables(tab) ; };}; // initialisation du résidu residu->Zero(); Vecteur poids =(doCo->segment).TaWi(); // poids d'interpolation = 2 poids(1) *= donnee_specif.section; ElemMeca::Cal_explicit ( doCo->tab_ddl,d_epsBB,nombre_V.nbi,poids,pa,atdt); return residu; }; // Calcul du residu local et de la raideur locale, // pour le schema implicite Element::ResRaid BielletteC1::Calcul_implicit (const ParaAlgoControle & pa) { BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture Tableau & d_epsBB = (doCo->d_epsBB);// " Tableau & d_sigHH = (doCo->d_sigHH);// " bool cald_Dvirtuelle = false; if (unefois.CalimpPrem == 0) { 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; doCo->met_BielletteC1.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(); Vecteur poids =(doCo->segment).TaWi(); // poids d'intégration = 2 poids(1) *= donnee_specif.section; ElemMeca::Cal_implicit (doCo->tab_ddl, d_epsBB,(doCo->d2_epsBB),d_sigHH, nombre_V.nbi,poids,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 * BielletteC1::CalculMatriceMasse (Enum_calcul_masse type_calcul_masse) { BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture // dimensionement de la métrique si nécessaire if (unefois.CalDynamique == 0) { 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; doCo->met_BielletteC1.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 = doCo->tab_ddl.NbDdl(); (doCo->matrice_masse).Initialise (nbddl,nbddl,0.); } }; Vecteur poids =(doCo->segmentMas).TaWi(); // poids d'intégration = 2 poids *= donnee_specif.section; // prise en compte de la section // appel de la routine générale ElemMeca::Cal_Mat_masse (doCo->tab_ddl,type_calcul_masse, nombre_V.nbiMas,(doCo->segmentMas).TaPhi(),nombre_V.nbne ,poids); return mat_masse; }; //------- calcul d'erreur, remontée des contraintes ------------------- // 1) calcul du résidu et de la matrice de raideur pour le calcul d'erreur Element::Er_ResRaid BielletteC1::ContrainteAuNoeud_ResRaid() { if(( unefois.CalResPrem_t == 0)|| (unefois.CalimpPrem == 0)) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; }; // appel du programme général ElemMeca:: SigmaAuNoeud_ResRaid(tab_noeud.Taille() ,(doCo->segment).TaPhi() ,(doCo->segment).TaWi() ,doCo-> resErr,doCo->raidErr ,(doCo->segmentEr).TaPhi() ,(doCo->segmentEr).TaWi()); return (Element::Er_ResRaid( &(doCo-> resErr),&(doCo->raidErr))); }; // 2) remontée aux erreurs aux noeuds Element::Er_ResRaid BielletteC1::ErreurAuNoeud_ResRaid() { if(( unefois.CalResPrem_t == 0)|| (unefois.CalimpPrem == 0)) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; }; // appel du programme général ElemMeca::Cal_ErrAuxNoeuds(tab_noeud.Taille(), (doCo->segment).TaPhi(), (doCo->segment).TaWi(),doCo-> resErr ); return (Element::Er_ResRaid( &(doCo-> resErr),&(doCo->raidErr))); }; // actualisation des ddl et des grandeurs actives de t+dt vers t void BielletteC1::TdtversT() { lesPtMecaInt.TdtversT(); // contrainte if (tabSaveDon(1) != NULL) tabSaveDon(1)->TdtversT(); if (tabSaveTP(1) != NULL) tabSaveTP(1)->TdtversT(); if (tabSaveDefDon(1) != NULL) tabSaveDefDon(1)->TdtversT(); ElemMeca::TdtversT_(); // appel de la procédure mère }; // actualisation des ddl et des grandeurs actives de t vers tdt void BielletteC1::TversTdt() { lesPtMecaInt.TversTdt(); // contrainte if (tabSaveDon(1) != NULL) tabSaveDon(1)->TversTdt(); if (tabSaveTP(1) != NULL) tabSaveTP(1)->TversTdt(); if (tabSaveDefDon(1) != NULL)tabSaveDefDon(1)->TversTdt(); ElemMeca::TversTdt_(); // appel de la procédure mère }; // 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 // type indique le type de calcul d'erreur : void BielletteC1::ErreurElement(int type,double& errElemRelative ,double& numerateur, double& denominateur) { if(( unefois.CalResPrem_t == 0)|| (unefois.CalimpPrem == 0)) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; }; // appel du programme général ElemMeca::Cal_ErrElem(type,errElemRelative,numerateur,denominateur, tab_noeud.Taille(),(doCo->segment).TaPhi(), (doCo->segment).TaWi(), (doCo->segmentEr).TaPhi(),(doCo->segmentEr).TaWi()); }; //============= 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 BielletteC1::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 la BielletteC1 switch (cas) { case 1 : // ------- on récupère tout ------------------------- { // construction du tableau de ddl des noeuds du triangle 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 << "BielletteC1::Lecture_base_info(ofstream& sort,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 BielletteC1::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 la BielletteC1 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 << "BielletteC1::Ecriture_base_info(ofstream& sort,int cas)" << " cas= " << cas << endl; Sortie(1); } } }; // Calcul de la matrice géométrique et initiale ElemMeca::MatGeomInit BielletteC1::MatricesGeometrique_Et_Initiale (const ParaAlgoControle & pa) { BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture Tableau & d_epsBB = (doCo->d_epsBB);// " Tableau & d_sigHH = (doCo->d_sigHH);// " bool cald_Dvirtuelle = false; if (unefois.CalimpPrem == 0) { 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; doCo->met_BielletteC1.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; }; // Par simplicité Mat_pleine & matGeom = doCo->matGeom; Mat_pleine & matInit = doCo->matInit; // mise à zéro de la matrice géométrique matGeom.Initialise(); Vecteur poids =(doCo->segment).TaWi(); // poids d'interpolation = 2 poids(1) *= donnee_specif.section; ElemMeca::Cal_matGeom_Init (matGeom,matInit,doCo->tab_ddl,d_epsBB, doCo->d2_epsBB,d_sigHH,1,poids,pa,cald_Dvirtuelle); return MatGeomInit(&matGeom,&matInit); } ; // retourne les tableaux de ddl associés aux noeuds, gere par l'element // ce tableau et specifique a l'element const DdlElement & BielletteC1::TableauDdl() const { return doCo->tab_ddl; }; // liberation de la place pointee void BielletteC1::Libere () {Element::Libere (); // liberation de residu et raideur LibereTenseur() ; // liberation des tenseur intermediaires }; // acquisition ou modification d'une loi de comportement void BielletteC1::DefLoi (LoiAbstraiteGeneral * NouvelleLoi) { // verification du type de loi if ((NouvelleLoi->Dimension_loi() != 1) && (NouvelleLoi->Dimension_loi() != 4)) { cout << "\n Erreur, la loi de comportement a utiliser avec des BielletteC1s"; cout << " doit etre de type 1D, \n ici est de type = " << (NouvelleLoi->Dimension_loi()) << " !!! " << endl; Sortie(1); } // cas d'une loi mécanique if (GroupeMecanique(NouvelleLoi->Id_categorie())) {loiComp = (Loi_comp_abstraite *) NouvelleLoi; // initialisation du stockage particulier, ici 1 pt d'integ tabSaveDon(1) = loiComp->New_et_Initialise(); // idem pour le type de déformation mécanique associé tabSaveDefDon(1) = 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, ici 1 pt d'integ tabSaveTP(1) = 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 BielletteC1::TestComplet() { int res = ElemMeca::TestComplet(); // test dans la fonction mere if ( donnee_specif.section == section_defaut) { cout << "\n la section de la BielletteC1 n'est pas defini \n"; res = 0; } if ( tab_noeud(1) == NULL) { cout << "\n les noeuds de la BielletteC1 ne sont pas defini \n"; res = 0; } else { int testi =1; int posi = Id_nom_ddl("X1") -1; int dim = ParaGlob::Dimension(); for (int i =1; i<= dim; i++) for (int j=1;j<=2;j++) if(!(tab_noeud(j)->Existe_ici(Enum_ddl(posi+i)))) testi = 0; if(testi == 0) { cout << "\n les ddls X1,X2 etc des noeuds de la BielletteC1 ne sont pas defini \n"; cout << " \n utilisez BielletteC1::ConstTabDdl() pour completer " ; res = 0; } } return res; }; // ajout du tableau de ddl des noeuds de BielletteC1 void BielletteC1::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 tab_noeud(1)->PlusTabDdl(ta); tab_noeud(2)->PlusTabDdl(ta); }; // procesure permettant de completer l'element apres // sa creation avec les donnees du bloc transmis // peut etre appeler plusieurs fois Element* BielletteC1::Complete(BlocGen & bloc,LesFonctions_nD* lesFonctionsnD) { // complétion avec bloc if (bloc.Nom(1) == "sections") { donnee_specif.section = bloc.Val(1); return this; } else return ElemMeca::Complete_ElemMeca(bloc,lesFonctionsnD); }; // Compléter pour la mise en place de la gestion de l'hourglass Element* BielletteC1::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 BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; ElemMeca::Init_hourglass_comp(*(doCo->segmentHourg),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)) &&(doCo->segmentHourg == NULL)) type_stabHourglass = STABHOURGLASS_NON_DEFINIE; return this; }; // affichage dans la sortie transmise, des variables duales "nom" // dans le cas ou nom est vide, affichage de "toute" les variables void BielletteC1::AfficheVarDual(ofstream& sort, Tableau& nom) {// affichage de l'entête de l'element sort << "\n******************************************************************"; sort << "\n Element bielette (2 noeuds 1 point d'integration) "; sort << "\n******************************************************************"; // appel de la procedure de elem meca if (!(unefois.dualSortbiel) && (unefois.CalimpPrem)) { VarDualSort(sort,nom,1,1); unefois.dualSortbiel += 1; } else if ((unefois.dualSortbiel) && (unefois.CalimpPrem)) VarDualSort(sort,nom,1,11); else if (!(unefois.dualSortbiel) && (unefois.CalResPrem_tdt)) { VarDualSort(sort,nom,1,2); unefois.dualSortbiel += 1; } else if ((unefois.dualSortbiel ) && (unefois.CalResPrem_tdt)) VarDualSort(sort,nom,1,12); // sinon on ne fait rien }; // récupération des valeurs au numéro d'ordre = iteg pour // les grandeur enu // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière Tableau BielletteC1::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.dualSortbiel) && (unefois.CalimpPrem)) { cas=1;unefois.dualSortbiel += 1; } else if ((unefois.dualSortbiel) && (unefois.CalimpPrem)) { cas = 11;} else if (!(unefois.dualSortbiel) && (unefois.CalResPrem_tdt)) { cas=2;unefois.dualSortbiel += 1; } else if ((unefois.dualSortbiel) && (unefois.CalResPrem_tdt)) { cas = 12;} // sinon pour l'instant pb, car il faut définir des variable dans la métrique 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.dualSortbiel= " << unefois.dualSortbiel << " unefois.CalimpPrem= " << unefois.CalimpPrem << "\n BielletteC1::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 BielletteC1::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.dualSortbiel == 0) && (unefois.CalimpPrem != 0)) { cas=1;unefois.dualSortbiel = 1; } else if ((unefois.dualSortbiel == 1) && (unefois.CalimpPrem != 0)) { cas = 11;} else if ((unefois.dualSortbiel == 0) && (unefois.CalResPrem_tdt != 0)) { cas=2;unefois.dualSortbiel = 1; } else if ((unefois.dualSortbiel == 1) && (unefois.CalResPrem_tdt != 0)) { cas = 12;} // sinon pour l'instant pb, car il faut définir des variable dans la métrique 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.dualSortbiel= " << unefois.dualSortbiel << " unefois.CalimpPrem= " << unefois.CalimpPrem << "\n BielletteC1::ValTensorielle_a_diff_temps(Enum_dure enu_t..."; Sortie(1); }; ElemMeca::Valeurs_Tensorielles(absolue,enu_t,enu,iteg,cas); }; // cas d'un chargement volumique, // force indique la force volumique appliquée // retourne le second membre résultant // ici on considère la section de la BielletteC1 pour constituer le volume // -> explicite à t ou tdt en fonction de la variable booléenne atdt Vecteur BielletteC1::SM_charge_volumique_E (const Coordonnee& force,Fonction_nD* pt_fonct,bool atdt,const ParaAlgoControle & pa,bool sur_volume_finale_) { BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if (!atdt) {if(( unefois.CalResPrem_t == 0) && (unefois.CalimpPrem == 0) && (unefois.CalSMvol_t == 0)) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; };} else {if(( unefois.CalResPrem_tdt == 0) && (unefois.CalimpPrem == 0) && (unefois.CalSMvol_tdt == 0)) { unefois.CalSMvol_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; doCo->met_BielletteC1.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 (donnee_specif.section * ElemMeca::SM_charge_vol_E (doCo->tab_ddl,(doCo->segment).TaPhi() ,tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa,sur_volume_finale_)); }; // 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 BielletteC1::SMR_charge_volumique_I (const Coordonnee& force,Fonction_nD* pt_fonct,const ParaAlgoControle & pa,bool sur_volume_finale_) { BielletteC1::DonneeCommune* doCo = 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.CalimpPrem == 0) && (unefois.CalSMvol_tdt == 0)) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; }; // appel du programme général d'elemmeca ElemMeca::SMR_charge_vol_I (doCo->tab_ddl,(doCo->segment).TaPhi() ,tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa,sur_volume_finale_); // prise en compte de la section (*residu) *= donnee_specif.section; (*raideur) *= donnee_specif.section; Element::ResRaid el; el.res = residu; el.raid = raideur; return el; }; // cas d'un chargement lineique, sur l'aretes de la BielletteC1 // force indique la force lineique appliquée // numarete indique le numéro de l'arete chargée ici 1 // retourne le second membre résultant // -> explicite à t ou tdt en fonction de la variable booléenne atdt Vecteur BielletteC1::SM_charge_lineique_E(const Coordonnee& force,Fonction_nD* pt_fonct,int ,bool atdt,const ParaAlgoControle & pa) { // initialisation du vecteur résidu ((*res_extA)(1))->Zero(); // on récupère ou on crée la frontière arrête Frontiere_lineique(1,true); BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if (!atdt) {if( unefois.CalSMlin_t == 0) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; };} else {if( unefois.CalSMlin_tdt == 0) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; };}; // on définit la déformation a doc si elle n'existe pas déjà if (defArete(1) == NULL) defArete(1) = def; // a priori idem que la BielletteC1 // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SM_charge_line_E (doCo->tab_ddl,1,(doCo->segment).TaPhi() ,tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa); }; // 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 ici 1 par défaut // 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 BielletteC1::SMR_charge_lineique_I (const Coordonnee& force,Fonction_nD* pt_fonct,int ,const ParaAlgoControle & pa) { // initialisation du vecteur résidu et de la raideur ((*res_extA)(1))->Zero(); ((*raid_extA)(1))->Zero(); // on récupère ou on crée la frontière arrête Frontiere_lineique(1,true); BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if( unefois.CalSMRlin == 0) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; }; // on définit la déformation a doc si elle n'existe pas déjà if (defArete(1) == NULL) defArete(1) = def; // a priori idem que la BielletteC1 // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SMR_charge_line_I (doCo->tab_ddl,1 ,(doCo->segment).TaPhi(),tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa); }; // cas d'un chargement lineique suiveuse, sur l'aretes frontière de la BielletteC1 (2D uniquement) // force indique la force lineique appliquée // numarete indique le numéro de l'arete chargée // retourne le second membre résultant // -> explicite à t ou tdt en fonction de la variable booléenne atdt Vecteur BielletteC1::SM_charge_lineique_Suiv_E(const Coordonnee& force,Fonction_nD* pt_fonct,int ,bool atdt,const ParaAlgoControle & pa) { // initialisation du vecteur résidu ((*res_extA)(1))->Zero(); // on récupère ou on crée la frontière arrête Frontiere_lineique(1,true); BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if (!atdt) {if( unefois.CalSMlin_t == 0) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; };} else {if( unefois.CalSMlin_tdt == 0) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; };}; // on définit la déformation a doc si elle n'existe pas déjà if (defArete(1) == NULL) defArete(1) = def; // a priori idem que la BielletteC1 // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SM_charge_line_Suiv_E (doCo->tab_ddl,1,(doCo->segment).TaPhi() ,tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa); }; // cas d'un chargement lineique suiveuse, sur l'aretes frontière de la BielletteC1 (2D uniquement) // force indique la force lineique appliquée // numarete indique le numéro de l'arete 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 BielletteC1::SMR_charge_lineique_Suiv_I (const Coordonnee& force,Fonction_nD* pt_fonct,int ,const ParaAlgoControle & pa) { // initialisation du vecteur résidu et de la raideur ((*res_extA)(1))->Zero(); ((*raid_extA)(1))->Zero(); // on récupère ou on crée la frontière arrête Frontiere_lineique(1,true); BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture // dimensionnement de la metrique if( unefois.CalSMRlin == 0) { 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; doCo->met_BielletteC1.PlusInitVariables(tab) ; }; // on définit la déformation a doc si elle n'existe pas déjà if (defArete(1) == NULL) defArete(1) = def; // a priori idem que la BielletteC1 // appel du programme général d'elemmeca et retour du vecteur second membre return ElemMeca::SMR_charge_line_Suiv_I (doCo->tab_ddl,1 ,(doCo->segment).TaPhi(),tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,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 BielletteC1::SM_charge_hydrodynamique_E( Courbe1D* frot_fluid,const double& poidvol , Courbe1D* coef_aero_n,int num,const double& coef_mul , Courbe1D* coef_aero_t,bool atdt,const ParaAlgoControle & pa) { // initialisation du vecteur résidu ((*res_extN)(num))->Zero(); // ici les frontières sont des points BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture // on récupère ou on crée la frontière points Frontiere_points(num,true); // on utilise l'élément géométique de l'élément, pour TaPhi et TaWi pour les passages de paramètres, // mais ces infos ne sont pas utilisées car pour les frontières points, pas de points d'integ // appel du programme général d'ElemMeca et retour du vecteur second membre return ElemMeca::SM_charge_hydrodyn_E (poidvol,doCo->segment.TaPhi(),1 ,frot_fluid,doCo->segment.TaWi() ,coef_aero_n,num,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 BielletteC1::SMR_charge_hydrodynamique_I( Courbe1D* frot_fluid,const double& poidvol , Courbe1D* coef_aero_n,int num,const double& coef_mul , Courbe1D* coef_aero_t,const ParaAlgoControle & pa) { // initialisation du vecteur résidu et de la raideur ((*res_extN)(num))->Zero(); ((*raid_extN)(num))->Zero(); BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture // on récupère ou on crée la frontière points Frontiere_points(num,true); // on utilise l'élément géométique de l'élément, pour TaPhi et TaWi pour les passages de paramètres, // mais ces infos ne sont pas utilisées car pour les frontières points, pas de points d'integ // appel du programme général d'ElemMeca et retour du vecteur second membre return ElemMeca::SM_charge_hydrodyn_I (poidvol,doCo->segment.TaPhi(),1 ,frot_fluid,doCo->segment.TaWi(),tabb(posi_tab_front_point+num)->DdlElem() ,coef_aero_n,num,coef_mul,coef_aero_t,pa); }; // Calcul des frontieres 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 & BielletteC1::Frontiere(bool force) { int cas = 6; // on veut des lignes et des points 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_point == 0) || force || (ind_front_point == 2)) // { // dimensionnement des tableaux intermediaires // Tableau tab(1); // les noeuds des points frontieres // DdlElement ddelem(1); // les ddlelements des points frontieres // int tail; // if ((ParaGlob::Dimension() == 1) && (ind_front_lin > 0)) // tail = 3; // deux points et une ligne // else if (ParaGlob::Dimension() == 1) // cas sans ligne // tail = 2; // 2 points // else // cas d'une dimension 2 et 3 // { tail = 3; // deux points et une ligne // ind_front_lin = 1; // } // tabb.Change_taille(tail); // le tableau total de frontières // // // premier point // tab(1) = tab_noeud(1); // ddelem.Change_un_ddlNoeudElement(1,doCo->tab_ddl(1)); // if (tabb(1+posi_tab_front_point) == NULL) // tabb(1+posi_tab_front_point) = new FrontPointF (tab,ddelem); // // second point // tab(1) = tab_noeud(2); // ddelem.Change_un_ddlNoeudElement(1,doCo->tab_ddl(2)); // if (tabb(2+posi_tab_front_point) == NULL) // tabb(2+posi_tab_front_point) = new FrontPointF (tab,ddelem); // // 3 ieme cote eventuelle // if (ind_front_lin > 0) // // cas où il y a une ligne, c'est forcément le premier élément // if (tabb(1) == NULL) // tabb(1) = new FrontSegLine(tab_noeud,doCo->tab_ddl); // // // mise à jour des indicateurs // ind_front_point = 1; // } // // return 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 BielletteC1::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 > 2)||(num <=0)) // { cout << "\n *** erreur, pour les BielletteC1s il n'y a que deux frontieres point ! " // << "\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 // if ((ind_front_lin > 0) && (ind_front_point==0)) // // cas où il y a une frontière ligne (mais pas de point) // { tabb.Change_taille(3); posi_tab_front_point = 1;} // else if ((ind_front_lin == 0) && (ind_front_point==0)) // // cas où aucune frontière existe, on crée pour les points // { tabb.Change_taille(2); posi_tab_front_point = 0;}; // // dans les deux autres cas ((ind_front_lin == 0) && (ind_front_point>0)) et // // ((ind_front_lin > 0) && (ind_front_point>0)), les points existant déjà on n'a rien n'a faire // // on définit les deux points par simplicité // // premier point // tab(1) = tab_noeud(1); // ddelem.Change_un_ddlNoeudElement(1,doCo->tab_ddl(1)); // if (tabb(1+posi_tab_front_point) == NULL) // tabb(1+posi_tab_front_point) = new FrontPointF (tab,ddelem); // // second point // tab(1) = tab_noeud(2); // ddelem.Change_un_ddlNoeudElement(1,doCo->tab_ddl(2)); // if (tabb(2+posi_tab_front_point) == NULL) // tabb(2+posi_tab_front_point) = new FrontPointF (tab,ddelem); // ind_front_point=1; // mise à jour de l'indicateur // }; // 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 BielletteC1::Frontiere_lineique(int num,bool ) // { // 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 != 1) // { cout << "\n *** erreur, pour les BielletteC1s il n'y a qu'une frontière ligne ! " // << "\n Frontiere_lineique(int num,bool force)"; // Sortie(1); // } // #endif // // // on regarde si les frontières linéiques existent sinon on les crée // if (ind_front_lin == 1) // return (ElFrontiere*)tabb(1); // else if ( ind_front_lin == 2) // // cas où certaines frontières existent // if (tabb(1) != NULL) // return (ElFrontiere*)tabb(1); // // dans tous les autres cas on construit la frontière ligne // // on commence par dimensionner le tableau de frontière // if (ind_front_point > 0) // // cas où il y a des frontières points (mais pas ligne) // // on décale le tableau // { tabb.Change_taille(3); // tabb(3) = tabb(2); // tabb(2) = tabb(1); // posi_tab_front_point = 1; // } // else // // cas d'une frontière linéique // tabb.Change_taille(1); // // on définit la ligne // tabb(1) = new FrontSegLine(tab_noeud,doCo->tab_ddl); // ind_front_lin = 1; // mise à jour de l'indicateur // // et normalement posi_tab_front_ligne = 0, car jamais changé // return (ElFrontiere*)tabb(num); // }; // // //// ramène la frontière surfacique //// éventuellement création des 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) //// ici normalement la fonction ne doit pas être appelée //ElFrontiere* const BielletteC1::Frontiere_surfacique(int ,bool ) // { cout << "\n *** erreur, pour les BielletteC1s il n'y a pas de frontiere surface ! " // << "\n Frontiere_surfacique(int ,bool force = false)"; // Sortie(1); // return NULL; // }; // =====>>>> methodes privées appelees par les classes dérivees <<<<===== // fonction d'initialisation servant au niveau du constructeur void BielletteC1::Init (Donnee_specif donnee_spec,bool sans_init_noeud) { // bien que la grandeur donnee_specif est défini dans la classe generique // le fait de le passer en paramètre permet de tout initialiser dans Init // et ceci soit avec les valeurs par défaut soit avec les bonnes valeurs donnee_specif =donnee_spec; // le fait de mettre les pointeurs a null permet // de savoir que l'element n'est pas complet tab_noeud.Change_taille(nombre_V.nbne); // 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_V.nbne;i++) tab_noeud(i) = NULL; // definition des donnees communes aux BielletteC1xxx // a la premiere definition d'une instance if (unefois.doCoMemb == NULL) BielletteC1::Def_DonneeCommune(); unefois.doCoMemb = doCo ; met = &(doCo->met_BielletteC1); // met est defini dans elemeca // def pointe sur la deformation specifique a l'element pour le calcul mecanique def = new Deformation(*met,tab_noeud,(doCo->segment).TaDphi(),(doCo->segment).TaPhi()); // idem pour la remontee aux contraintes et le calcul d'erreur defEr = new Deformation(*met,tab_noeud,(doCo->segmentEr).TaDphi(),(doCo->segmentEr).TaPhi()); // idem pour la remontee aux contraintes et le calcul d'erreur defMas = new Deformation(*met,tab_noeud,(doCo->segmentMas).TaDphi(),(doCo->segmentMas).TaPhi()); // idem pour le calcul de second membre defArete.Change_taille(1); // 1 arrête utilisée pour le second membre // la déformation sera construite si nécessaire au moment du calcul de second membre defArete(1) = NULL; //dimensionnement des deformations et contraintes etc.. int dimtens = 1; lesPtMecaInt.Change_taille_PtIntegMeca(nombre_V.nbi,dimtens); // attribution des numéros de référencement dans le conteneur for (int ni = 1; ni<= nombre_V.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 au point d'integ tabSaveDon.Change_taille(nombre_V.nbi); tabSaveTP.Change_taille(nombre_V.nbi); tabSaveDefDon.Change_taille(nombre_V.nbi); tab_energ.Change_taille(nombre_V.nbi); tab_energ_t.Change_taille(nombre_V.nbi); // initialisation des pointeurs définis dans la classe Element concernant les résidus et // raideur // --- cas de la puissance interne --- residu = &(doCo->residu_interne); // residu local raideur = &(doCo->raideur_interne); // raideur locale // --- cas de la dynamique ----- mat_masse = &(doCo->matrice_masse); // --- cas des efforts externes concernant les noeuds ------ res_extN = &(doCo->residus_externeN); // pour les résidus et second membres raid_extN= &(doCo->raideurs_externeN);// pour les raideurs // --- cas des efforts externes concernant les aretes ------ res_extA = &(doCo->residus_externeA); // pour les résidus et second membres raid_extA= &(doCo->raideurs_externeA);// pour les raideurs }; // fonction privee // dans cette fonction il ne doit y avoir que les données communes !! void BielletteC1::Def_DonneeCommune() { int nbn = nombre_V.nbne; // interpollation : element geometrique correspondant: 1 pt integ, 2 noeuds GeomSeg seg(nombre_V.nbi,nbn) ; // degre de liberte int dim = ParaGlob::Dimension(); DdlElement tab_ddl(nbn,dim); int posi = Id_nom_ddl("X1") -1; for (int i =1; i<= dim; i++) for (int j=1; j<= nbn; 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 int nbcomposante = ParaGlob::NbCompTens (); DdlElement tab_ddlErr(nbn,nbcomposante); posi = Id_nom_ddl("SIG11") -1; for (int j=1; j<= nbn; j++) {// on definit le nombre de composante de sigma en absolu if (nbcomposante == 3) // dans l'énumération les composantes ne sont pas a suivre {//tab_ddlErr (j,1) = Enum_ddl(1+posi); // cas de SIG11 //tab_ddlErr (j,2) = Enum_ddl(2+posi); // cas de SIG22 //tab_ddlErr (j,3) = Enum_ddl(4+posi); // cas de SIG12 tab_ddlErr.Change_Enum(j,1,Enum_ddl(1+posi)); // cas de SIG11 tab_ddlErr.Change_Enum(j,2,Enum_ddl(2+posi)); // cas de SIG22 tab_ddlErr.Change_Enum(j,3,Enum_ddl(4+posi)); // cas de SIG12 } else if (nbcomposante == 6) // les composantes sont a suivre dans l'enumération for (int i= 1;i<= nbcomposante; i++) // tab_ddlErr (j,i) = Enum_ddl(i+posi); tab_ddlErr.Change_Enum(j,i,Enum_ddl(i+posi)); else // tab_ddlErr (j,1) = Enum_ddl(1+posi); // uniquement SIG11 tab_ddlErr.Change_Enum(j,1,Enum_ddl(1+posi)); // uniquement SIG11 } // egalement pour tab_Err1Sig11, def d'un tableau de un ddl : enum SIG11 // par noeud DdlElement tab_Err1Sig11(nbn,DdlNoeudElement(SIG11)); // toujours pour le calcul d'erreur definition des fonctions d'interpolation // pour le calcul du hession de la fonctionnelle : // 2 points d'integration et 2 noeuds GeomSeg segEr(nombre_V.nbiEr,nbn) ; // pour le calcul de la matrice masse definition des fonctions d'interpolation // nbiMas points d'integration et nbn noeuds, en particulier pour le calcul de la masse consistante GeomSeg segMa(nombre_V.nbiMas,nbn) ; // pour le calcul relatifs à la stabilisation d'hourglass GeomSeg* segmentHourg = NULL; if (nombre_V.nbiHour > 0) segmentHourg = new GeomSeg(nombre_V.nbiHour,nombre_V.nbne); // def metrique // on definit les variables a priori toujours utiles Tableau tab(24); 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 ; tab(16)= id_gijBB_tdt; tab(17)= id_giH_tdt; tab(18)= id_gijHH_tdt; tab(19)= idMtdt ; tab(20)= id_jacobien_tdt;tab(21)= id2_gijBB_tdt; tab(22)= igradVBB_tdt; tab(23) = iVtdt; tab(24)= idVtdt; // dim du pb , nb de vecteur de la base , tableau de ddl et la def de variables Met_BielletteC1 metri(ParaGlob::Dimension(),tab_ddl,tab,nbn) ; // ---- cas du calcul d'erreur sur sigma ou epsilon // les tenseurs sont exprimees en absolu donc nombre de composante fonction // de la dimension absolue Tableau resEr(nbcomposante); for (int i = 1;i<= nbcomposante; i++) resEr(i)=new Vecteur (nbn); // une composante par noeud Mat_pleine raidEr(nbn,nbn); // la raideur pour l'erreur // dimensionnement des différents résidus et raideurs pour le calcul mécanique int nbddl = tab_ddl.NbDdl(); 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 // il y a deux extrémités mais identiques Tableau residus_extN(2); residus_extN(1) = new Vecteur(dim); residus_extN(2) = residus_extN(1); int nbddlA = nombre_V.nbneA * dim; int nbA = 1; // 1 arêtes Tableau residus_extA(nbA); residus_extA(1) = new Vecteur(nbddlA); Tableau raideurs_extA(nbA); raideurs_extA(1) = new Mat_pleine(nbddlA,nbddlA); Tableau raideurs_extN(2);raideurs_extN(1) = new Mat_pleine(dim,dim); raideurs_extN(2) = raideurs_extN(1); // definition de la classe static contenant toute les variables communes aux BielletteC1s doCo = new DonneeCommune(seg,tab_ddl,tab_ddlErr,tab_Err1Sig11,metri,resEr,raidEr,segEr, residu_int,raideur_int,residus_extN,raideurs_extN,residus_extA,raideurs_extA ,matmasse,segMa,nombre_V.nbi,segmentHourg); }; // destructions de certaines grandeurs pointées, créées au niveau de l'initialisation void BielletteC1::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 { BielletteC1::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture int resErrTaille = doCo->resErr.Taille(); for (int i=1;i<= resErrTaille;i++) delete doCo->resErr(i); delete doCo->residus_externeN(1); delete doCo->raideurs_externeN(1); delete doCo->residus_externeA(1); delete doCo->raideurs_externeA(1); if (doCo->segmentHourg != NULL) delete doCo->segmentHourg; } };