// 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 "Charge.h" #include "ConstMath.h" #include "MathUtil.h" #include "ReferenceNE.h" #include "ReferenceAF.h" // mise en place du chargement impose sur le second membre global // cas explicite // la variable relative à la classe assemblage est passé en argument ce qui permet d'utiliser // éventuellement plusieurs type d'assemblage // retour false si problème de calcul bool Charge::ChargeSecondMembre_Ex_mecaSolid (Assemblage & assemb,LesMaillages * lesMail,LesReferences* lesRef,Vecteur& vecglob ,const ParaAlgoControle & pa,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD) { temps_cpu_chargement.Mise_en_route_du_comptage(); // temps cpu bool retour = true; // init par défaut try { #ifndef UTILISATION_MPI // pour tous les maillages, pour tous les éléments on effectue une initialisation éventuelle int nb_mail = lesMail->NbMaillage(); for (int imail = 1;imail<=nb_mail;imail++) { int nb_ele = lesMail->Nombre_element(imail); for (int iele =1;iele<=nb_ele;iele++) lesMail->Element_LesMaille(imail,iele).Initialisation_avant_chargement(); }; #else // ParaGlob::Monde()->barrier(); // synchronisation ici de tous les process // cas d'un calcul multi-CPU , si cpu 0, idem cas non parallèle if (ParaGlob::Monde()->rank() == 0) {int nb_mail = lesMail->NbMaillage(); for (int imail = 1;imail<=nb_mail;imail++) { int nb_ele = lesMail->Nombre_element(imail); for (int iele =1;iele<=nb_ele;iele++) lesMail->Element_LesMaille(imail,iele).Initialisation_avant_chargement(); }; } else {// on va parcourir les éléments associés au cpu // on récupère la distribution d'éléments concernant le cpu en cours // si c'est le cpu 0 on initialise tous les éléments const Tableau < list >& tab_list_elem_cpu = ParaGlob::param->const_List_element_CPU_en_cours(); // la taille est identique à nbMailMax, sauf si c'est le cpu 0, là le tableau est vide et il n'y // aura pas de boucle, (ou plutôt on sortira directement de la boucle // le cas CPU 0 est traité à part int nb_mail_distrib = tab_list_elem_cpu.Taille(); for (int imail =1; imail<= nb_mail_distrib; imail++) {const list & list_elem_cpu= tab_list_elem_cpu(imail); // pour simplifier // on balaie les éléments nécessaires list ::const_iterator il,ilfin=list_elem_cpu.end(); for (il = list_elem_cpu.begin();il != ilfin;il++) { int iele = (*il); lesMail->Element_LesMaille(imail,iele).Initialisation_avant_chargement(); }; }; }; #endif int posi = Id_nom_ddl("X1") -1; int dim = ParaGlob::Dimension(); // on regarde si l'espace de travail est axi-symétrique, si oui on utilise une dimension réduite int dima=dim; // c'est dima qui est utiliser pour la mise en place des efforts ponctuelles if (ParaGlob::AxiSymetrie()) dima--; // $$$ --- cas des forces ponctuelles --- $$$ #ifdef UTILISATION_MPI // cas d'un calcul //, les chargements aux noeuds sont calculés par le cpu 0 uniquement if (ParaGlob::Monde()->rank() == 0) #endif // on parcours le tableau tabPonctuel for (int i=1;i<= tabPonctuel.Taille();i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabPonctuel_i = tabPonctuel(i); const ReferenceNE & ref = ((ReferenceNE &) lesRef->Trouve(tabPonctuel_i.NomRef(),tabPonctuel_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabPonctuel_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabPonctuel_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabPonctuel_i.Pas_a_prendre_en_compte(Temps_courant()))) {int reftaille = ref.Taille(); int cas_assemblage = (assemb.Nb_cas_assemb()).n; // récup du cas d'assemblage for (int nn =1; nn<= reftaille;nn++) // nn = position des num de noeud { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn)); Coordonnee force; // init // on traite en fonction du fait que c'est un champ ou non switch (tabPonctuel_i.Champ()) { case 1: // c'est un champ de valeurs force = (tabPonctuel_i.Vect_de_coordonnee(nn)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_nn_3 = lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,3)); force(3) = pt_courbe_force_nn_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_nn_2 = lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,2)); force(2) = pt_courbe_force_nn_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_nn_1 = lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,1)); force(1) = pt_courbe_force_nn_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ force = (tabPonctuel_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales // ou au point + grandeurs globales éventuelles const string & nomF_charge = tabPonctuel_i.NomF_charge(); if(nomF_charge.length()) {Fonction_nD * pt_fonct = lesFonctionsnD->Trouve(nomF_charge); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales { // on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {force *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : force(3) *= tava(3); case 2 : force(2) *= tava(2); case 1 : force(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; } else if (pt_fonct->Depend_M() != 0) // signifie qu'il y a une dépendance au point M { if ((pt_fonct->Depend_M()<0) && pt_fonct->Nom_variables().Taille() == dima) { // on vérifie également qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace const Coordonnee& M_tdt = noe.Coord2(); // comme pt_fonct->Depend_M()<0, cela signifie que la fonction dépend que de M // donc dans les paramètres d'appel on ne transmet que M Tableau tab_M(1,M_tdt); Tableau & tava = pt_fonct->Val_FnD_Evoluee(NULL,&tab_M,NULL); // pour simplifier if (pt_fonct->NbComposante() == 1) {force *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : force(3) *= tava(3); case 2 : force(2) *= tava(2); case 1 : force(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef() << " en retour, le nombre de composante est different de 1 ou de la dimension de l'espace : " << dima << endl; Sortie(1); }; } else { if ((pt_fonct->Depend_M()>0)|| pt_fonct->Nom_variables().Taille() != dima) {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD utilisant les coordonnees du point" << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef() << " actuellement la fonction doit avoir uniquement comme parametres, les coordonnees du point et eventuellement des variables globales "; cout << "\n fonction nD: "; pt_fonct->Affiche(); Sortie(1); }; }; } }; if ((noe.Pointeur_assemblage(X1,cas_assemblage) != -1) && (noe.En_service(X1)) && (noe.UneVariable(X1))) // les forces ponctuelles correspondent aux ddl X1 X2 ... for (int di=1;di<= dima;di++) {// maintenant deux cas: soit avec ou sans courbes de charge int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage); if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD // cas d'une grandeur fixe { vecglob (pt_adres) += coeff * (force(di));} else // cas d'une fonction nD sans courbe { vecglob (pt_adres) += (force(di))* tabPonctuel_i.Echelle_courbe();} } else { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabPonctuel_i.Echelle_courbe(); vecglob (pt_adres) += coeff_charge * (force(di)); }; }; }; }; }; // $$$ --- cas des torseurs de forces ponctuelles --- $$$ #ifdef UTILISATION_MPI // cas d'un calcul //, les chargements aux noeuds sont calculés par le cpu 0 uniquement if (ParaGlob::Monde()->rank() == 0) #endif // on parcours le tableau tabTorseurPonct for (int i=1;i<= tabTorseurPonct.Taille();i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabTorseurPonct_i = tabTorseurPonct(i); const ReferenceNE & ref = ((ReferenceNE &) lesRef->Trouve(tabTorseurPonct_i.NomRef(),tabTorseurPonct_i.NomMaillage())); // dans le cas d'une courbe de charge globale on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabTorseurPonct_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabTorseurPonct_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabTorseurPonct_i.Pas_a_prendre_en_compte(Temps_courant()))) {int reftaille = ref.Taille(); int cas_assemblage = (assemb.Nb_cas_assemb()).n; // récup du cas d'assemblage // on renseigne les positions des points où vont s'appliquer les forces équivalentes au torseur Tableau < Coordonnee >& tab_P_i = tab_P(i); // pour simplier Tableau < Coordonnee >& t_force_i = t_force(i); // idem for (int nn =1; nn<= reftaille;nn++) // nn = position des num de noeud { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn)); tab_P_i(nn) = noe.Coord2(); }; // on calcul les efforts ponctuelles correspondant au torseur tabTorseurPonct_i.Force(tab_P_i,t_force_i); // maintenant on calcule l'action de chaque force ponctuelle for (int nn =1; nn<= reftaille;nn++) // nn = position des num de noeud { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn)); Coordonnee& force = t_force_i(nn); // pour simplifier // on regarde si jamais on a une dépendance globale à des grandeurs globales // ou au point + grandeurs globales éventuelles const string & nomF_charge = tabTorseurPonct_i.NomF_charge(); if(nomF_charge.length()) {Fonction_nD * pt_fonct = lesFonctionsnD->Trouve(nomF_charge); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales { // on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {force *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : force(3) *= tava(3); case 2 : force(2) *= tava(2); case 1 : force(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des torseurs de forces ponctuelles sur la ref: " << tabTorseurPonct_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; } else if (pt_fonct->Depend_M() != 0) // signifie qu'il y a une dépendance au point M { if ((pt_fonct->Depend_M()<0) && pt_fonct->Nom_variables().Taille() == dima) { // on vérifie également qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace const Coordonnee& M_tdt = noe.Coord2(); // comme pt_fonct->Depend_M()<0, cela signifie que la fonction dépend que de M // donc dans les paramètres d'appel on ne transmet que M Tableau tab_M(1,M_tdt); Tableau & tava = pt_fonct->Val_FnD_Evoluee(NULL,&tab_M,NULL); // pour simplifier if (pt_fonct->NbComposante() == 1) {force *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : force(3) *= tava(3); case 2 : force(2) *= tava(2); case 1 : force(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des torseurs de forces ponctuelles sur la ref: " << tabTorseurPonct_i.NomRef() << " en retour, le nombre de composante est different de 1 ou de la dimension de l'espace : " << dima << endl; Sortie(1); }; } else { if ((pt_fonct->Depend_M()>0)|| pt_fonct->Nom_variables().Taille() != dima) {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD utilisant les coordonnees du point" << " cas des forces ponctuelles sur la ref: "<< tabTorseurPonct_i.NomRef() << " actuellement la fonction doit avoir uniquement comme parametres, les coordonnees du point et eventuellement des variables globales "; cout << "\n fonction nD: "; pt_fonct->Affiche(); Sortie(1); }; }; } }; if ((noe.Pointeur_assemblage(X1,cas_assemblage) != -1) && (noe.En_service(X1)) && (noe.UneVariable(X1))) // les forces ponctuelles correspondent aux ddl X1 X2 ... for (int di=1;di<= dima;di++) {// maintenant deux cas: soit avec ou sans courbes de charge int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage); if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD // cas d'une grandeur fixe { vecglob (pt_adres) += coeff * (force(di));} else // cas d'une fonction nD sans courbe { vecglob (pt_adres) += (force(di))* tabTorseurPonct_i.Echelle_courbe();} } else { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabTorseurPonct_i.Echelle_courbe(); vecglob (pt_adres) += coeff_charge * (force(di)); }; }; }; }; }; #ifdef UTILISATION_MPI // on définit les conteneurs de passage d'info int proc_en_cours = ParaGlob::Monde()->rank(); int index_transfert; // index pour incrémenter dans les tableaux bool premier_passage = true; if (proc_en_cours != 0) { #endif // --- cas des forces surfaciques ---$$$ int tabsurfactaille = tabFsurfac.Taille(); if ( tabsurfactaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau tabFsurfac dans le repère absolu for (int i=1;i<= tabsurfactaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabFsurfac_i = tabFsurfac(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(tabFsurfac_i.NomRef(),tabFsurfac_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabFsurfac_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabFsurfac_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabFsurfac_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les surfaces associes Coordonnee vforce; // init par défaut int nbref = ref.Taille(); for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // récupération de la force de référence // on traite en fonction du fait que c'est un champ ou non switch (tabFsurfac_i.Champ()) { case 1: // c'est un champ de valeurs vforce = (tabFsurfac_i.Vect_de_coordonnee(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_ns_3 = lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,3)); vforce(3) = pt_courbe_force_ns_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_ns_2 = lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,2)); vforce(2) = pt_courbe_force_ns_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_ns_1 = lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,1)); vforce(1) = pt_courbe_force_ns_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ vforce = (tabFsurfac_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabFsurfac_i.NomF_charge(); Fonction_nD * pt_fonct=NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(tabFsurfac_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {vforce *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : vforce(3) *= tava(3); case 2 : vforce(2) *= tava(2); case 1 : vforce(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des forces surfaciques sur la ref: "<< tabFsurfac_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { vforce *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { vforce *= tabFsurfac_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFsurfac_i.Echelle_courbe(); vforce *= coeff_charge; }; #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression // appel du calcul du second membre correspondant à la charge surfacique Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier int num_face = ref.NumeroFA(ns); SM_transfert = elem.SM_charge_surfacique_E_tdt(vforce,pt_fonct,num_face,pa); int tyfront = 5; // init par défaut du type de frontière if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {//enum Enum_type_geom { POINT_G = 1 , LIGNE , SURFACE, VOLUME, RIEN_TYPE_GEOM }; tyfront = 2;} else // cas normale {tyfront = 3; };// la surface temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_SM(premier_passage,tyfront, num_face,SM_transfert, index_transfert,num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // appel du calcul du second membre correspondant à la charge surfacique Vecteur& SM = elem.SM_charge_surfacique_E_tdt(vforce,pt_fonct,ref.NumeroFA(ns),pa); // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée const ElFrontiere* elfront = NULL; // init if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns));} else // cas normale {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); };// la surface frontière assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage } #endif }; }; }; }; // --- cas des pressions ---$$$ int tabPresUniftaille = tabPresUnif.Taille(); if ( tabPresUniftaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau tabPresUnif for (int i=1;i<= tabPresUniftaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabPresUnif_i = tabPresUnif(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(tabPresUnif_i.NomRef(),tabPresUnif_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabPresUnif_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabPresUnif_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabPresUnif_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les surfaces associes int nbref = ref.Taille(); for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); double press_ac = 0.; // init // récupération de la pression de référence // on traite en fonction du fait que c'est un champ ou non switch (tabPresUnif_i.Champ()) { case 1: // c'est un champ de valeurs press_ac = (tabPresUnif_i.Val(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps { Courbe1D * pt_courbe_force_ns = lesCourbes1D->Trouve(tabPresUnif_i.Nom_val(ns)); press_ac = pt_courbe_force_ns->Valeur(Temps_courant()); break; } case 0: default:// ce n'est pas un champ press_ac = (tabPresUnif_i.Val(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabPresUnif_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(tabPresUnif_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {press_ac *= tava(1);} else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas de pression sur la ref: "<< tabPresUnif_i.NomRef() << " en retour, le nombre de composante est different de 1 : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { press_ac *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { press_ac *= tabPresUnif_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabPresUnif_i.Echelle_courbe(); press_ac *= coeff_charge; }; #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier int num_face = ref.NumeroFA(ns); SM_transfert = elem.SM_charge_pression_E_tdt(press_ac,pt_fonct,num_face,pa); int tyfront = 5; // init par défaut du type de frontière if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {//enum Enum_type_geom { POINT_G = 1 , LIGNE , SURFACE, VOLUME, RIEN_TYPE_GEOM }; tyfront = 2;} else // cas normale {tyfront = 3; };// la surface temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_SM(premier_passage,tyfront, num_face,SM_transfert, index_transfert,num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // appel du calcul du second membre correspondant à la charge de type pression // ici on récupère une adresse car il n'y a pas de crainte d'écraser un vecteur déjà en cours Vecteur& SM = elem.SM_charge_pression_E_tdt(press_ac,pt_fonct,ref.NumeroFA(ns),pa); // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée const ElFrontiere* elfront = NULL; // init if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns));} else // cas normale {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); };// la surface frontière assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage } #endif }; }; } }; // --- cas des pression unidirectionnelle --- // c'est-à-dire des forces qui sont appliquée sur la surface en gardant la même // direction localement pendant toute la transformation int PresUniDirtaille = PresUniDir.Taille(); if ( PresUniDirtaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau PresUniDir for (int i=1;i<= PresUniDirtaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& PresUniDir_i = PresUniDir(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(PresUniDir_i.NomRef(),PresUniDir_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(PresUniDir_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(PresUniDir_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(PresUniDir_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les surfaces associes int nbref = ref.Taille(); Coordonnee press_uni; // init par défaut for (int ns=1;ns<= nbref;ns++) {int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // récupération de la force de référence // on traite en fonction du fait que c'est un champ ou non switch (PresUniDir_i.Champ()) { case 1: // c'est un champ de valeurs press_uni = (PresUniDir_i.Vect_de_coordonnee(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_ns_3 = lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,3)); press_uni(3) = pt_courbe_force_ns_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_ns_2 = lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,2)); press_uni(2) = pt_courbe_force_ns_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_ns_1 = lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,1)); press_uni(1) = pt_courbe_force_ns_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ press_uni = (PresUniDir_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = PresUniDir_i.NomF_charge(); Fonction_nD * pt_fonct=NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(PresUniDir_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {press_uni *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : press_uni(3) *= tava(3); case 2 : press_uni(2) *= tava(2); case 1 : press_uni(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des des pression unidirectionnelle sur la ref: "<< PresUniDir_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { press_uni *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { press_uni *= PresUniDir_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * PresUniDir_i.Echelle_courbe(); press_uni *= coeff_charge; }; // appel du calcul du second membre correspondant à la charge de type pression #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier int num_face = ref.NumeroFA(ns); SM_transfert = elem.SM_charge_presUniDir_E_tdt(press_uni,pt_fonct,num_face,pa); int tyfront = 5; // init par défaut du type de frontière if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {//enum Enum_type_geom { POINT_G = 1 , LIGNE , SURFACE, VOLUME, RIEN_TYPE_GEOM }; tyfront = 2;} else // cas normale {tyfront = 3; };// la surface temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_SM(premier_passage,tyfront, num_face,SM_transfert, index_transfert,num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // appel du calcul du second membre correspondant à la charge de type pression // ici on récupère une adresse car il n'y a pas de crainte d'écraser un vecteur déjà en cours Vecteur& SM = elem.SM_charge_presUniDir_E_tdt(press_uni,pt_fonct,ref.NumeroFA(ns),pa); // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée const ElFrontiere* elfront = NULL; // init if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns));} else // cas normale {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); };// la surface frontière assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage } #endif }; }; }; }; // --- cas des forces lineique ---$$$ int tablineiquetaille = tabFlineique.Taille(); if ( tablineiquetaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau tabFlineique dans le repère absolu for (int i=1;i<= tablineiquetaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabFlineique_i = tabFlineique(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(tabFlineique_i.NomRef(),tabFlineique_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabFlineique_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabFlineique_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabFlineique_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les aretes associes int nbref = ref.Taille(); Coordonnee f_lin; // init par défaut for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // récupération de la force de référence // on traite en fonction du fait que c'est un champ ou non switch (tabFlineique_i.Champ()) { case 1: // c'est un champ de valeurs f_lin = (tabFlineique_i.Vect_de_coordonnee(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_ns_3 = lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,3)); f_lin(3) = pt_courbe_force_ns_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_ns_2 = lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,2)); f_lin(2) = pt_courbe_force_ns_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_ns_1 = lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,1)); f_lin(1) = pt_courbe_force_ns_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ f_lin = (tabFlineique_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabFlineique_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(tabFlineique_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {f_lin *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : f_lin(3) *= tava(3); case 2 : f_lin(2) *= tava(2); case 1 : f_lin(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des des forces lineiques sur la ref: "<< tabFlineique_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { f_lin *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { f_lin *= tabFlineique_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFlineique_i.Echelle_courbe(); f_lin *= coeff_charge; }; #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier int num_ligne = ref.NumeroFA(ns); SM_transfert =elem.SM_charge_lineique_E_tdt(f_lin,pt_fonct,num_ligne,pa); int tyfront = 2; // type de frontière ligne // si on est en axi c'est interprété comme une force surfacique // on ne change pas le type de frontière temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_SM(premier_passage,tyfront, num_ligne,SM_transfert, index_transfert,num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // appel du calcul du second membre correspondant à la charge lineique Vecteur& SM = elem.SM_charge_lineique_E_tdt(f_lin,pt_fonct,ref.NumeroFA(ns),pa); // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la ligne chargée const ElFrontiere* elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // l'arete frontière assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage } #endif }; }; } }; // --- cas des forces lineiques suiveuses ---$$$ int tablineiqueSuivtaille = tabFlineiqueSuiv.Taille(); if ( tablineiqueSuivtaille != 0) {// on parcours le tableau tabFlineiqueSuiv des forces exprimées dans le repère global #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif for (int i=1;i<= tablineiqueSuivtaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim > & tabFlineiqueSuiv_i = tabFlineiqueSuiv(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(tabFlineiqueSuiv_i.NomRef(),tabFlineiqueSuiv_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabFlineiqueSuiv_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabFlineiqueSuiv_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabFlineiqueSuiv_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les aretes associes int nbref = ref.Taille(); Coordonnee f_lin; // init par défaut for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // récupération de la force de référence // on traite en fonction du fait que c'est un champ ou non switch (tabFlineiqueSuiv_i.Champ()) { case 1: // c'est un champ de valeurs f_lin = (tabFlineiqueSuiv_i.Vect_de_coordonnee(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_ns_3 = lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,3)); f_lin(3) = pt_courbe_force_ns_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_ns_2 = lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,2)); f_lin(2) = pt_courbe_force_ns_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_ns_1 = lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,1)); f_lin(1) = pt_courbe_force_ns_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ f_lin = (tabFlineiqueSuiv_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabFlineiqueSuiv_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(tabFlineiqueSuiv_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {f_lin *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : f_lin(3) *= tava(3); case 2 : f_lin(2) *= tava(2); case 1 : f_lin(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des des forces lineiques suiveuses sur la ref: "<< tabFlineiqueSuiv_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { f_lin *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { f_lin *= tabFlineiqueSuiv_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFlineiqueSuiv_i.Echelle_courbe(); f_lin *= coeff_charge; } #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier int num_ligne = ref.NumeroFA(ns); SM_transfert = elem.SM_charge_lineique_Suiv_E_tdt(f_lin,pt_fonct,num_ligne,pa); int tyfront = 2; // type de frontière ligne // si on est en axi c'est interprété comme une force surfacique // on ne change pas le type de frontière temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_SM(premier_passage,tyfront, num_ligne,SM_transfert, index_transfert,num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // appel du calcul du second membre correspondant à la charge lineique // appel du calcul du second membre correspondant à la charge lineique // Vecteur f_lin = (tabFlineiqueSuiv_i.Coord()) * coeff ; Vecteur SM& = elem.SM_charge_lineique_Suiv_E_tdt(f_lin,pt_fonct,ref.NumeroFA(ns),pa); // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la ligne chargée const ElFrontiere* elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // l'arete frontière assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage }; #endif }; }; } }; // --- cas des forces volumiques ---$$$ int tabvoltaille = tabFvol.Taille(); if ( tabvoltaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau tabFvol dans le repère absolu for (int i=1;i<= tabvoltaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabFvol_i = tabFvol(i); const ReferenceNE & ref = ((ReferenceNE &) lesRef->Trouve(tabFvol_i.NomRef(),tabFvol_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabFvol_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabFvol_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabFvol_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments int nbref = ref.Taille(); Coordonnee f_vol; // init par défaut for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.Numero(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // récupération de la force de référence // on traite en fonction du fait que c'est un champ ou non switch (tabFvol_i.Champ()) { case 1: // c'est un champ de valeurs f_vol = (tabFvol_i.Vect_de_coordonnee(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_ns_3 = lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,3)); f_vol(3) = pt_courbe_force_ns_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_ns_2 = lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,2)); f_vol(2) = pt_courbe_force_ns_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_ns_1 = lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,1)); f_vol(1) = pt_courbe_force_ns_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ f_vol = (tabFvol_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabFvol_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(tabFvol_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {f_vol *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : f_vol(3) *= tava(3); case 2 : f_vol(2) *= tava(2); case 1 : f_vol(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des des forces lineiques sur la ref: "<< tabFvol_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { f_vol *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { f_vol *= tabFvol_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFvol_i.Echelle_courbe(); f_vol *= coeff_charge; }; // on indique si éventuellement on veut se référer au volume initial bool volume_finale = !(tabFvol_i.Attribut() == "sur_volume_initial_"); #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier int num_inter = 0; // ne sert pas ici SM_transfert = elem.SM_charge_volumique_E_tdt(f_vol,pt_fonct,pa,volume_finale); int tyfront = 4; // type d'un volume temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_SM(premier_passage,tyfront, num_inter,SM_transfert, index_transfert,num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // {// appel du calcul du second membre correspondant à la charge volumique Vecteur SM& = elem.SM_charge_volumique_E_tdt(f_vol,pt_fonct,pa,volume_finale); // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à l'élément assemb.AssemSM (vecglob,SM,elem.TableauDdl(),elem.Tab_noeud()); // assemblage } #endif }; }; } } // --- cas des des pressions hydrostatiques ---$$$ int PresHydrotaille = PresHydro.Taille(); if ( PresHydrotaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau PresHydro for (int i=1;i<= PresHydrotaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim > & PresHydro_i = PresHydro(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(PresHydro_i.NomRef(),PresHydro_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(PresHydro_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(PresHydro_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(PresHydro_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les surfaces associes int nbref = ref.Taille(); // on récupère le vecteur normal et le point de la surface libre const Coordonnee& N = PresHydro_i.Direction_N(); // récup de la direction N normée const Coordonnee& A = PresHydro_i.Point_M(); // récup du point de la surface libre double p_hydro_ref = (PresHydro_i.Val(1)); for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); double p_hydro = p_hydro_ref; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) { p_hydro *= coeff;} // cas d'une grandeur fixe else { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * PresHydro_i.Echelle_courbe(); p_hydro *= coeff_charge; } // appel du calcul du second membre correspondant à la charge de type pression bool sans_limitation=false; if (PresHydro_i.Attribut() == "sans_limitation_") sans_limitation= true; #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier int num_face = ref.NumeroFA(ns); SM_transfert = elem.SM_charge_hydrostatique_E_tdt (N,p_hydro,num_face,A,pa,sans_limitation) ; int tyfront = 5; // init par défaut du type de frontière if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {//enum Enum_type_geom { POINT_G = 1 , LIGNE , SURFACE, VOLUME, RIEN_TYPE_GEOM }; tyfront = 2;} else // cas normale {tyfront = 3; };// la surface temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_SM(premier_passage,tyfront, num_face,SM_transfert, index_transfert,num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // {Vecteur SM& = elem.SM_charge_hydrostatique_E_tdt (N,p_hydro,ref.NumeroFA(ns),A,pa,sans_limitation) ; // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée const ElFrontiere* elfront = NULL; // init if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns));} else // cas normale {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); };// la surface frontière assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage }; #endif }; }; } }; // --- cas des des efforts hydrodynamiques ---$$$ int coefHydroDynataille = coefHydroDyna.Taille(); if ( coefHydroDynataille != 0) {if (ParaGlob::AxiSymetrie()) // pas pris en compte {cout << "\n *** erreur : le chargement hydrodynamique en axisymetrique n'est pas pris en compte ! "; Sortie(1); }; #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau coefHydroDyna for (int i=1;i<= coefHydroDynataille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim > & coefHydroDyna_i = coefHydroDyna(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(coefHydroDyna_i.NomRef(),coefHydroDyna_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(coefHydroDyna_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(coefHydroDyna_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(coefHydroDyna_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les noeuds ou arretes ou surfaces associes suivant la dimension int nbref = ref.Taille(); // récup du poids volumique poids volumique = ro * g, qui n'est pas soumis à la courbe de charge double poidvol = (coefHydroDyna_i.Val(1)); for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // maintenant deux cas: soit avec ou sans courbes de charge double coeff_charge=1.; if(pt_courbe != NULL) { coeff_charge = (pt_courbe->Valeur(Temps_courant())) * coefHydroDyna_i.Echelle_courbe();}; #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier int num_front = ref.NumeroFA(ns); SM_transfert = elem.SM_charge_hydrodynamique_E_tdt (coefHydroDyna_i.Frot_fluid(),poidvol ,coefHydroDyna_i.Coef_aero_n(),num_front,coeff_charge ,coefHydroDyna_i.Coef_aero_t(),pa) ; int tyfront = 5; // init par défaut du type de frontière switch (ParaGlob::Dimension()) { case 3: {tyfront=3; break;};// la surface frontière case 2: {tyfront=2; break;};// la ligne frontière case 1: {tyfront=1; break;};// le point frontière }; temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_SM(premier_passage,tyfront, num_front,SM_transfert, index_transfert,num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // {// appel du calcul du second membre correspondant aux efforts hydrodynamique Vecteur SM& = elem.SM_charge_hydrodynamique_E_tdt(coefHydroDyna_i.Frot_fluid(),poidvol ,coefHydroDyna_i.Coef_aero_n(),ref.NumeroFA(ns),coeff_charge ,coefHydroDyna_i.Coef_aero_t(),pa) ; // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée const ElFrontiere* elfront; switch (ParaGlob::Dimension()) { case 3: {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); break;};// la surface frontière case 2: {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); break;};// la ligne frontière case 1: {elfront = elem.Frontiere_points(ref.NumeroFA(ns)); break;};// le point frontière }; assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage } #endif }; }; }; }; //-- fin de hydrodyna temps_cpu_chargement.Arret_du_comptage(); // fin comptage cpu // maintenant on s'occupe du cpu 0 #ifdef UTILISATION_MPI // arrivée ici on va indiquer au CPU 0 que c'est fini int tyfront = - proc_en_cours; // init par défaut du type de frontière // le signe - indique une fin // appel du transfert et mise à jour de l'index des tableaux Vecteur& SM_transfert = tabV_transfert(index_transfert); // ne sert pas ici Transfert_SM(premier_passage,tyfront, 0,SM_transfert, index_transfert,0,0 ); } else // donc ici cas ou on a le rank == 0 et on récolte { // récup d'un second membre temps_cpu_chargement.Mise_en_route_du_comptage(); // comptage cpu int nb_proc_terminer = 0; // permettra de terminer // on va boucler sur les éléments récupérés des différents process // jusqu'à ce que tous les éléments aient finis premier_passage = true; int index_transfert = 1; // on se sert du premier élément de tableau int ne,nbMail,tyfront,num_front,taille_SM; // inter pour la persistance int source; // le nb de la source émettante // int index_local_de_tab = 0; mpi::request & reqs1 = tab_reqs1(index_transfert); // pour simplifier mpi::request & reqs2 = tab_reqs2(index_transfert); // pour simplifier Vecteur& six_faux_entiers = tab_six_faux_entiers(index_transfert); // pour simplifier Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier // int i_proc = 1; // essai de réception avec le num du proc // Tableau < int > non_fin_proc(nb_proc-1,1); // essai while (nb_proc_terminer < (ParaGlob::Monde()->size()-1)) // gérer par les valeurs de tyfront { // on récupère un résultat de calcul temps_cpu_chargement.Arret_du_comptage(); // fin comptage cpu temps_transfert_court.Mise_en_route_du_comptage(); // comptage cpu reqs1 = six_faux_entiers.Irecup_MPI(mpi::any_source, 34); // if (non_fin_proc(i_proc)) // essai // {reqs1 = six_faux_entiers.Irecup_MPI(i_proc, 34); // essai // temps_transfert_court.Arret_du_comptage(); // comptage cpu // temps_attente.Mise_en_route_du_comptage(); // comptage cpu mpi::status stat = reqs1.wait(); temps_transfert_court.Arret_du_comptage(); // comptage cpu // temps_attente.Arret_du_comptage(); // comptage cpu temps_cpu_chargement.Mise_en_route_du_comptage(); // comptage cpu ne = (int) six_faux_entiers(1); nbMail = (int) six_faux_entiers(2); tyfront = (int) six_faux_entiers(3); num_front = (int) six_faux_entiers(4); taille_SM = (int) six_faux_entiers(5); index_transfert = six_faux_entiers(6); source = stat.source(); ////-------debug //cout << "\n debug Charge::ChargeSecondMembre_Ex_mecaSolid proc= "<< ParaGlob::Monde()->rank() //// << " recup CinqEntiers " << num_el_mail_tyfront_nufront << flush; // << " recup six_faux_entiers " << six_faux_entiers << flush; ////-------fin debug // int tyfront = num_el_mail_tyfront_nufront.trois; if (tyfront > 0) {Element& elem = lesMail->Element_LesMaille(nbMail,ne); // def de SM avec la bonne dimension SM_transfert.Change_taille(taille_SM);//num_el_mail_tyfront_nufront.cinq); // récupération des conteneurs ad hoc vecteur et raideur temps_cpu_chargement.Arret_du_comptage(); // fin comptage cpu temps_transfert_long.Mise_en_route_du_comptage(); // comptage cpu // reqs2 = SM_transfert.Irecup_MPI(source, 350+index_transfert); reqs2 = SM_transfert.Irecup_MPI(source, 35); reqs2.wait(); // on attend que le conteneur soit rempli temps_transfert_long.Arret_du_comptage(); // fin comptage cpu ////-------debug //cout << "\n debug Charge::ChargeSecondMembre_Ex_mecaSolid proc= "<< ParaGlob::Monde()->rank() // << " SM= " << SM << flush; ////-------fin debug // assemblage temps_cpu_chargement.Mise_en_route_du_comptage(); // comptage cpu switch (tyfront) { case 1: // cas d'un point frontière {const ElFrontiere* elfront =elem.Frontiere_points(num_front);// num_el_mail_tyfront_nufront.quatre); assemb.AssemSM (vecglob,SM_transfert,elfront->DdlElem_const(),elfront->TabNoeud_const()); break; } case 2: // cas d'une ligne {const ElFrontiere* elfront = elem.Frontiere_lineique(num_front);// num_el_mail_tyfront_nufront.quatre); assemb.AssemSM (vecglob,SM_transfert,elfront->DdlElem_const(),elfront->TabNoeud_const()); break; } case 3: // cas d'une surface {const ElFrontiere* elfront = elem.Frontiere_surfacique(num_front);//num_el_mail_tyfront_nufront.quatre); assemb.AssemSM (vecglob,SM_transfert,elfront->DdlElem_const(),elfront->TabNoeud_const()); break; } case 4: // cas d'un volume {assemb.AssemSM (vecglob,SM_transfert,elem.TableauDdl(),elem.Tab_noeud()); break; } default: {cout << "\n erreur*** Charge::ChargeSecondMembre_Ex_mecaSolid proc= "<< proc_en_cours << " le cas tyfront = " << tyfront << " n'est pas pris en compte pour l'instant !! "; Sortie(1); } ; }; } else if (tyfront < 0) {// on est dans le cas où un proc a terminé nb_proc_terminer++; // non_fin_proc(i_proc)=0; // essai } else {cout << "\n erreur*** Charge::ChargeSecondMembre_Ex_mecaSolid proc= "<< proc_en_cours << " tyfront ne doit jamais etre nul !! "; Sortie(1); } ; // }; // essai // i_proc++; // essai // if (i_proc > nb_proc-1) // essai // i_proc = 1; // essai }; temps_cpu_chargement.Arret_du_comptage(); // fin comptage cpu ////-------debug //cout << "\n debug Charge::ChargeSecondMembre_Ex_mecaSolid proc= "<< ParaGlob::Monde()->rank() // << " vecglob final = " << vecglob << flush; ////------fin debug }; #endif //#ifdef UTILISATION_MPI ////-------debug //cout << "\n debug Charge::ChargeSecondMembre_Ex_mecaSolid poc= "<< proc_en_cours // << " SM= " << vecglob << flush; //#else // cout << "\n debug Charge::ChargeSecondMembre_Ex_mecaSolid " // << " SM= "; vecglob.Affiche(); ////------fin debug //#endif // affichage éventuelle du second membre if (ParaGlob::NiveauImpression() >= 10) { string entete = " affichage du second membre apres chargement"; vecglob.Affichage_ecran(entete); }; } catch (ErrSortieFinale) // cas d'une direction voulue vers la sortie // on relance l'interuption pour le niveau supérieur { ErrSortieFinale toto; throw (toto); } catch ( ... ) { if (ParaGlob::NiveauImpression() >= 1) {cout << "\n\n ****!!** erreur dans le calcul du chargement impose ***"; if (ParaGlob::NiveauImpression() >= 4) cout << "\n Charge::ChargeSecondMembre_Ex_mecaSolid(.."; cout << flush; }; retour = false; }; temps_cpu_chargement.Arret_du_comptage(); #ifdef UTILISATION_MPI temps_attente.Mise_en_route_du_comptage(); ParaGlob::Monde()->barrier(); // synchronisation ici de tous les process temps_attente.Arret_du_comptage(); #endif return retour; }; // cas implicite avec également modification de la raideur quand le chargement dépend // de la géométrie, en fonction des paramètres de pa // la variable relative à la classe assemblage est passé en argument ce qui permet d'utiliser // éventuellement plusieurs type d'assemblage // de plus il faut prévoir un pointeur d'assemblage pour séparer le cas d'assemblage symétrique ou non // d'où la variable pointeur de fonction membre d'assemblage // retour false si problème de calcul bool Charge::ChargeSMembreRaideur_Im_mecaSolid (Assemblage & assemb,LesMaillages * lesMail,LesReferences* lesRef ,Vecteur& vecglob, Mat_abstraite& matglob ,void (Assemblage::* assembMat) // le pointeur de fonction (Mat_abstraite & matglob,const Mat_abstraite & matloc, const DdlElement& tab_ddl,const Tableau&tab_noeud) ,const ParaAlgoControle & pa ,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD) { temps_cpu_chargement.Mise_en_route_du_comptage(); bool retour = true; // init par défaut try { #ifndef UTILISATION_MPI // pour tous les maillages, pour tous les éléments on effectue une initialisation éventuelle int nb_mail = lesMail->NbMaillage(); for (int imail = 1;imail<=nb_mail;imail++) { int nb_ele = lesMail->Nombre_element(imail); for (int iele =1;iele<=nb_ele;iele++) lesMail->Element_LesMaille(imail,iele).Initialisation_avant_chargement(); }; #else // cas d'un calcul multi-CPU , si cpu 0, idem cas non parallèle if (ParaGlob::Monde()->rank() == 0) {int nb_mail = lesMail->NbMaillage(); for (int imail = 1;imail<=nb_mail;imail++) { int nb_ele = lesMail->Nombre_element(imail); for (int iele =1;iele<=nb_ele;iele++) lesMail->Element_LesMaille(imail,iele).Initialisation_avant_chargement(); }; } else {// on va parcourir les éléments associés au cpu // on récupère la distribution d'éléments concernant le cpu en cours // si c'est le cpu 0 on initialise tous les éléments const Tableau < list >& tab_list_elem_cpu = ParaGlob::param->const_List_element_CPU_en_cours(); // la taille est identique à nbMailMax, sauf si c'est le cpu 0, là le tableau est vide et il n'y // aura pas de boucle, (ou plutôt on sortira directement de la boucle // le cas CPU 0 est traité à part int nb_mail_distrib = tab_list_elem_cpu.Taille(); for (int imail =1; imail<= nb_mail_distrib; imail++) {const list & list_elem_cpu= tab_list_elem_cpu(imail); // pour simplifier // on balaie les éléments nécessaires list ::const_iterator il,ilfin=list_elem_cpu.end(); for (il = list_elem_cpu.begin();il != ilfin;il++) { int iele = (*il); lesMail->Element_LesMaille(imail,iele).Initialisation_avant_chargement(); }; }; }; #endif // récup bool avec_raid = pa.Var_charge_externe(); int cas_assemblage = (assemb.Nb_cas_assemb()).n; // récup du cas d'assemblage int posi = Id_nom_ddl("X1") -1; int dim = ParaGlob::Dimension(); // on regarde si l'espace de travail est axi-symétrique, si oui on utilise une dimension réduite int dima=dim; // c'est dima qui est utiliser pour la mise en place des efforts ponctuelles if (ParaGlob::AxiSymetrie()) dima--; // $$$--- cas des forces ponctuelles ---$$$ #ifdef UTILISATION_MPI // cas d'un calcul //, les chargements aux noeuds sont calculés par le cpu 0uniquement if (ParaGlob::Monde()->rank() == 0) #endif // on parcours le tableau tabPonctuel // ici il n'y a pas d'implication sur la raideur for (int i=1;i<= tabPonctuel.Taille();i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabPonctuel_i = tabPonctuel(i); //tabPonctuel_i.Affiche();cout << endl; const ReferenceNE & ref = ((ReferenceNE &) lesRef->Trouve(tabPonctuel_i.NomRef(),tabPonctuel_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabPonctuel_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabPonctuel_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabPonctuel_i.Pas_a_prendre_en_compte(Temps_courant()))) {int reftaille = ref.Taille(); for (int nn =1; nn<= reftaille;nn++) // nn = position des num de noeud { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn)); Coordonnee force; // init // on traite en fonction du fait que c'est un champ ou non switch (tabPonctuel_i.Champ()) { case 1: // c'est un champ de valeurs force = (tabPonctuel_i.Vect_de_coordonnee(nn)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_nn_3 = lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,3)); force(3) = pt_courbe_force_nn_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_nn_2 = lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,2)); force(2) = pt_courbe_force_nn_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_nn_1 = lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,1)); force(1) = pt_courbe_force_nn_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ force = (tabPonctuel_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabPonctuel_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(nomF_charge); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales { // on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {force *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : force(3) *= tava(3); case 2 : force(2) *= tava(2); case 1 : force(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; } else if (pt_fonct->Depend_M() != 0) // signifie qu'il y a une dépendance au point M { if ((pt_fonct->Depend_M()<0) && pt_fonct->Nom_variables().Taille() == dima) { // on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace const Coordonnee& M_tdt = noe.Coord2(); Tableau tab_M(1,M_tdt); Tableau & tava = pt_fonct->Val_FnD_Evoluee(NULL,&tab_M,NULL); // pour simplifier if (pt_fonct->NbComposante() == 1) {force *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : force(3) *= tava(3); case 2 : force(2) *= tava(2); case 1 : force(1) *= tava(1); ////----- debug --- //cout << "\n debug Charge::ChargeSMembreRaideur_Im_mecaSolid ** "; //cout << "\n tava= " << tava << endl; // ////--- fin debug --- }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef() << " en retour, le nombre de composante est different de 1 ou de la dimension de l'espace : " << dima << endl; Sortie(1); }; } else { if ((pt_fonct->Depend_M()>0)|| pt_fonct->Nom_variables().Taille() != dima) {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD utilisant les coordonnees du point" << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef() << " actuellement la fonction doit avoir uniquement comme parametres, les coordonnees du point et eventuellement des variables globales "; cout << "\n fonction nD: "; pt_fonct->Affiche(); Sortie(1); }; }; }; }; // les efforts s'appliques sur les ddl de position, il faut donc que X1 soit actif // pour le cas d'assemblage considéré et en activité if ((noe.Pointeur_assemblage(X1,cas_assemblage) != -1) && (noe.En_service(X1)) && (noe.UneVariable(X1))) // les forces ponctuelles correspondent aux ddl 1 2 ... for (int di=1;di<= dima;di++) {// maintenant deux cas: soit avec ou sans courbes de charge int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage); if (pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD // cas d'une grandeur fixe { vecglob (pt_adres) += coeff * (force(di));} else // cas d'une fonction nD sans courbe { vecglob (pt_adres) += (force(di))* tabPonctuel_i.Echelle_courbe();} } else // cas avec courbe, avec ou sans fonction c'est pareil { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabPonctuel_i.Echelle_courbe(); vecglob (pt_adres) += coeff_charge * (force(di)); }; }; }; }; }; // $$$--- cas des torseurs de forces ponctuelles ---$$$ #ifdef UTILISATION_MPI // cas d'un calcul //, les chargements aux noeuds sont calculés par le cpu 0 uniquement if (ParaGlob::Monde()->rank() == 0) #endif for (int i=1;i<= tabTorseurPonct.Taille();i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabTorseurPonct_i = tabTorseurPonct(i); //tabTorseurPonct_i.Affiche();cout << endl; const ReferenceNE & ref = ((ReferenceNE &) lesRef->Trouve(tabTorseurPonct_i.NomRef(),tabTorseurPonct_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabTorseurPonct_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabTorseurPonct_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabTorseurPonct_i.Pas_a_prendre_en_compte(Temps_courant()))) {int reftaille = ref.Taille(); // on renseigne les positions des points où vont s'appliquer les forces équivalentes au torseur Tableau < Coordonnee >& tab_P_i = tab_P(i); // pour simplier Tableau < Coordonnee >& t_force_i = t_force(i); // idem for (int nn =1; nn<= reftaille;nn++) // nn = position des num de noeud { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn)); tab_P_i(nn) = noe.Coord2(); }; // on calcul les efforts ponctuelles correspondant au torseur tabTorseurPonct_i.Force(tab_P_i,t_force_i); // maintenant on calcule l'action de chaque force ponctuelle for (int nn =1; nn<= reftaille;nn++) // nn = position des num de noeud { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn)); Coordonnee& force = t_force_i(nn); // pour simplifier // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabTorseurPonct_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(nomF_charge); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales { // on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {force *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : force(3) *= tava(3); case 2 : force(2) *= tava(2); case 1 : force(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des torseurs de forces ponctuelles sur la ref: "<< tabTorseurPonct_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; } else if (pt_fonct->Depend_M() != 0) // signifie qu'il y a une dépendance au point M { if ((pt_fonct->Depend_M()<0) && pt_fonct->Nom_variables().Taille() == dima) { // on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace const Coordonnee& M_tdt = noe.Coord2(); Tableau tab_M(1,M_tdt); Tableau & tava = pt_fonct->Val_FnD_Evoluee(NULL,&tab_M,NULL); // pour simplifier if (pt_fonct->NbComposante() == 1) {force *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : force(3) *= tava(3); case 2 : force(2) *= tava(2); case 1 : force(1) *= tava(1); ////----- debug --- //cout << "\n debug Charge::ChargeSMembreRaideur_Im_mecaSolid ** "; //cout << "\n tava= " << tava << endl; // ////--- fin debug --- }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des torseurs de forces ponctuelles sur la ref: "<< tabTorseurPonct_i.NomRef() << " en retour, le nombre de composante est different de 1 ou de la dimension de l'espace : " << dima << endl; Sortie(1); }; } else { if ((pt_fonct->Depend_M()>0)|| pt_fonct->Nom_variables().Taille() != dima) {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD utilisant les coordonnees du point" << " cas des des torseurs de ponctuelles sur la ref: "<< tabTorseurPonct_i.NomRef() << " actuellement la fonction doit avoir uniquement comme parametres, les coordonnees du point et eventuellement des variables globales "; cout << "\n fonction nD: "; pt_fonct->Affiche(); Sortie(1); }; }; }; }; // les efforts s'appliques sur les ddl de position, il faut donc que X1 soit actif // pour le cas d'assemblage considéré et en activité if ((noe.Pointeur_assemblage(X1,cas_assemblage) != -1) && (noe.En_service(X1)) && (noe.UneVariable(X1))) // les forces ponctuelles correspondent aux ddl 1 2 ... for (int di=1;di<= dima;di++) {// maintenant deux cas: soit avec ou sans courbes de charge int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage); if (pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD // cas d'une grandeur fixe { vecglob (pt_adres) += coeff * (force(di));} else // cas d'une fonction nD sans courbe { vecglob (pt_adres) += (force(di))* tabTorseurPonct_i.Echelle_courbe();} } else // cas avec courbe, avec ou sans fonction c'est pareil { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabTorseurPonct_i.Echelle_courbe(); vecglob (pt_adres) += coeff_charge * (force(di)); }; }; }; }; }; #ifdef UTILISATION_MPI // on définit les conteneurs de passage d'info int proc_en_cours = ParaGlob::Monde()->rank(); int index_transfert; // index pour incrémenter dans les tableaux bool premier_passage = true; if (proc_en_cours != 0) { #endif // $$$--- cas des forces surfaciques ---$$$ // il y a une implication éventuelle sur la raideur int tabsurfactaille = tabFsurfac.Taille(); if ( tabsurfactaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau tabFsurfac dans le repère absolu for (int i=1;i<= tabsurfactaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabFsurfac_i = tabFsurfac(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(tabFsurfac_i.NomRef(),tabFsurfac_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabFsurfac_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabFsurfac_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabFsurfac_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les surfaces associes int nbref = ref.Taille(); Coordonnee vforce; // init par défaut for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // récupération de la force de référence // on traite en fonction du fait que c'est un champ ou non switch (tabFsurfac_i.Champ()) { case 1: // c'est un champ de valeurs vforce = (tabFsurfac_i.Vect_de_coordonnee(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_ns_3 = lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,3)); vforce(3) = pt_courbe_force_ns_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_ns_2 = lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,2)); vforce(2) = pt_courbe_force_ns_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_ns_1 = lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,1)); vforce(1) = pt_courbe_force_ns_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ vforce = (tabFsurfac_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabFsurfac_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(tabFsurfac_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {vforce *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : vforce(3) *= tava(3); case 2 : vforce(2) *= tava(2); case 1 : vforce(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des forces surfaciques sur la ref: "<< tabFsurfac_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { vforce *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { vforce *= tabFsurfac_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFsurfac_i.Echelle_courbe(); vforce *= coeff_charge; }; // -- appel du calcul du second membre et de la raideur éventuelle // correspondant à la charge surfacique Element::ResRaid resu = elem.SMR_charge_surfacique_I(vforce,pt_fonct,ref.NumeroFA(ns),pa); #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression // appel du calcul du second membre correspondant à la charge surfacique Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier Mat_pleine& MAT_transfert = tabMat_transfert(index_transfert); Mat_pleine* pt_MAT = NULL; int num_face = ref.NumeroFA(ns); // on change de conteneur pour permettre de calculer plusieurs chargements même si le transfert n'est pas terminé SM_transfert = *(resu.res); if (avec_raid) {MAT_transfert = *(resu.raid); pt_MAT = &MAT_transfert; }; int tyfront = 5; // init par défaut du type de frontière if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {//enum Enum_type_geom { POINT_G = 1 , LIGNE , SURFACE, VOLUME, RIEN_TYPE_GEOM }; tyfront = 2;} else // cas normale {tyfront = 3; };// la surface temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_MatSm(premier_passage,tyfront, num_face,SM_transfert, pt_MAT,index_transfert, num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée const ElFrontiere* elfront = NULL; if (!ParaGlob::AxiSymetrie()) elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); // la surface frontière: cas normal else // cas axisymétrique elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // la surface arête générant une surface en rotation assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage du second membre // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); } #endif }; }; }; }; // $$$--- cas des pressions ---$$$ int tabPresUniftaille = tabPresUnif.Taille(); if ( tabPresUniftaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau tabPresUnif for (int i=1;i<= tabPresUniftaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabPresUnif_i = tabPresUnif(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(tabPresUnif_i.NomRef(),tabPresUnif_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabPresUnif_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabPresUnif_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabPresUnif_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les surfaces associes int nbref = ref.Taille(); for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); double press_ac = 0.; // init // récupération de la pression de référence // on traite en fonction du fait que c'est un champ ou non switch (tabPresUnif_i.Champ()) { case 1: // c'est un champ de valeurs press_ac = (tabPresUnif_i.Val(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps { Courbe1D * pt_courbe_force_ns = lesCourbes1D->Trouve(tabPresUnif_i.Nom_val(ns)); press_ac = pt_courbe_force_ns->Valeur(Temps_courant()); break; } case 0: default:// ce n'est pas un champ press_ac = (tabPresUnif_i.Val(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabPresUnif_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(tabPresUnif_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales { // on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {press_ac *= tava(1);} else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas de pression sur la ref: "<< tabPresUnif_i.NomRef() << " en retour, le nombre de composante est different de 1 : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { press_ac *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { press_ac *= tabPresUnif_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabPresUnif_i.Echelle_courbe(); press_ac *= coeff_charge; } // appel du calcul du second membre et de la raideur éventuelle // correspondant à la charge de type pression Element::ResRaid resu = elem.SMR_charge_pression_I(press_ac,pt_fonct,ref.NumeroFA(ns),pa); #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression // appel du calcul du second membre correspondant à la charge surfacique Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier Mat_pleine& MAT_transfert = tabMat_transfert(index_transfert); Mat_pleine* pt_MAT = NULL; int num_face = ref.NumeroFA(ns); // on change de conteneur pour permettre de calculer plusieurs chargements même si le transfert n'est pas terminé SM_transfert = *(resu.res); if (avec_raid) {MAT_transfert = *(resu.raid); pt_MAT = &MAT_transfert; }; int tyfront = 5; // init par défaut du type de frontière if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {//enum Enum_type_geom { POINT_G = 1 , LIGNE , SURFACE, VOLUME, RIEN_TYPE_GEOM }; tyfront = 2;} else // cas normale {tyfront = 3; };// la surface temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_MatSm(premier_passage,tyfront, num_face,SM_transfert, pt_MAT,index_transfert, num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée const ElFrontiere* elfront = NULL; if (!ParaGlob::AxiSymetrie()) elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); // la surface frontière: cas normal else // cas axisymétrique elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // la surface arête générant une surface en rotation assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); } #endif }; }; } } // $$$--- cas des pression unidirectionnelle ---$$$ // c'est-à-dire des forces qui sont appliquée sur la surface en gardant la même // direction localement pendant toute la transformation int PresUniDirtaille = PresUniDir.Taille(); if ( PresUniDirtaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau PresUniDir for (int i=1;i<= PresUniDirtaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& PresUniDir_i = PresUniDir(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(PresUniDir_i.NomRef(),PresUniDir_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(PresUniDir_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(PresUniDir_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(PresUniDir_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les surfaces associes int nbref = ref.Taille(); Coordonnee press_uni; // init par défaut for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // récupération de la force de référence // on traite en fonction du fait que c'est un champ ou non switch (PresUniDir_i.Champ()) { case 1: // c'est un champ de valeurs press_uni = (PresUniDir_i.Vect_de_coordonnee(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_ns_3 = lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,3)); press_uni(3) = pt_courbe_force_ns_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_ns_2 = lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,2)); press_uni(2) = pt_courbe_force_ns_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_ns_1 = lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,1)); press_uni(1) = pt_courbe_force_ns_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ press_uni = (PresUniDir_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = PresUniDir_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(PresUniDir_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {press_uni *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : press_uni(3) *= tava(3); case 2 : press_uni(2) *= tava(2); case 1 : press_uni(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des des pression unidirectionnelle sur la ref: "<< PresUniDir_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { press_uni *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { press_uni *= PresUniDir_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * PresUniDir_i.Echelle_courbe(); press_uni *= coeff_charge; } // appel du calcul du second membre et de la raideur éventuelle // correspondant à la charge de type pression // Vecteur pres_u = (PresUniDir_i.Coord()) * coeff; Element::ResRaid resu = elem.SMR_charge_presUniDir_I(press_uni,pt_fonct,ref.NumeroFA(ns),pa); #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression // appel du calcul du second membre correspondant à la charge surfacique Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier Mat_pleine& MAT_transfert = tabMat_transfert(index_transfert); Mat_pleine* pt_MAT = NULL; int num_face = ref.NumeroFA(ns); // on change de conteneur pour permettre de calculer plusieurs chargements même si le transfert n'est pas terminé SM_transfert = *(resu.res); if (avec_raid) {MAT_transfert = *(resu.raid); pt_MAT = &MAT_transfert; }; int tyfront = 5; // init par défaut du type de frontière if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {//enum Enum_type_geom { POINT_G = 1 , LIGNE , SURFACE, VOLUME, RIEN_TYPE_GEOM }; tyfront = 2;} else // cas normale {tyfront = 3; };// la surface temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_MatSm(premier_passage,tyfront, num_face,SM_transfert, pt_MAT,index_transfert, num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée const ElFrontiere* elfront = NULL; if (!ParaGlob::AxiSymetrie()) elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); // la surface frontière: cas normal else // cas axisymétrique elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // la surface arête générant une surface en rotation assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); } #endif }; }; } } // $$$--- cas des forces lineique ---$$$ int tablineiquetaille = tabFlineique.Taille(); if ( tablineiquetaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau tabFlineique dans le repère absolu for (int i=1;i<= tablineiquetaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabFlineique_i = tabFlineique(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(tabFlineique_i.NomRef(),tabFlineique_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabFlineique_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabFlineique_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabFlineique_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les aretes associes int nbref = ref.Taille(); Coordonnee f_lin; // init par défaut for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // récupération de la force de référence // on traite en fonction du fait que c'est un champ ou non switch (tabFlineique_i.Champ()) { case 1: // c'est un champ de valeurs f_lin = (tabFlineique_i.Vect_de_coordonnee(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_ns_3 = lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,3)); f_lin(3) = pt_courbe_force_ns_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_ns_2 = lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,2)); f_lin(2) = pt_courbe_force_ns_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_ns_1 = lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,1)); f_lin(1) = pt_courbe_force_ns_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ f_lin = (tabFlineique_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabFlineique_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(tabFlineique_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {f_lin *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : f_lin(3) *= tava(3); case 2 : f_lin(2) *= tava(2); case 1 : f_lin(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des des forces lineiques sur la ref: "<< tabFlineique_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { f_lin *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { f_lin *= tabFlineique_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFlineique_i.Echelle_courbe(); f_lin *= coeff_charge; } // appel du calcul du second membre et de la raideur éventuelle // correspondant à la charge lineique // Vecteur f_lin = (tabFlineique_i.Coord()) * coeff; Element::ResRaid resu = elem.SMR_charge_lineique_I(f_lin,pt_fonct,ref.NumeroFA(ns),pa) ; #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression // appel du calcul du second membre correspondant à la charge surfacique Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier Mat_pleine& MAT_transfert = tabMat_transfert(index_transfert); Mat_pleine* pt_MAT = NULL; int num_ligne = ref.NumeroFA(ns); // on change de conteneur pour permettre de calculer plusieurs chargements même si le transfert n'est pas terminé SM_transfert = *(resu.res); if (avec_raid) {MAT_transfert = *(resu.raid); pt_MAT = &MAT_transfert; }; int tyfront = 2; // type de frontière ligne // si on est en axi c'est interprété comme une force surfacique // on ne change pas le type de frontière temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_MatSm(premier_passage,tyfront, num_ligne,SM_transfert, pt_MAT,index_transfert, num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la ligne chargée const ElFrontiere* elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // l'arrête frontière assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); } #endif }; }; }; }; // $$$--- cas des forces lineiques suiveuses ---$$$ int tablineiqueSuivtaille = tabFlineiqueSuiv.Taille(); if ( tablineiqueSuivtaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau tabFlineiqueSuiv dans le repère global for (int i=1;i<= tablineiqueSuivtaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim > & tabFlineiqueSuiv_i = tabFlineiqueSuiv(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(tabFlineiqueSuiv_i.NomRef(),tabFlineiqueSuiv_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabFlineiqueSuiv_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabFlineiqueSuiv_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabFlineiqueSuiv_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les aretes associes int nbref = ref.Taille(); Coordonnee f_lin; // init par défaut for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // récupération de la force de référence // on traite en fonction du fait que c'est un champ ou non switch (tabFlineiqueSuiv_i.Champ()) { case 1: // c'est un champ de valeurs f_lin = (tabFlineiqueSuiv_i.Vect_de_coordonnee(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_ns_3 = lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,3)); f_lin(3) = pt_courbe_force_ns_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_ns_2 = lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,2)); f_lin(2) = pt_courbe_force_ns_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_ns_1 = lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,1)); f_lin(1) = pt_courbe_force_ns_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ f_lin = (tabFlineiqueSuiv_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabFlineiqueSuiv_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(tabFlineiqueSuiv_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {f_lin *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : f_lin(3) *= tava(3); case 2 : f_lin(2) *= tava(2); case 1 : f_lin(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des des forces lineiques suiveuses sur la ref: "<< tabFlineiqueSuiv_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { f_lin *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { f_lin *= tabFlineiqueSuiv_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFlineiqueSuiv_i.Echelle_courbe(); f_lin *= coeff_charge; } // appel du calcul du second membre et de la raideur éventuelle // correspondant à la charge lineique // Vecteur f_lin = (tabFlineiqueSuiv_i.Coord()) * coeff; Element::ResRaid resu = elem.SMR_charge_lineique_Suiv_I(f_lin,pt_fonct,ref.NumeroFA(ns),pa) ; #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression // appel du calcul du second membre correspondant à la charge surfacique Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier Mat_pleine& MAT_transfert = tabMat_transfert(index_transfert); Mat_pleine* pt_MAT = NULL; int num_ligne = ref.NumeroFA(ns); // on change de conteneur pour permettre de calculer plusieurs chargements même si le transfert n'est pas terminé SM_transfert = *(resu.res); if (avec_raid) {MAT_transfert = *(resu.raid); pt_MAT = &MAT_transfert; }; int tyfront = 2; // type de frontière ligne // si on est en axi c'est interprété comme une force surfacique // on ne change pas le type de frontière temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_MatSm(premier_passage,tyfront, num_ligne,SM_transfert, pt_MAT,index_transfert, num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la ligne chargée const ElFrontiere* elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // l'arrête frontière assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); } #endif }; }; }; }; // $$$--- cas des forces volumiques ---$$$ int tabvoltaille = tabFvol.Taille(); if ( tabvoltaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau tabFvol dans le repère absolu for (int i=1;i<= tabvoltaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim >& tabFvol_i = tabFvol(i); const ReferenceNE & ref = ((ReferenceNE &) lesRef->Trouve(tabFvol_i.NomRef(),tabFvol_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(tabFvol_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(tabFvol_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(tabFvol_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments int nbref = ref.Taille(); Coordonnee f_vol; // init par défaut for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.Numero(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // récupération de la force de référence // on traite en fonction du fait que c'est un champ ou non switch (tabFvol_i.Champ()) { case 1: // c'est un champ de valeurs f_vol = (tabFvol_i.Vect_de_coordonnee(ns)); break; case 2: // c'est un champ de fonctions // pour l'instant, le seul cas est celui des fonctions du temps // on récupère les courbes pour les dim composantes { switch (dima) { case 3 : {Courbe1D * pt_courbe_force_ns_3 = lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,3)); f_vol(3) = pt_courbe_force_ns_3->Valeur(Temps_courant()); } case 2 : {Courbe1D * pt_courbe_force_ns_2 = lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,2)); f_vol(2) = pt_courbe_force_ns_2->Valeur(Temps_courant()); } case 1 : {Courbe1D * pt_courbe_force_ns_1 = lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,1)); f_vol(1) = pt_courbe_force_ns_1->Valeur(Temps_courant()); } }; break; } case 0: default:// ce n'est pas un champ f_vol = (tabFvol_i.Vect_de_coordonnee(1)); break; }; // on regarde si jamais on a une dépendance à des grandeurs globales const string & nomF_charge = tabFvol_i.NomF_charge(); Fonction_nD * pt_fonct = NULL; if(nomF_charge.length()) {pt_fonct = lesFonctionsnD->Trouve(tabFvol_i.NomF_charge()); if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales {// on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace Tableau & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier if (pt_fonct->NbComposante() == 1) {f_vol *= tava(1);} else if (pt_fonct->NbComposante() == dima) {switch (dima) { case 3 : f_vol(3) *= tava(3); case 2 : f_vol(2) *= tava(2); case 1 : f_vol(1) *= tava(1); }; } else {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)" << " cas des des forces lineiques suiveuses sur la ref: "<< tabFvol_i.NomRef() << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : " << dima << endl; Sortie(1); }; // on désactive la fonction nD car elle a déjà été utilisée pt_fonct = NULL; }; }; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD { f_vol *= coeff;} // cas d'une grandeur fixe else // sinon il ne reste plus qu'à intégrer l'échelle { f_vol *= tabFvol_i.Echelle_courbe();}; } else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFvol_i.Echelle_courbe(); f_vol *= coeff_charge; } // appel du calcul du second membre et de la raideur éventuelle // correspondant à la charge volumique // on indique si éventuellement on veut se référer au volume initial bool volume_finale = !(tabFvol_i.Attribut() == "sur_volume_initial_"); // Vecteur f_vol = (tabFvol_i.Coord()) * coeff; Element::ResRaid resu = elem.SMR_charge_volumique_I(f_vol,pt_fonct,pa,volume_finale); #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression // appel du calcul du second membre correspondant à la charge surfacique Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier Mat_pleine& MAT_transfert = tabMat_transfert(index_transfert); Mat_pleine* pt_MAT = NULL; int num_inter = 0; // ne sert pas ici // on change de conteneur pour permettre de calculer plusieurs hargements même si le transfert n'est pas terminé SM_transfert = *(resu.res); if (avec_raid) {MAT_transfert = *(resu.raid); pt_MAT = &MAT_transfert; }; int tyfront = 4; // type d'un volume temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_MatSm(premier_passage,tyfront, num_inter,SM_transfert, pt_MAT,index_transfert, num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // { // assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à l'élément assemb.AssemSM (vecglob,*(resu.res),elem.TableauDdl(),elem.Tab_noeud()); // assemblage // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,*(resu.raid),elem.TableauDdl(),elem.Tab_noeud()); } #endif }; }; } }; // $$$--- cas des des pressions hydrostatiques ---$$$ int PresHydrotaille = PresHydro.Taille(); if ( PresHydrotaille != 0) { #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau PresHydro for (int i=1;i<= PresHydrotaille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim > & PresHydro_i = PresHydro(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(PresHydro_i.NomRef(),PresHydro_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(PresHydro_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(PresHydro_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(PresHydro_i.Pas_a_prendre_en_compte(Temps_courant()))) {// on récupère le vecteur normal et le point de la surface libre const Coordonnee& N = PresHydro_i.Direction_N(); // récup de la direction N normée const Coordonnee& A = PresHydro_i.Point_M(); // récup du point de la surface libre double p_hydro_ref = (PresHydro_i.Val(1)); // Maintenant les éléments et les surfaces associes int nbref = ref.Taille(); for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); double p_hydro = p_hydro_ref; // maintenant deux cas: soit avec ou sans courbes de charge if(pt_courbe == NULL) { p_hydro *= coeff;} // cas d'une grandeur fixe else { // cas d'une courbe de charge calcul de la valeur du coeff de charge double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * PresHydro_i.Echelle_courbe(); p_hydro *= coeff_charge; } // appel du calcul du second membre et de la raideur éventuelle // correspondant à la charge de type pression bool sans_limitation=false; if (PresHydro_i.Attribut() == "sans_limitation_") sans_limitation= true; Element::ResRaid resu = elem.SMR_charge_hydrostatique_I (N,p_hydro,ref.NumeroFA(ns),A,pa,sans_limitation) ; #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression // appel du calcul du second membre correspondant à la charge surfacique Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier int num_face = ref.NumeroFA(ns); Mat_pleine& MAT_transfert = tabMat_transfert(index_transfert); Mat_pleine* pt_MAT = NULL; int tyfront = 5; // init par défaut du type de frontière if (ParaGlob::AxiSymetrie()) // dans le cas axiSymétrique, les surfaces sont générées par les lignes, {//enum Enum_type_geom { POINT_G = 1 , LIGNE , SURFACE, VOLUME, RIEN_TYPE_GEOM }; tyfront = 2;} else // cas normale {tyfront = 3; };// la surface // on change de conteneur pour permettre de calculer plusieurs hargements même si le transfert n'est pas terminé SM_transfert = *(resu.res); if (avec_raid) {MAT_transfert = *(resu.raid); pt_MAT = &MAT_transfert; }; temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_MatSm(premier_passage,tyfront, num_face,SM_transfert, pt_MAT,index_transfert, num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // {// assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée const ElFrontiere* elfront = NULL; if (!ParaGlob::AxiSymetrie()) elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); // la surface frontière: cas normal else // cas axisymétrique elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // la surface arête générant une surface en rotation assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); }; #endif }; }; } }; // $$$--- cas des des efforts hydrodynamiques ---$$$ int coefHydroDynataille = coefHydroDyna.Taille(); if ( coefHydroDynataille != 0) {if (ParaGlob::AxiSymetrie()) // pas pris en compte {cout << "\n *** erreur : le chargement hydrodynamique en axisymetrique n'est pas pris en compte ! "; Sortie(1); }; #ifdef UTILISATION_MPI index_transfert=1; // index pour incrémenter dans les tableaux #endif // on parcours le tableau coefHydroDyna for (int i=1;i<= coefHydroDynataille;i++) { // recup de la reference correspondant au mot cle BlocCharge< BlocDdlLim > & coefHydroDyna_i = coefHydroDyna(i); const ReferenceAF & ref = ((ReferenceAF &) lesRef->Trouve(coefHydroDyna_i.NomRef(),coefHydroDyna_i.NomMaillage())); // dans le cas d'une courbe de charge on récupère l'adresse Courbe1D * pt_courbe = NULL; if(coefHydroDyna_i.NomCourbeCharge() != "") pt_courbe = lesCourbes1D->Trouve(coefHydroDyna_i.NomCourbeCharge()); // pas de calcul, si le temps est inférieur au temps actif // ou si d'une part le temps n'est pas actif et qu'au pas précédent // il n'était pas également actif on ne fait rien if (!(coefHydroDyna_i.Pas_a_prendre_en_compte(Temps_courant()))) {// Maintenant les éléments et les surfaces associes int nbref = ref.Taille(); // récup du poids volumique poids volumique = ro * g, qui n'est pas soumis à la courbe de charge double poidvol = (coefHydroDyna_i.Val(1)); for (int ns=1;ns<= nbref;ns++) { int num_elem = ref.NumeroElem(ns); int num_mail = ref.Nbmaille(); #ifdef UTILISATION_MPI // on ne continue que si l'élément est concerné if (ParaGlob::param->Element_concerner(num_mail,num_elem) ) {if (proc_en_cours != 0) { #endif // récupération de l'élément fini Element& elem = lesMail->Element_LesMaille(num_mail, num_elem); // maintenant deux cas: soit avec ou sans courbes de charge double coeff_charge=1.; if(pt_courbe != NULL) { coeff_charge = (pt_courbe->Valeur(Temps_courant())) * coefHydroDyna_i.Echelle_courbe();}; // appel du calcul du second membre correspondant aux efforts hydrodynamique Element::ResRaid resu = elem.SMR_charge_hydrodynamique_I(coefHydroDyna_i.Frot_fluid(),poidvol ,coefHydroDyna_i.Coef_aero_n(),ref.NumeroFA(ns),coeff_charge ,coefHydroDyna_i.Coef_aero_t(),pa) ; #ifdef UTILISATION_MPI // appel du calcul du second membre correspondant à la charge de type pression // appel du calcul du second membre correspondant à la charge surfacique Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier int num_front = ref.NumeroFA(ns); Mat_pleine& MAT_transfert = tabMat_transfert(index_transfert); Mat_pleine* pt_MAT = NULL; int tyfront = 5; // init par défaut du type de frontière switch (ParaGlob::Dimension()) { case 3: {tyfront=3; break;};// la surface frontière case 2: {tyfront=2; break;};// la ligne frontière case 1: {tyfront=1; break;};// le point frontière }; // on change de conteneur pour permettre de calculer plusieurs hargements même si le transfert n'est pas terminé SM_transfert = *(resu.res); if (avec_raid) {MAT_transfert = *(resu.raid); pt_MAT = &MAT_transfert; }; temps_cpu_chargement.Arret_du_comptage(); // appel du transfert et mise à jour de l'index des tableaux Transfert_MatSm(premier_passage,tyfront, num_front,SM_transfert, pt_MAT,index_transfert, num_mail,num_elem ); temps_cpu_chargement.Mise_en_route_du_comptage(); } }; #else // cas non // {// assemblage du second membre // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée const ElFrontiere* elfront; switch (ParaGlob::Dimension()) { case 3: {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); break;};// la surface frontière case 2: {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); break;};// la ligne frontière case 1: {elfront = elem.Frontiere_points(ref.NumeroFA(ns)); break;};// le point frontière }; assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); } #endif }; }; }; }; //-- fin de hydrodyna temps_cpu_chargement.Arret_du_comptage(); // fin comptage cpu // maintenant on s'occupe du cpu 0 #ifdef UTILISATION_MPI // arrivée ici on va indiquer au CPU 0 que c'est fini int tyfront = - proc_en_cours; // init par défaut du type de frontière // le signe - indique une fin // appel du transfert et mise à jour de l'index des tableaux Vecteur& SM_transfert = tabV_transfert(index_transfert); // ne sert pas ici Transfert_SM(premier_passage,tyfront, 0,SM_transfert, index_transfert,0,0 ); } else // donc ici cas ou on a le rank == 0 et on récolte { // récup d'un second membre et de la matrice temps_cpu_chargement.Mise_en_route_du_comptage(); // comptage cpu int nb_proc_terminer = 0; // permettra de terminer // on va boucler sur les éléments récupérés des différents process // jusqu'à ce que tous les éléments aient finis premier_passage = true; int index_transfert = 1; // on se sert du premier élément de tableau int ne,nbMail,tyfront,num_front,taille_SM; // inter pour la persistance int source; // le nb de la source émettante mpi::request & reqs1 = tab_reqs1(index_transfert); // pour simplifier mpi::request & reqs2 = tab_reqs2(index_transfert); // pour simplifier mpi::request & reqs3 = tab_reqs3(index_transfert); // pour simplifier Vecteur& six_faux_entiers = tab_six_faux_entiers(index_transfert); // pour simplifier Vecteur& SM_transfert = tabV_transfert(index_transfert); // pour simplifier Mat_pleine& MAT_transfert = tabMat_transfert(index_transfert); // idem while (nb_proc_terminer < (ParaGlob::Monde()->size()-1)) // gérer par les valeurs de tyfront { // on récupère un résultat de calcul temps_cpu_chargement.Arret_du_comptage(); // fin comptage cpu temps_transfert_court.Mise_en_route_du_comptage(); // comptage cpu reqs1 = six_faux_entiers.Irecup_MPI(mpi::any_source, 34); mpi::status stat = reqs1.wait(); temps_transfert_court.Arret_du_comptage(); // comptage cpu temps_cpu_chargement.Mise_en_route_du_comptage(); // comptage cpu ne = (int) six_faux_entiers(1); nbMail = (int) six_faux_entiers(2); tyfront = (int) six_faux_entiers(3); num_front = (int) six_faux_entiers(4); taille_SM = (int) six_faux_entiers(5); index_transfert = six_faux_entiers(6); source = stat.source(); ////-------debug //cout << "\n debug Charge::ChargeSMembreRaideur_Im_mecaSolid proc= "<< ParaGlob::Monde()->rank() //// << " recup CinqEntiers " << num_el_mail_tyfront_nufront << flush; // << " recup six_faux_entiers " << six_faux_entiers << flush; ////-------fin debug // int tyfront = num_el_mail_tyfront_nufront.trois; if (tyfront > 0) {Element& elem = lesMail->Element_LesMaille(nbMail,ne); // def de SM et MAT avec la bonne dimension SM_transfert.Change_taille(taille_SM); if (MAT_transfert.Nb_ligne () != taille_SM) MAT_transfert.Initialise(taille_SM,taille_SM); // récupération des conteneurs ad hoc vecteur et raideur temps_cpu_chargement.Arret_du_comptage(); // fin comptage cpu temps_transfert_long.Mise_en_route_du_comptage(); // comptage cpu reqs2 = SM_transfert.Irecup_MPI(source, 35); reqs2.wait(); // on attend que le conteneur soit rempli // idem pour la matrice reqs3 = MAT_transfert.Irecup_MPI(source,36); temps_transfert_long.Arret_du_comptage(); // fin comptage cpu ////-------debug //cout << "\n debug Charge::ChargeSMembreRaideur_Im_mecaSolid proc= "<< ParaGlob::Monde()->rank() // << " SM= " << SM << flush; ////-------fin debug // assemblage temps_cpu_chargement.Mise_en_route_du_comptage(); // comptage cpu switch (tyfront) { case 1: // cas d'un point frontière {const ElFrontiere* elfront =elem.Frontiere_points(num_front);// num_el_mail_tyfront_nufront.quatre); assemb.AssemSM (vecglob,SM_transfert,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,MAT_transfert,elfront->DdlElem_const(),elfront->TabNoeud_const()); break; } case 2: // cas d'une ligne {const ElFrontiere* elfront = elem.Frontiere_lineique(num_front);// num_el_mail_tyfront_nufront.quatre); assemb.AssemSM (vecglob,SM_transfert,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,MAT_transfert,elfront->DdlElem_const(),elfront->TabNoeud_const()); break; } case 3: // cas d'une surface {const ElFrontiere* elfront = elem.Frontiere_surfacique(num_front);//num_el_mail_tyfront_nufront.quatre); assemb.AssemSM (vecglob,SM_transfert,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,MAT_transfert,elfront->DdlElem_const(),elfront->TabNoeud_const()); break; } case 4: // cas d'un volume {assemb.AssemSM (vecglob,SM_transfert,elem.TableauDdl(),elem.Tab_noeud()); // assemblage de la raideur si nécessaire if (avec_raid) // appel de la fonction adoc, via le pointeur de fonction (assemb.*assembMat) (matglob,MAT_transfert,elem.TableauDdl(),elem.Tab_noeud()); break; } default: {cout << "\n erreur*** Charge::ChargeSMembreRaideur_Im_mecaSolid proc= "<< proc_en_cours << " le cas tyfront = " << tyfront << " n'est pas pris en compte pour l'instant !! "; Sortie(1); } ; }; } else if (tyfront < 0) {// on est dans le cas où un proc a terminé nb_proc_terminer++; } else {cout << "\n erreur*** Charge::ChargeSMembreRaideur_Im_mecaSolid proc= "<< proc_en_cours << " tyfront ne doit jamais etre nul !! "; Sortie(1); } ; }; temps_cpu_chargement.Arret_du_comptage(); // fin comptage cpu ////-------debug //cout << "\n debug Charge::ChargeSMembreRaideur_Im_mecaSolid proc= "<< ParaGlob::Monde()->rank() // << " vecglob final = " << vecglob << flush; ////------fin debug }; #endif //#ifdef UTILISATION_MPI ////-------debug //cout << "\n debug Charge::ChargeSecondMembre_Ex_mecaSolid poc= "<< proc_en_cours // << " SM= " << vecglob << flush; //#else // cout << "\n debug Charge::ChargeSecondMembre_Ex_mecaSolid " // << " SM= "; vecglob.Affiche(); ////------fin debug //#endif // affichage éventuelle de la matrice de raideur et du second membre if (ParaGlob::NiveauImpression() >= 10) { string entete = " affichage de la matrice de raideur apres chargement d'efforts externes "; matglob.Affichage_ecran(entete); entete = " affichage du second membre apres chargement d'efforts externes "; vecglob.Affichage_ecran(entete); }; } catch (ErrSortieFinale) // cas d'une direction voulue vers la sortie // on relance l'interuption pour le niveau supérieur { ErrSortieFinale toto; throw (toto); } catch ( ... ) { if (ParaGlob::NiveauImpression() >= 1) {cout << "\n *** erreur dans le calcul du chargement impose "; if (ParaGlob::NiveauImpression() >= 4) cout << "\n Charge::ChargeSMembreRaideur_Im_mecaSolid(.."; cout << flush; }; retour = false; }; temps_cpu_chargement.Arret_du_comptage(); return retour; };