// 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 "ElContact.h" #include "Droite.h" #include "ConstMath.h" #include "MathUtil.h" #include "ElemMeca.h" #include "Util.h" #include #include "TypeConsTens.h" #include "TypeQuelconqueParticulier.h" // --- calcul des puissances virtuelles développées par les efforts de contact ---------- // et eventuellement calcul de la raideur associé // -> explicite à tdt Vecteur* ElContact::SM_charge_contact() { if (Permet_affichage() > 5) {cout << "\n -- ElContact::SM_charge_contact: ";this->Affiche(1);} int contactType = ElContact::Recup_et_mise_a_jour_type_contact(); //cout << "\n -- debug ElContact::SM_charge_contact: contactType= " << contactType; // def de quelques dimensions pour simplifier int dima = ParaGlob::Dimension(); int tab_taille = tabNoeud.Taille(); int nbddl = ddlElement_assemblage->NbDdl(); energie_penalisation = 0.; energie_frottement.Inita(0.); // récup d'une place pour le résidu local et mise à jour de list_SM éventuellement RecupPlaceResidu(nbddl); Vecteur& SM_res = *residu; // pour simplifier Mise_a_jour_ddlelement_cas_solide_assemblage(); // mise à jour de ddl element et de cas_solide // la suite de la méthode ne fonctionne que pour les type 2 et 4 de contact, // dans le cas contraire on revient directement if (!((contactType == 2) || (contactType == 4) || (contactType == 41)|| (contactType == 42))) return residu; // ---- sinon .... cas du contact par pénalisation ----- // dans le cas 2: on calcul un facteur de pénalisation // dans le cas 4: on déduit le facteur de pénalisation à partir des forces internes // et de la condition cinématique de contact // on considère ici que le point de contact est déterminé et on s'occupe de déterminer les // efforts éventuels dus au contact ainsi que les énergies échangées sur le pas de temps // a priori on suppose que le contact est collant --> calcul des efforts de contact Coordonnee Noe_atdt = noeud->Coord2(); // récup de la position actuelle du noeud projectile // recup des dernières différentes informations calculées au niveau de la cinématique de contact Coordonnee M_impact; // le point d'impact, la normale Vecteur phii; RecupInfo(N,M_impact,phii,false); // en fait due aux différents algos de recherche de contact, normalement la position de contact est sauvegardée // dans l'élément, donc on utilise cette info, par contre on utilise RecupInfo pour calculer les autres infos : normale et phi // changement éventuel par une moyenne glissante des positions successive du noeud esclave // Noe_atdt est modifié éventuellement ChangeEnMoyGlissante(Noe_atdt); // calcul de la pénétration normale: deltaX=de la surface au point pénétré Coordonnee deltaX(dima); switch (contactType) { case 2: // cas d'une pénalisation classique case 41: case 42:// idem après le basculement du cas 4 { deltaX = Noe_atdt - M_impact; break; } case 4: // cas ou le noeud à été projeté sur la surface // normalement Noe_atdt est identique à M_impact { deltaX = M_noeud_tdt_avant_projection - Noe_atdt; break; } default: cout << "\n erreur, le cas contactType = " << contactType << " ne devrait pas exister ici !! "<< endl; Sortie(1); break; }; // le N sort de la surface maître, donc rentre dans la surface esclave double gap= N * deltaX; // gap est négatif quand il y a pénétration if ((Permet_affichage() > 7) && (contactType==4)) {double a = (deltaX - gap*N).Norme(); // on calcule pour vérifier cout << " diff_projX= "<< a; // l'erreur de projection normale }; // récupération du facteur de pénalisation adaptée au cas de calcul et à la géométrie double d_beta_gap=0.; // init de la dérivée du facteur de pénalisation / au gap // calcul du facteur de pénalisation si le contact est de type 2 if (contactType == 2) { penalisation = CalFactPenal(gap,d_beta_gap,contactType);} // sinon dans le cas 41 on utilise la valeur qui est en cours pondéré dans une moyenne glissante // avec le type 4 else if (contactType == 41) // après basculement du 4 on utilise la pénalisation { double essai_penalisation = CalFactPenal(gap,d_beta_gap,contactType); // maintenant on va calculer la moyenne glissante, de manière à lisser les variations ElContact::Moyenne_glissante_penalisation(penalisation, essai_penalisation); }; // si on est en 42, on conserve le facteur précédemment calculé gap_tdt = gap; // sauvegarde /* // debug //cout << "\n noeud: " << noeud->Num_noeud() <<" mailage:"<Num_Mail()<< " : gap= " << gap <<" \n Noe_atdt ";Noe_atdt.Affiche(); cout << " Noe_a0:"; noeud->Coord0().Affiche(); //cout << " Noe_at:"; noeud->Coord1().Affiche(); //cout << "\n M_impact:"; M_impact.Affiche(); cout << "\n deltaX:"; deltaX.Affiche(); //cout << "\n N:"; N.Affiche(); cout << " debug ElContact::SM_charge_contact()"<ParaAlgoControleActifs().Penetration_contact_maxi(); */ int typ_cal= ParaGlob::param->ParaAlgoControleActifs().TypeCalculPenalisationPenetration(); double borne_regularisation; if (fct_nD_contact.fct_nD_penetration_borne_regularisation != NULL) {borne_regularisation = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_borne_regularisation); } else {borne_regularisation = ParaGlob::param->ParaAlgoControleActifs().Penetration_borne_regularisation();} // pour l'instant je supprime le cas de la limitation du déplacement, car je pense qu'il faudrait y adjoindre les forces // de réaction associées sinon l'équilibre est faux !! donc on commente pour l'instant // méthode de limitation de la pénétration: opérationnelle si max_pene > 0 // si max_pene < 0 on intervient sur l'incrément de temps !! /* if (max_pene > 0) {if ((gap < (- max_pene)) && (Dabs(gap) > ConstMath::trespetit)) // on modifie la position du noeud pour tenir compte de la pénalisation // { Coordonnee nevez_posi(M_impact - 2 * max_pene * N); { Coordonnee nevez_posi(M_impact + (max_pene / deltaX.Norme()) * deltaX); noeud->Change_coord2(nevez_posi); gap = - max_pene; }; };*/ // --- calcul de la force normale // avec limitation éventuelle de l'intensité de la force de contact double intens_force = 0. ; // init bool force_limiter=false; Coordonnee F_N(dima); // init switch (contactType) { case 2: // cas d'une pénalisation classique case 41: case 42: // idem après basculement du cas 4 { intens_force = gap * penalisation; // négatif quand on a pénétration double max_force_noeud; if (fct_nD_contact.fct_nD_force_contact_noeud_maxi != NULL) {max_force_noeud = Valeur_fct_nD(fct_nD_contact.fct_nD_force_contact_noeud_maxi ,ElContact::Fct_nD_contact::tqi_fct_nD_force_contact_noeud_maxi ,ElContact::Fct_nD_contact::tqi_const_fct_nD_force_contact_noeud_maxi ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_force_contact_noeud_maxi); } else {max_force_noeud = ParaGlob::param->ParaAlgoControleActifs().Force_contact_noeud_maxi();} if (Dabs(intens_force) > max_force_noeud) {intens_force = max_force_noeud * Signe(intens_force); force_limiter = true; } // on recalcul la pénalisation associée if ((Dabs(gap) > ConstMath::trespetit) && force_limiter) {penalisation = intens_force / gap; if (Permet_affichage() > 5) {cout << " ** limitation penalisation= " << penalisation << " intensite force = " << intens_force << flush;}; }; // F_N c'est la force qui appuie sur le noeud esclave d'où le - car intens_force est négatif F_N = - N * intens_force; //(-gap * penalisation); break; } case 4: // cas ou le noeud à été projeté sur la surface // on récupère la force de contact via la puissance interne des éléments contenants // le noeud exclave (le calcul des forces internes doit déjà avoir été effectué) { TypeQuelconque_enum_etendu enu(FORCE_GENE_INT); // récupération du conteneur pour lecture uniquement const TypeQuelconque& typquel = noeud->Grandeur_quelconque(enu); const Grandeur_coordonnee& cofo = *((Grandeur_coordonnee*) typquel.Const_Grandeur_pointee()); // signe - car le contact doit s'opposer exactement à l'effort interne F_N = - cofo.ConteneurCoordonnee_const(); // récup de - la réaction if (Permet_affichage() > 6) cout << " F_N(force_interne)= " << F_N << " "; // en fait on ne va garder que la partie normale à la facette, ce qui laissera le glissement possible intens_force = - F_N * N; // 1) projection de la réaction F_N = - N * intens_force; // 2) recalcul uniquement de la partie normale // on va également tenter un calcul d'une pénalisation équivalente // penalisation = Dabs(intens_force / (Dabs(gap)));// + ConstMath::unpeupetit)) ; double max_pene; if (fct_nD_contact.fct_nD_penetration_contact_maxi != NULL) {max_pene = Dabs(Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_contact_maxi)); } else {max_pene = Dabs(ParaGlob::param->ParaAlgoControleActifs().Penetration_contact_maxi());} // double essai_penalisation; // variable intermédiaire // if (Dabs(gap) > max_pene) // {essai_penalisation = Dabs(intens_force) / Dabs(gap);} //(Dabs(gap) + ConstMath::unpeupetit)) ;} // else // sinon la pénétration est vraiment petite, on va cibler la valeur de max_pene // {essai_penalisation = Dabs(intens_force) / max_pene; // if (Permet_affichage() > 5) cout << "\n (penalisation limite) "; // }; // // on tient compte d'un facteur multiplicatif éventuelle // a priori c'est cette méthode qui est la plus performante (idem implicite) double essai_penalisation = Dabs(intens_force) / max_pene; // maintenant on va calculer la moyenne glissante ElContact::Moyenne_glissante_penalisation(penalisation, essai_penalisation); // if (Dabs(gap) > max_pene) // {penalisation = Dabs(intens_force) / Dabs(gap);} //(Dabs(gap) + ConstMath::unpeupetit)) ;} // else // sinon la pénétration est vraiment petite, on va cibler la valeur de max_pene // {penalisation = Dabs(intens_force) / max_pene; // if (Permet_affichage() > 5) cout << "\n (penalisation limite) "; // }; if (Permet_affichage() > 6) cout << " F_N(pur)= " << F_N << " "; break; } default: break; // déjà vu au dessus }; F_N_max = Dabs(intens_force); // sauvegarde if (Permet_affichage() > 5) {cout << " noeud: " << noeud->Num_noeud(); cout << "\n deltaX " << deltaX <<"\n Noe_atdt= " << Noe_atdt << " M_impact= " << M_impact <<"\n gap= " << gap << " N= " << N << " penalisation= " << penalisation << " intensite force = " << intens_force << flush; }; //cout << "\n SM_charge_contact: penalisation= "<< penalisation; // -- ici on a un traitement différent suivant typ_cal --- // Dans le cas du contact 4 on ne s'occupe pas du type de calcul de pénalisation if (contactType != 4) { switch (typ_cal) {case 1: case 2: case 3: // cas le plus simple : pas de calcul particulier if (gap > 0) // non ce n'est pas bizarre c'est le cas où le noeud décolle // c'est aussi le cas où on a un type de calcul de la pénalisation = 4 (il y aura une très légère force de collage pour // gap > 0 mais < borne_regularisation // par contre on annulle effectivement le résidu { //cout << "\n *** bizarre la penetration devrait etre negative ?? , on annule les forces de contact " // << "\n ElContact::SM_charge_contact(..."; force_contact=F_N; // on sauvegarde pour le traitement du décollement éventuel // au cas où on met à 0 les forces de réaction sur la facette int nb_n_facette = tabForce_cont.Taille(); for (int i=1; i<= nb_n_facette; i++) tabForce_cont(i).Zero(); if (Permet_affichage() > 6) { cout << "\n noeud: " << noeud->Num_noeud() << ": force collage non prise en compte " << force_contact << " gap(positif)= " << gap << ", Normale= " << N; cout << "\n deltaX " << deltaX <<"\n Noe_atdt= " << Noe_atdt << " M_impact= " << M_impact << " penalisation= " << penalisation << " intensite force = " << intens_force ; }; // on ne contribue pas aux réactions aux niveaux des noeuds car la valeur que l'on devrait ajoutée est nulle return residu; // les résidu + raideur sont nul à ce niveau de la méthode }; break; case 4: case 5:case 6: if (gap > Dabs(borne_regularisation)) // on est sortie de la zone d'accostage donc on neutralise le contact { force_contact=F_N; // on sauvegarde pour le traitement du décollement éventuel // au cas où on met à 0 les forces de réaction sur la facette int nb_n_facette = tabForce_cont.Taille(); for (int i=1; i<= nb_n_facette; i++) tabForce_cont(i).Zero(); if (Permet_affichage() >= 6) { cout << "\n noeud: " << noeud->Num_noeud() << ": force collage non prise en compte " << force_contact << " gap(positif et > borne regularisation)= " << gap << ", Normale= " << N; cout << "\n deltaX " << deltaX <<"\n Noe_atdt= " << Noe_atdt << " M_impact= " << M_impact << " penalisation= " << penalisation << " intensite force = " << intens_force ; }; // on ne contribue pas aux réactions aux niveaux des noeuds car la valeur que l'on devrait ajoutée est nulle return residu; // les résidu + raideur sont nul à ce niveau de la méthode }; break; case 7: case 8: // par rapport aux autres cas, on neutralise rien break; default : cout << "\n *** cas pas pris en compte !! \n SM_charge_contact() " << endl; Sortie(1); }; }; // -- fin du traitement différent suivant typ_cal --- /* ////////debug // if (ParaGlob::NiveauImpression() >= 6) //// if (noeud->Num_noeud()==113) // { //// cout << "\n debug : ElContact::SM_charge_contact(): noeud 113"; // cout << "\n noeud: " << noeud->Num_noeud(); // cout << "\n deltaX " << deltaX // <<"\n Noe_atdt= " << Noe_atdt << " M_impact= " << M_impact // <<"\n gap= " << gap << " N= " << N << " penalisation= " << penalisation << endl // ; // }; ////// fin debug //////debug // if (noeud->Num_noeud()==113) // { // cout << "\n debug : ElContact::SM_charge_contact(): noeud 113" // << "\n force " << F_N << " force_limiter "<PtEI(); // récup de l'élément fini CompFrotAbstraite* loiFrot = ((ElemMeca*) elemf)->LoiDeFrottement(); penalisation_tangentielle = 0.; // init bool force_T_limiter = false; // init if (loiFrot != NULL) { // --- cas de l'existance d'une loi de frottement // calcul du déplacement tangentiel Coordonnee M_tatdt; if (Mt.Dimension() != 0) {M_tatdt = (Mtdt-Mt);} // déplacement total // sinon Mt n'est pas encore définit donc on ne peut pas calculer le déplacement // il est donc laissé à 0. dep_tangentiel = M_tatdt- N * (M_tatdt*N); // on retire le deplacement normal éventuel // calcul de la force tangentielle en supposant tout d'abord un contact collant double penalisation_tangentielle; if (fct_nD_contact.fct_nD_penalisationTangentielle != NULL) {penalisation_tangentielle = Valeur_fct_nD(fct_nD_contact.fct_nD_penalisationTangentielle ,ElContact::Fct_nD_contact::tqi_fct_nD_penalisationTangentielle ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penalisationTangentielle ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penalisationTangentielle); } else {penalisation_tangentielle = ParaGlob::param->ParaAlgoControleActifs().PenalisationTangentielleContact();} // F_T =dep_tangentiel*(-ParaGlob::param->ParaAlgoControleActifs().PenalisationTangentielleContact()); F_T =dep_tangentiel*(-penalisation_tangentielle); // -- vérification en appelant le comportement en frottement ---- // tout d'abord on récupère le delta t pour le calcul de la vitesse double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant(); double unSurDeltat=0; if (Abs(deltat) >= ConstMath::trespetit) {unSurDeltat = 1./deltat;} else // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand { // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb if (unSurDeltat < 0) { cout << "\n le pas de temps est négatif !! " << "\n ElContact::SM_charge_contact(..."; }; unSurDeltat = ConstMath::tresgrand; }; // calcul des vitesses Coordonnee vit_T = dep_tangentiel * unSurDeltat; // appel de la loi de comportement: calcul de la force de frottement et des énergies échangées Coordonnee nevez_force_frot;EnergieMeca energie; bool glisse = loiFrot->Cal_explicit_contact_tdt(vit_T,N,F_N,F_T,energie,deltat,nevez_force_frot); // if (vit_T * nevez_force_frot > 0) // cout << "\n bizarre puissance de frottement positive !"; if (glisse) // on remplace la force de frottement si ça glisse sinon on conserve F_T F_T = nevez_force_frot; // --- calcul des énergies échangées energie_frottement += energie; if (Permet_affichage() > 6) { cout << "\n deltaX_T: " << dep_tangentiel <<"\n vit_T= " << vit_T << " glisse= " << glisse << " penalisation_tangentielle= " << penalisation_tangentielle; }; } else if ( cas_collant) { // --- cas d'un contact collant sans loi de frottement explicite // calcul du déplacement tangentiel par rapport au premier point de contact ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite const Coordonnee & M_C0 = elfro.Metrique()->PointM_tdt // ( stockage a tdt) (elfro.TabNoeud(),phi_theta_0); Coordonnee M_0atdt = (Noe_atdt-M_C0); // déplacement total à partir du point d'impact initial // transportée via l'interpolation: c-a-d les phi_theta_0 qui restent constant dep_tangentiel = M_0atdt- N * (M_0atdt*N); // on retire le deplacement normal éventuel double dep_T = dep_tangentiel.Norme(); // calcul de la pénalisation double d_beta_dep_T=0.; // init de la dérivée du facteur de pénalisation / au déplacement // pour l'instant ne sert pas penalisation_tangentielle = CalFactPenalTangentiel(dep_T,d_beta_dep_T); dep_T_tdt = dep_T; // sauvegarde // calcul de la force tangentielle en supposant un contact collant F_T = -dep_tangentiel * penalisation_tangentielle; // limitation éventuelle de la force double max_force_T_noeud; if (fct_nD_contact.fct_nD_force_tangentielle_noeud_maxi != NULL) {max_force_T_noeud = Valeur_fct_nD(fct_nD_contact.fct_nD_force_tangentielle_noeud_maxi ,ElContact::Fct_nD_contact::tqi_fct_nD_force_tangentielle_noeud_maxi ,ElContact::Fct_nD_contact::tqi_const_fct_nD_force_tangentielle_noeud_maxi ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_force_tangentielle_noeud_maxi); } else {max_force_T_noeud = ParaGlob::param->ParaAlgoControleActifs().Force_tangentielle_noeud_maxi();} double n_F_T = Dabs(penalisation_tangentielle * dep_T); if (n_F_T > max_force_T_noeud) {F_T *= max_force_T_noeud/n_F_T; force_T_limiter = true; }; if (Permet_affichage() > 6) { cout << "\n deltaX_T: " << dep_tangentiel << " penalisation_tangentielle= " << penalisation_tangentielle; }; // --- calcul des énergies échangées // ici il s'agit d'une énergie élastique, que l'on ajoute à la pénalisation normale // énergie élastique de pénalisation: énergie élastique = 1/2 k x^2 = 1/2 F x energie_penalisation += Dabs(0.5 * F_T * dep_tangentiel); // pour que le résultat soit > 0 energie_frottement.ChangeEnergieElastique(energie_penalisation); }; // --- calcul des puissances virtuelles F_T_max = F_T.Norme(); // sauvegarde force_contact = F_T + F_N; if (Permet_affichage() >= 5) { cout << "\n noeud: " << noeud->Num_noeud() << ": force contact " << force_contact ; cout << " ||F_N||= "<< F_N.Norme() << " ||F_T||= "<< F_T_max; if (force_limiter) cout << " force normale en limite " ; if (force_T_limiter) cout << " force tangentielle en limite " ; }; if (Permet_affichage() > 6) cout << "\n ==> force contact imposee au residu (noeud:" << noeud->Num_noeud() << " maillage:"<< noeud->Num_Mail() << ") " << force_contact; /* ////debug // if (noeud->Num_noeud()==207) // { // cout << "\n debug : ElContact::SM_charge_contact(): noeud 113"; // cout << "\n noeud: " << noeud->Num_noeud() << ": force contact " << force_contact ; // if (force_limiter) cout << " force en limite " ; // }; //// fin debug // -- debug //double inten = force_contact.Norme(); //if (inten > 0) // cout << "\n force_contact = " << force_contact << " (ElContact::SM_charge_contact())" << endl; // -- fin debug */ int posi = Id_nom_ddl("R_X1") -1; // préparation pour les réactions Ddl_enum_etendu& ddl_reaction_normale=Ddl_enum_etendu::Tab_FN_FT()(1); Ddl_enum_etendu& ddl_reaction_tangentielle=Ddl_enum_etendu::Tab_FN_FT()(2); // une petite manip pour le cas axisymétrique: dans ce cas les ddl à considérer pour le résidu et la raideur // doivent être uniquement x et y donc 2 ddl au lieu de 3. (rien ne dépend de z) int nb_ddl_en_var = dima; // init bool axi = false; if (ParaGlob::AxiSymetrie()) { nb_ddl_en_var -= 1; axi = true;}; // a- cas du noeud projectile int ism=1; if ((cas_solide == 0 ) || (cas_solide == 1)) // cas où le noeud est libre {for (int ix=1;ix<=nb_ddl_en_var;ix++,ism++) // pour le noeud projectile la vitesse virtuelle associé est { SM_res(ism)= force_contact(ix); // directement la vitesse du noeud // on abonde au niveau de la réaction au noeud noeud->Ajout_val_tdt(Enum_ddl(ix+posi),force_contact(ix)); }; } else // on s'occupe quand même des réactions, ce qui nous permettra de les récupérer même sur un solide {for (int ix=1;ix<=dima;ix++) // pour le noeud projectile la vitesse virtuelle associé est // on abonde au niveau de la réaction au noeud noeud->Ajout_val_tdt(Enum_ddl(ix+posi),force_contact(ix)); }; // dans le cas de l'axi il faut décaler ism pour la suite, ce qui correspond au z // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus // if (axi) ism++; // qui n'est pas pris en compte // on abonde au niveau de la réaction normale et tangentielle pour le noeud esclave noeud->ModifDdl_etendu(ddl_reaction_normale).Valeur() += -intens_force; double intensite_tangentielle = F_T.Norme(); noeud->ModifDdl_etendu(ddl_reaction_tangentielle).Valeur() += intensite_tangentielle; // b- cas de la surface cible if ((cas_solide == 0 ) || (cas_solide == 2)) // cas où la facette est libre {for (int ir=1;ir<=tab_taille-1;ir++) { Noeud* noe = tabNoeud(ir+1); // pour simplifier tabForce_cont(ir) = -force_contact*phii(ir); for (int ix=1;ix<=nb_ddl_en_var;ix++,ism++) { double force_imp = -force_contact(ix)*phii(ir); SM_res(ism)= force_imp; // on abonde au niveau de la réaction au noeud noe->Ajout_val_tdt(Enum_ddl(ix+posi),force_imp); }; // dans le cas de l'axi il faut décaler ism pour la suite, ce qui correspond au z // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus // if (axi) ism++; // qui n'est pas pris en compte // on abonde au niveau de la réaction normale et tangentielle pour les noeuds maîtres noe->ModifDdl_etendu(ddl_reaction_normale).Valeur() += -intens_force*phii(ir); noe->ModifDdl_etendu(ddl_reaction_tangentielle).Valeur() += intensite_tangentielle*phii(ir); }; } else // on s'occupe quand même des réactions, ce qui nous permettra de les récupérer même sur un solide {for (int ir=1;ir<=tab_taille-1;ir++) { Noeud* noe = tabNoeud(ir+1); // pour simplifier tabForce_cont(ir) = -force_contact*phii(ir); for (int ix=1;ix<=dima;ix++) { double force_imp = -force_contact(ix)*phii(ir); // on abonde au niveau de la réaction au noeud noe->Ajout_val_tdt(Enum_ddl(ix+posi),force_imp); }; // on abonde au niveau de la réaction normale et tangentielle pour les noeuds maîtres noe->ModifDdl_etendu(ddl_reaction_normale).Valeur() += -intens_force*phii(ir); noe->ModifDdl_etendu(ddl_reaction_tangentielle).Valeur() += intensite_tangentielle*phii(ir); }; }; if (Permet_affichage() >= 7) { cout << "\n =-=-=- bilan sur le calcul de l'element de contact: =-=-=- "; this->Affiche(); }; /* // retour ////debug // cout << "\n debug : ElContact::SM_charge_contact()"; //noeud->Affiche(); //for (int ir=1;ir<=tab_taille-1;ir++) // tabNoeud(ir+1)->Affiche(); // // cout << "\n residu " << *residu < implicite, Element::ResRaid ElContact::SM_K_charge_contact() { if (Permet_affichage() > 5) {cout << "\n -- SM_K_charge_contact: ";this->Affiche(1);}; // le contact 4 passe en 41 ou 42 quand on dépasse un nombre donné d'itérations par exemple int contactType = ElContact::Recup_et_mise_a_jour_type_contact(); // préparation du retour Element::ResRaid el;el.res=NULL;el.raid=NULL; energie_penalisation = 0.; energie_frottement.Inita(0.); Mise_a_jour_ddlelement_cas_solide_assemblage(); // mise à jour de ddl element et de cas_solide // dans le cas du modèle de contact cinématique, on n'a pas de contribution (pour l'instant) au second membre et à la raideur // cela pourrait changer si l'on considérait le frottement, mais pour l'instant ce n'est pas le cas // on retourne un pointeur null if (contactType == 1) return el; // def de quelques dimensions pour simplifier int dima = ParaGlob::Dimension(); int tab_taille = tabNoeud.Taille(); int nbddl = ddlElement_assemblage->NbDdl(); // récup d'une place pour le résidu local et mise à jour de list_SM éventuellement RecupPlaceResidu(nbddl); Vecteur& SM_res = *residu; // pour simplifier // idem pour la raideur RecupPlaceRaideur(nbddl); // ---- initialisation de SM et de la raideur SM_res.Zero();raideur->Initialise(0.0); el.res = residu; el.raid = raideur; // la suite de la méthode ne fonctionne que pour les type 2 et 4 de contact, // dans le cas contraire on revient directement if (!((contactType == 2) || (contactType == 4) || (contactType == 41) || (contactType == 42))) return el; // ---- sinon .... cas du contact par pénalisation ----- // dans le cas 2: on calcul un facteur de pénalisation // dans le cas 4: on déduit le facteur de pénalisation à partir des forces internes // et de la condition cinématique de contact // on considère ici que le point de contact est déterminé et on s'occupe de déterminer les // efforts éventuels dus au contact ainsi que les énergies échangées sur le pas de temps // a priori on suppose que le contact est collant --> calcul des efforts de contact Coordonnee Noe_atdt = noeud->Coord2(); // récup de la position actuelle du noeud projectile // recup du dernier des différentes informations Coordonnee M_impact; // le point d'impact, la normale Vecteur phii; ////debug // if (noeud->Num_noeud()== 2) // { // cout << "\n debug : ElContact::SM_K_charge_contact():" // << " noeud: " << noeud->Num_noeud() <<" mailage:"<Num_Mail() // << endl; // }; //// fin debug // //debug // if (Permet_affichage() > 4) // {cout << "\n debug ElContact::SM_K_charge_contact avant RecupInfo " ; // }; // //fin debug // d_T_pt: représente la variation du vecteur normale Tableau * d_T_pt = (RecupInfo(N,M_impact,phii,true)); // en fait due aux différents algos de recherche de contact, normalement la position de contact est sauvegardée // dans l'élément, donc on utilise cette info, par contre on utilise RecupInfo pour calculer les autres infos : normale et phi // M_impact et Mtdt sont différents uniquement dans le cas d'un restart // Mtdt est mis à jour pendant la recherche du contact M_impact = Mtdt; // changement éventuel par une moyenne glissante des positions successive du noeud esclave // Noe_atdt est modifié éventuellement ChangeEnMoyGlissante(Noe_atdt); // calcul de la pénétration normale Coordonnee deltaX(dima); switch (contactType) { case 2: // cas d'une pénalisation classique case 41 : case 42: // idem après basculement du 4 { deltaX = Noe_atdt - M_impact; break; } case 4: // cas ou le noeud à été projeté sur la surface // normalement Noe_atdt est identique à M_impact { deltaX = M_noeud_tdt_avant_projection - Noe_atdt; break; } default: cout << "\n *** erreur, le cas contactType = " << contactType << " ne devrait pas exister ici !! "<< endl; Sortie(1); break; }; double gap= N * deltaX; // ce qui constitue la fonction g(ddl) à annuler voir à minimiser if ((Permet_affichage() > 7) && (contactType==4)) {double a = (deltaX - gap*N).Norme(); // on calcule pour vérifier cout << " diff_projX= "<< a; // l'erreur de projection normale }; ////debug // if (noeud->Num_noeud()==56) // { // cout << "\n debug : ElContact::SM_K_charge_contact(): noeud 409" // << "\n deltaX " << deltaX // <<"\n Noe_atdt= " << Noe_atdt << " M_impact= " << M_impact // <<"\n gap= " << gap << " N= " << N << endl // ; // }; //// fin debug // récupération du facteur de pénalisation adaptée au cas de calcul et à la géométrie double d_beta_gap=0.; // init de la dérivée du facteur de pénalisation / au gap // calcul du facteur de pénalisation si le contact est de type 2 if (contactType == 2) {penalisation = CalFactPenal(gap,d_beta_gap,contactType);} // sinon dans le cas 41 ou 42 on utilise la valeur qui est en cours // il n'y a pas de différence entre 41 et 42 en implicite gap_tdt = gap; // sauvegarde // limitation au niveau du gap à 2 fois la pénétration maxi // double max_pene = ParaGlob::param->ParaAlgoControleActifs().Penetration_contact_maxi(); double borne_regularisation; if (fct_nD_contact.fct_nD_penetration_borne_regularisation != NULL) {borne_regularisation = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_borne_regularisation); } else {borne_regularisation = ParaGlob::param->ParaAlgoControleActifs().Penetration_borne_regularisation();} int typ_cal= ParaGlob::param->ParaAlgoControleActifs().TypeCalculPenalisationPenetration(); //debug // if (noeud->Num_noeud()==7) // {double max_pene = ParaGlob::param->ParaAlgoControleActifs().Penetration_contact_maxi(); // if ((gap_t > max_pene) && (gap_tdt > max_pene)) // {cout << "\n debug : ElContact::SM_K_charge_contact(): noeud 7" // << "\n gap_t " << gap_t << " max_pene= " << max_pene // <<"\n gap_tdt= " << gap_tdt << " M_impact= " << M_impact // <<"\n gap= " << gap << " N= " << N << endl; // }; // }; // fin debug // pour l'instant je supprime le cas de la limitation du déplacement, car je pense qu'il faudrait y adjoindre les forces // de réaction associées sinon l'équilibre est faux !! donc on commente pour l'instant /* // méthode de limitation de la pénétration: opérationnelle si max_pene > 0 // si max_pene < 0 on intervient sur l'incrément de temps !! if (max_pene > 0) {if ((typ_cal == 5) && (gap < (- 2 * max_pene))) // on modifie la position du noeud pour tenir compte de la pénalisation // { Coordonnee nevez_posi(M_impact - 2 * max_pene * N); { Coordonnee nevez_posi(M_impact + (2. * max_pene / deltaX.Norme()) * deltaX); noeud->Change_coord2(nevez_posi); }; };*/ // --- calcul de la force normale // limitation éventuelle de l'intensité de la force de contact double intens_force = 0. ; // init bool force_limiter=false; Coordonnee F_N(dima); // init double max_force_noeud; if (fct_nD_contact.fct_nD_force_contact_noeud_maxi != NULL) {max_force_noeud = Valeur_fct_nD(fct_nD_contact.fct_nD_force_contact_noeud_maxi ,ElContact::Fct_nD_contact::tqi_fct_nD_force_contact_noeud_maxi ,ElContact::Fct_nD_contact::tqi_const_fct_nD_force_contact_noeud_maxi ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_force_contact_noeud_maxi); } else {max_force_noeud = ParaGlob::param->ParaAlgoControleActifs().Force_contact_noeud_maxi();} switch (contactType) { case 2: // cas d'une pénalisation classique case 41: case 42: // idem après basculement du cas 4 { intens_force = gap * penalisation; // négatif quand on a pénétration if (Dabs(intens_force) > max_force_noeud) {intens_force = max_force_noeud * Signe(intens_force); force_limiter = true; }; // on recalcul la pénalisation associée if ((Dabs(gap) > ConstMath::trespetit) && force_limiter) {penalisation = intens_force / gap; if (Permet_affichage() > 5) {cout << " ** limitation penalisation= " << penalisation << " intensite force = " << intens_force << flush;}; }; // F_N c'est la force qui appuie sur le noeud esclave d'où le - car intens_force est négatif F_N = - N * intens_force; //(-gap * penalisation); break; } case 4: // cas ou le noeud à été projeté sur la surface // on récupère la force de contact via la puissance interne des éléments contenants // le noeud exclave (le calcul des forces internes doit déjà avoir été effectué) { TypeQuelconque_enum_etendu enu(FORCE_GENE_INT); // récupération du conteneur pour lecture uniquement const TypeQuelconque& typquel = noeud->Grandeur_quelconque(enu); const Grandeur_coordonnee& cofo = *((Grandeur_coordonnee*) typquel.Const_Grandeur_pointee()); // signe - car le contact doit s'opposer exactement à l'effort interne F_N = - cofo.ConteneurCoordonnee_const(); // récup de - la réaction if (Permet_affichage() > 6) cout << " F_N(force_interne)= " << F_N << " "; // en fait on ne va garder que la partie normale à la facette, ce qui laissera le glissement possible intens_force = - F_N * N; // 1) projection de la réaction F_N = - N * intens_force; // 2) recalcul uniquement de la partie normale // on va également tenter un calcul d'une pénalisation équivalente double max_pene; if (fct_nD_contact.fct_nD_penetration_contact_maxi != NULL) {max_pene = Dabs(Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_contact_maxi)); } else {max_pene = Dabs(ParaGlob::param->ParaAlgoControleActifs().Penetration_contact_maxi());} // l'idée est de viser la pénétration max_pene, car en pénalisation, le fait // de conserver la pénétration et la force constatées, ne va pas diminuer la pénétration // or on aimerait une pénétration petite !! a priori celle = max_pene // // si cette grandeur est plus petite que gap // // si par contre max_pene est grand, (ou bien plus grand que gap) on utilise gap // penalisation = Dabs(intens_force) / (MiN(max_pene, Dabs(gap))); // if (Dabs(gap) > max_pene) // {penalisation = Dabs(intens_force) / max_pene; } // else // {penalisation = Dabs(intens_force) / Dabs(gap);} // a priori c'est cette méthode qui est la plus performante penalisation = Dabs(intens_force) / max_pene; // if (Dabs(gap) > max_pene) // {penalisation = Dabs(intens_force) / Dabs(gap);} //(Dabs(gap) + ConstMath::unpeupetit)) ;} // else // sinon la pénétration est vraiment petite, on va cibler la valeur de max_pene // {penalisation = Dabs(intens_force) / max_pene; // if (Permet_affichage() > 5) cout << "\n (penalisation limite) "; // }; if (Permet_affichage() > 6) cout << " F_N(pur)= " << F_N << " "; break; } default: break; // déjà vu au dessus }; F_N_max = Dabs(intens_force); // sauvegarde if (Permet_affichage() > 5) {cout << " noeud: " << noeud->Num_noeud(); cout << "\n deltaX " << deltaX <<"\n Noe_atdt= " << Noe_atdt << " M_impact= " << M_impact <<"\n gap= " << gap << " N= " << N << " penalisation= " << penalisation << " intensite force = " << intens_force << flush; }; //cout << "\n SM_K_charge_contact: penalisation= "<< penalisation << " gap= " << gap << " FN= " << F_N ; //debug // if ((F_N_max > 5000000) || (noeud->Num_noeud()==2)) // if (penalisation == 0.) // if (noeud->Num_noeud()==2) // { // cout << "\n debug : ElContact::SM_K_charge_contact(): noeud " << noeud->Num_noeud() // << " zone: " << num_zone_contact // << "\n force " << F_N << " force_limiter "< 6) { cout << " noeud: " << noeud->Num_noeud() << ": force collage non prise en compte " << force_contact << " gap(positif)= " << gap << ", Normale= " << N; cout << "\n deltaX " << deltaX <<"\n Noe_atdt= " << Noe_atdt << " M_impact= " << M_impact << " penalisation= " << penalisation << " intensite force = " << intens_force ; }; // on ne contribue pas aux réactions aux niveaux des noeuds car la valeur que l'on devrait ajoutée est nulle return el; // les résidu + raideur sont nul à ce niveau de la méthode }; break; case 4: case 5: case 6: if (gap > Dabs(borne_regularisation)) // on est sortie de la zone d'accostage donc on neutralise le contact { force_contact=F_N; // on sauvegarde pour le traitement du décollement éventuel // au cas où on met à 0 les forces de réaction sur la facette int nb_n_facette = tabForce_cont.Taille(); for (int i=1; i<= nb_n_facette; i++) tabForce_cont(i).Zero(); if (Permet_affichage() >= 6) { cout << " noeud: " << noeud->Num_noeud() << ": force collage non prise en compte " << force_contact << " gap(positif et > borne regularisation)= " << gap << ", Normale= " << N; cout << "\n deltaX " << deltaX <<"\n Noe_atdt= " << Noe_atdt << " M_impact= " << M_impact << " penalisation= " << penalisation << " intensite force = " << intens_force ; }; // on ne contribue pas aux réactions aux niveaux des noeuds // car la valeur que l'on devrait ajoutée est nulle return el; // les résidu + raideur sont nul à ce niveau de la méthode }; break; case 7: case 8: // par rapport aux autres cas, on neutralise rien break; default : cout << "\n *** cas pas pris en compte !! \n SM_K_charge_contact() " << endl; Sortie(1); }; }; // énergie élastique de pénalisation: énergie élastique = 1/2 k x^2 = 1/2 F x energie_penalisation = 0.5 * F_N_max * gap; // la partie déplacement normal // on initialise la partie frottement avec au moins la partie pénalisation, ensuite // ce sera éventuellement abondé par le glissement ? ou la force de collage .... en devenir energie_frottement.ChangeEnergieElastique(energie_penalisation); // -- maintenant on regarde s'il faut s'occuper du frottement Coordonnee F_T(dima); Element * elemf = elfront->PtEI(); // récup de l'élément fini CompFrotAbstraite* loiFrot = ((ElemMeca*) elemf)->LoiDeFrottement(); penalisation_tangentielle = 0.; // init bool force_T_limiter = false; // init if (loiFrot != NULL) { // --- cas de l'existance d'une loi de frottement // calcul du déplacement tangentiel par rapport au premier point de contact // mis à jour en fonction du glissement éventuelle ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite const Coordonnee & M_C0 = elfro.Metrique()->PointM_tdt // ( stockage a tdt) (elfro.TabNoeud(),phi_theta_0); ////---- debug ---- //{const Coordonnee & M_a0 = elfro.Metrique()->PointM_0 // ( stockage a tdt) // (elfro.TabNoeud(),phi_theta_0); // cout << "\n debug smK contact: M_a0: ";M_a0.Affiche_1(cout); // cout << " a tdt: ";M_C0.Affiche_1(cout); //} ////---- fin debug --- Coordonnee M_0atdt = (Noe_atdt-M_C0); // déplacement total à partir du point d'impact mis à jour // transportée via l'interpolation: c-a-d les phi_theta_0 qui restent constant // ici du au frottement, le point d'impact est mis à jour en fonction du glissement éventuel dep_tangentiel = M_0atdt- N * (M_0atdt*N); // on retire le deplacement normal éventuel ////---- debug ---- //{ cout << "\n dep_tangentiel: ";dep_tangentiel.Affiche_1(cout); // cout << "\n Noe_atdt: "; Noe_atdt.Affiche_1(cout); //} ////---- fin debug --- double dep_T = dep_tangentiel.Norme(); dep_T_tdt = dep_T; // sauvegarde // // calcul du déplacement tangentiel: là il s'agit du déplacement relatif entre l'ancien point projeté et le nouveau // // Coordonnee M_tatdt(dima); // if (Mt.Dimension() != 0) // {M_tatdt = (Mtdt-Mt);} // déplacement total // ; // sinon Mt n'est pas encore définit donc on ne peut pas calculer le déplacement // // il est donc laissé à 0. // dep_tangentiel = M_tatdt- N * (M_tatdt*N); // on retire le deplacement normal éventuel // calcul de la force tangentielle en supposant tout d'abord un contact collant double penalisation_tangentielle; if (fct_nD_contact.fct_nD_penalisationTangentielle != NULL) {penalisation_tangentielle = Valeur_fct_nD(fct_nD_contact.fct_nD_penalisationTangentielle ,ElContact::Fct_nD_contact::tqi_fct_nD_penalisationTangentielle ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penalisationTangentielle ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penalisationTangentielle); } else {penalisation_tangentielle = ParaGlob::param->ParaAlgoControleActifs().PenalisationTangentielleContact();} // F_T =dep_tangentiel*(-ParaGlob::param->ParaAlgoControleActifs().PenalisationTangentielleContact()); F_T =dep_tangentiel*(-penalisation_tangentielle); // -- verification en appelant le comportement en frottement ---- // tout d'abord on récupère le delta t pour le calcul de la vitesse double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant(); double unSurDeltat=0; if (Abs(deltat) >= ConstMath::trespetit) {unSurDeltat = 1./deltat;} else // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand { // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb if (unSurDeltat < 0) { cout << "\n le pas de temps est négatif !! " << "\n ElContact::SM_charge_contact(..."; }; unSurDeltat = ConstMath::tresgrand; }; // calcul des vitesses Coordonnee vit_T = dep_tangentiel * unSurDeltat; // appel de la loi de comportement: calcul de la force de frottement et des énergies échangées Coordonnee nevez_force_frot;EnergieMeca energie; bool glisse = loiFrot->Cal_explicit_contact_tdt(vit_T,N,F_N,F_T,energie,deltat,nevez_force_frot); if (Permet_affichage() > 6) { cout << "\n deltaX_T: " << dep_tangentiel <<"\n vit_T= " << vit_T << " glisse= " << glisse << " penalisation_tangentielle= " << penalisation_tangentielle; }; // *** à faire // mise à jour du point d'impact en fct du glissement , ou d'un point de référence qui tient compte // du glissement // if (vit_T * nevez_force_frot > 0) // cout << "\n bizarre puissance de frottement positive !"; if (glisse) // on remplace la force de frottement si ça glisse sinon on conserve F_T F_T = nevez_force_frot; // --- calcul des énergies échangées energie_frottement += energie; } else if ( cas_collant) { // --- cas d'un contact collant sans loi de frottement explicite // calcul du déplacement tangentiel par rapport au premier point de contact ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite const Coordonnee & M_C0 = elfro.Metrique()->PointM_tdt // ( stockage a tdt) (elfro.TabNoeud(),phi_theta_0); ////---- debug ---- //{const Coordonnee & M_a0 = elfro.Metrique()->PointM_0 // ( stockage a tdt) // (elfro.TabNoeud(),phi_theta_0); // cout << "\n debug smK contact: M_a0: ";M_a0.Affiche_1(cout); // cout << " a tdt: ";M_C0.Affiche_1(cout); //} ////---- debug --- Coordonnee M_0atdt = (Noe_atdt-M_C0); // déplacement total à partir du point d'impact initial // transportée via l'interpolation: c-a-d les phi_theta_0 qui restent constant dep_tangentiel = M_0atdt- N * (M_0atdt*N); // on retire le deplacement normal éventuel ////---- debug ---- //{ cout << "\n dep_tangentiel: ";dep_tangentiel.Affiche_1(cout); // cout << "\n Noe_atdt: "; Noe_atdt.Affiche_1(cout); //} ////---- debug --- double dep_T = dep_tangentiel.Norme(); // calcul de la pénalisation double d_beta_dep_T=0.; // init de la dérivée du facteur de pénalisation / au déplacement // pour l'instant ne sert pas penalisation_tangentielle = CalFactPenalTangentiel(dep_T,d_beta_dep_T); dep_T_tdt = dep_T; // sauvegarde // calcul de la force tangentielle en supposant un contact collant F_T = -dep_tangentiel * penalisation_tangentielle; // limitation éventuelle de la force double max_force_T_noeud; if (fct_nD_contact.fct_nD_force_tangentielle_noeud_maxi != NULL) {max_force_T_noeud = Dabs(Valeur_fct_nD(fct_nD_contact.fct_nD_force_tangentielle_noeud_maxi ,ElContact::Fct_nD_contact::tqi_fct_nD_force_tangentielle_noeud_maxi ,ElContact::Fct_nD_contact::tqi_const_fct_nD_force_tangentielle_noeud_maxi ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_force_tangentielle_noeud_maxi)); } else {max_force_T_noeud = Dabs(ParaGlob::param->ParaAlgoControleActifs().Force_tangentielle_noeud_maxi());} double n_F_T = Dabs(penalisation_tangentielle * dep_T); ////---- debug ---- //{ cout << "\n n_F_T= " << n_F_T ; // cout << " F_T: ";F_T.Affiche_1(cout); //} ////---- debug --- if (Permet_affichage() > 6) { cout << "\n deltaX_T: " << dep_tangentiel << " penalisation_tangentielle= " << penalisation_tangentielle; }; if ((n_F_T > max_force_T_noeud) && (max_force_T_noeud > ConstMath::petit)) {F_T *= max_force_T_noeud/n_F_T; force_T_limiter = true; // on recalcul la pénalisation associée if (Dabs(dep_T) > ConstMath::trespetit) {penalisation_tangentielle = max_force_T_noeud / dep_T;} if (Permet_affichage() >= 6) {cout << " ** limitation penalisation tangentielle= " << penalisation_tangentielle << " intensite force = " << max_force_T_noeud << flush;}; }; ////---- debug ---- //{ cout << "\n F_T_max= " << n_F_T ; // cout << " F_T: ";F_T.Affiche_1(cout); //} ////---- debug --- // --- calcul des énergies échangées // ici il s'agit d'une énergie élastique, que l'on ajoute à la pénalisation normale // énergie élastique de pénalisation: énergie élastique = 1/2 k x^2 = 1/2 F x energie_penalisation -= 0.5 * F_T * dep_tangentiel; // - pour que le résultat soit > 0 energie_frottement.ChangeEnergieElastique(energie_penalisation); }; F_T_max = F_T.Norme(); // sauvegarde // --- calcul des puissances virtuelles force_contact = F_T + F_N; if (Permet_affichage() > 5) { cout << " noeud: " << noeud->Num_noeud() << ": force contact " << force_contact ; cout << " ||F_N||= "<< F_N.Norme() << " ||F_T||= "<< F_T_max; if (force_limiter) cout << " force normale en limite " ; if (force_T_limiter) cout << " force tangentielle en limite " ; }; if (Permet_affichage() > 6) cout << "\n ==> force contact imposee au residu (noeud:" << noeud->Num_noeud() << " maillage:"<< noeud->Num_Mail() << ") " << force_contact; int posi = Id_nom_ddl("R_X1") -1; // préparation pour les réactions Ddl_enum_etendu& ddl_reaction_normale=Ddl_enum_etendu::Tab_FN_FT()(1); Ddl_enum_etendu& ddl_reaction_tangentielle=Ddl_enum_etendu::Tab_FN_FT()(2); // une petite manip pour le cas axisymétrique: dans ce cas les ddl à considérer pour le résidu et la raideur // doivent être uniquement x et y donc 2 ddl au lieu de 3. (rien ne dépend de z) int nb_ddl_en_var = dima; // init bool axi = false; if (ParaGlob::AxiSymetrie()) { nb_ddl_en_var -= 1; axi = true;}; // récup du kronecker qui va bien TenseurHB& IdHB = *Id_dim_HB(ParaGlob::Dimension()); // a- cas du noeud projectile int ism=1; // init cas normale if ((cas_solide == 0 ) || (cas_solide == 1)) // cas où le noeud est libre {for (int ix=1;ix<=nb_ddl_en_var;ix++,ism++) // pour le noeud projectile la vitesse virtuelle associé est {SM_res(ism) += force_contact(ix); // directement la vitesse du noeud // F_N = force_contact // = - N * intens_force = - N * (gap * penalisation) // = - N * ( (N * deltaX) * penalisation ); // donc force_contact(ix) = - N(ix) * ( (N * deltaX) * penalisation ) // et d_force_contact(ix) / d ddl_iy = - N(ix) * ( (N * d_deltaX_d_ddl_iy ) * penalisation ) // + d(- N(ix))_d_ddl_iy * ( (N * deltaX) * penalisation ) // + - N(ix) * ( (d_N_d_ddl_iy * deltaX) * penalisation ) // + - N(ix) * ( (N * deltaX) * d_penalisation_d_ddl_iy ) // avec d_penalisation_d_ddl_iy = d_penalisation_d_gap * d_gap_d_ddl_iy // et sachant que d_gap_d_ddl_iy = N(iy) pour le noeud central // on abonde au niveau de la réaction au noeud noeud->Ajout_val_tdt(Enum_ddl(ix+posi),force_contact(ix)); // tout d'abord les termes strictement relatifs au noeud esclave int jsm=1; if ( cas_collant) {for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++) { (*raideur)(ism,jsm) += penalisation * N(ix)*N(iy) // partie classique //+ N(ix) * ( gap * d_beta_gap * N(iy) ) // partie intégrant la variation de la pénalisation + penalisation_tangentielle *(IdHB(ix,iy) - N(ix)*N(iy)) // partie collant ; }; } else // sans frottement {for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++) { (*raideur)(ism,jsm) += penalisation * N(ix)*N(iy) // partie classique //+ N(ix) * ( gap * d_beta_gap * N(iy) ) ; // partie intégrant la variation de la pénalisation }; }; //--- fin raideur : noeud esclave - noeud esclave // dans le cas de l'axi il faut décaler jsm, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus // --- ensuite les termes de couplage noeud esclave - noeuds de la facette if (cas_solide == 0) // cas du bi_déformable, sinon pas de terme de couplage {int iddl_facette = 1; if (d_T_pt != NULL) // car même si on l'a demandé, il peut être null s'il n'existe pas (ex cas 1D ) {Tableau & d_T = *d_T_pt; // variation du vecteur normale: pour simplifier if ( cas_collant) {if (actif==1) // on doit prendre en compte que deltax n'est pas normal à la facette {for (int jr=1;jr<=tab_taille-1;jr++) // boucle sur les noeuds de la facette { for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) // boucle sur les ddl (*raideur)(ism,jsm) += penalisation * ( - phii(jr) * N(ix)*N(iy) + gap * d_T(iddl_facette)(ix) + N(ix) * (d_T(iddl_facette) * deltaX) ) + penalisation_tangentielle * phi_theta_0(jr) *(-IdHB(ix,iy)+N(ix)*N(iy)) // partie collant ; // dans le cas de l'axi il faut décaler jsm, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus } } else // sinon deltax est normal à la facette et d_N * deltaX = 0 {for (int jr=1;jr<=tab_taille-1;jr++) // boucle sur les noeuds de la facette { for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) // boucle sur les ddl (*raideur)(ism,jsm) += penalisation * ( - phii(jr) * N(ix)*N(iy) + gap * d_T(iddl_facette)(ix) ) + penalisation_tangentielle * phi_theta_0(jr) *(-IdHB(ix,iy)+N(ix)*N(iy)) // partie collant ; // dans le cas de l'axi il faut décaler jsm, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus } }; } else // sinon pas collant avec (d_T_pt != NULL) {if (actif==1) // on doit prendre en compte que deltax n'est pas normal à la facette {for (int jr=1;jr<=tab_taille-1;jr++) // boucle sur les noeuds de la facette { for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) // boucle sur les ddl (*raideur)(ism,jsm) += penalisation * ( - phii(jr) * N(ix)*N(iy) + gap * d_T(iddl_facette)(ix) + N(ix) * (d_T(iddl_facette) * deltaX) ); // dans le cas de l'axi il faut décaler jsm, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus } } else // sinon deltax est normal à la facette et d_N * deltaX = 0 {for (int jr=1;jr<=tab_taille-1;jr++) // boucle sur les noeuds de la facette { for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) // boucle sur les ddl (*raideur)(ism,jsm) += penalisation * ( - phii(jr) * N(ix)*N(iy) + gap * d_T(iddl_facette)(ix) ); // dans le cas de l'axi il faut décaler jsm, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus } }; }; } else // sinon cas d_T_pt == NULL et (cas_solide == 0) {if ( cas_collant) {for (int jr=1;jr<=tab_taille-1;jr++) // boucle sur les noeuds de la facette { for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) // boucle sur les ddl (*raideur)(ism,jsm) += penalisation * ( - phii(jr) * N(ix)*N(iy)) + penalisation_tangentielle * phi_theta_0(jr) *(-IdHB(ix,iy)+N(ix)*N(iy)); // partie collant // dans le cas de l'axi il faut décaler jsm, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus }; } else // sinon cas non collant {for (int jr=1;jr<=tab_taille-1;jr++) // boucle sur les noeuds de la facette { for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) // boucle sur les ddl (*raideur)(ism,jsm) += penalisation * ( - phii(jr) * N(ix)*N(iy)); // dans le cas de l'axi il faut décaler jsm, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus }; }; // fin du if (cas_collant) }; // fin cas d_T_pt == NULL et (cas_solide == 0) }; // fin du if (cas_solide == 0) }; // fin de la boucle ix sur le noeud esclave } // fin de la première partie de: if ((cas_solide == 0 ) || (cas_solide == 1)) else // sinon c-a-d: cas_solide différent de 0 et 1 c-a-d le noeud est bloqué // on s'occupe quand même des réactions, ce qui nous permettra de les récupérer même sur un solide {for (int ix=1;ix<=dima;ix++) // pour le noeud projectile la vitesse virtuelle associé est // on abonde au niveau de la réaction au noeud noeud->Ajout_val_tdt(Enum_ddl(ix+posi),force_contact(ix)); }; // dans le cas de l'axi il faut décaler ism pour la suite, ce qui correspond au z // if (axi) ism++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus // on abonde au niveau de la réaction normale et tangentielle pour le noeud esclave noeud->ModifDdl_etendu(ddl_reaction_normale).Valeur() += -intens_force; double intensite_tangentielle = F_T.Norme(); noeud->ModifDdl_etendu(ddl_reaction_tangentielle).Valeur() += intensite_tangentielle; // b- cas de la surface cible if ((cas_solide == 0 ) || (cas_solide == 2)) // cas où la facette est libre // pour les noeud de la facette maître, sachant que l'on a pour la force sur le noeud esclave // F_N = - N * intens_force = - N * (gap * penalisation) = - N * ( (N * deltaX) * penalisation ) // = force_contact // donc sur les noeuds "s" de la facette on a : // F_N^s = phi_s * N * ( (N * deltaX) * penalisation ) // d'autre part on a gap = N * (X(t+dt)-M_impact) = N * (X(t+dt)- phi_r * X(tdt)^r) // avec X(tdt)^r les coordonnées du noeud r de la facette // d'où la forme finale pour le noeud "s" de la facette : // F_N^s = phi_s * N * ( (N * (X(t+dt)- phi_r * X(tdt)^r)) * penalisation ) // = - phi_s * force_contact // // maintenant au niveau des raideurs: // les raideurs relativements au noeud esclave -> dF_N^s/dX(t+dt)(iy) // dF_N^s/dX(t+dt)(iy) = phi_s * N * (N(iy)) * penalisation // // les raideurs relativements entre noeuds de la facette // dF_N^s/dX(t+dt)(iy) = phi_s * N * ( (N(iy) * (- phi_r ) * penalisation ) // + phi_s * d_N/dX(t+dt)(iy) * ( (N * (X(t+dt)- phi_r * X(tdt)^r)) * penalisation ) // + phi_s * N * (d_N/dX(t+dt)(iy) * (X(t+dt)- phi_r * X(tdt)^r)) * penalisation ) // donc force_contact(ix) = - N(ix) * ( (N * deltaX) * penalisation ) // et d_force_contact(ix) / d ddl_iy = - N(ix) * ( (N * d_deltaX_d_ddl_iy ) * penalisation ) // + d(- N(ix))_d_ddl_iy * ( (N * deltaX) * penalisation ) // + - N(ix) * ( (d_N_d_ddl_iy * deltaX) * penalisation ) // + - N(ix) * ( (N * deltaX) * d_penalisation_d_ddl_iy ) // avec d_penalisation_d_ddl_iy = d_penalisation_d_gap * d_gap_d_ddl_iy // et sachant que d_gap_d_ddl_iy = N(iy) pour le noeud central {for (int is=1;is<=tab_taille-1;is++) { Noeud* noe = tabNoeud(is+1); // pour simplifier if (Permet_affichage() >= 7) cout << "\n reaction facette imposee au residu (noeud:" << noe->Num_noeud() << " maillage:"<< noe->Num_Mail() << ") " ; tabForce_cont(is) = -force_contact*phii(is); for (int ix=1;ix<=nb_ddl_en_var;ix++,ism++) { double force_imp = -force_contact(ix)*phii(is); SM_res(ism) += force_imp; if (Permet_affichage() > 6) cout << force_imp << " "; // on abonde au niveau de la réaction au noeud noe->Ajout_val_tdt(Enum_ddl(ix+posi),force_imp); int jsm = 1; // tout d'abord le terme de couplage if (cas_solide == 0) // on intervient que si c'est déformable-déformable { if ( cas_collant) {for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++) { (*raideur)(ism,jsm) += - phii(is) * (penalisation * (N(ix)*N(iy)) + penalisation_tangentielle * (IdHB(ix,iy) - N(ix)*N(iy)) // partie collant ); } } else // cas non collant {for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++) { (*raideur)(ism,jsm) += - phii(is) * penalisation * N(ix)*N(iy);} } }; // dans le cas de l'axi il faut décaler jsm pour la suite, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus // puis les termes facette - facette int iddl_facette = 1; if (d_T_pt != NULL) // car même si on l'a demandé, il peut être null s'il n'existe pas (ex cas 1D ) {Tableau & d_T = *d_T_pt; // pour simplifier if (cas_collant) {if (actif==1) // on doit prendre en compte que deltax n'est pas normal à la facette {for (int js=1;js<=tab_taille-1;js++) {for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) (*raideur)(ism,jsm) += - phii(is) * ( penalisation *(N(ix) * ( - phii(js) * N(iy)) + d_T(iddl_facette)(ix) * gap + N(ix) * (d_T(iddl_facette) * deltaX) ) + penalisation_tangentielle * phi_theta_0(js) *(-IdHB(ix,iy)+N(ix)*N(iy)) // partie collant ); // dans le cas de l'axi il faut décaler jsm pour la suite, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus } } //fin cas actif == 1 else // sinon: actif diff de 1, deltax est normal à la facette et d_N * deltaX = 0 {for (int js=1;js<=tab_taille-1;js++) {for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) (*raideur)(ism,jsm) += - phii(is) * ( penalisation * (N(ix) * ( - phii(js) * N(iy)) + d_T(iddl_facette)(ix) * gap ) + penalisation_tangentielle * phi_theta_0(js) *(-IdHB(ix,iy)+N(ix)*N(iy)) // partie collant ); // dans le cas de l'axi il faut décaler jsm pour la suite, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus } }; // fin cas actif != 1 } else // sinon cas non collant {if (actif==1) // on doit prendre en compte que deltax n'est pas normal à la facette {for (int js=1;js<=tab_taille-1;js++) {for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) (*raideur)(ism,jsm) += - phii(is) * penalisation * ( N(ix) * ( - phii(js) * N(iy)) + d_T(iddl_facette)(ix) * gap + N(ix) * (d_T(iddl_facette) * deltaX) ); // dans le cas de l'axi il faut décaler jsm pour la suite, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus } } //fin cas actif == 1 else // sinon: actif diff de 1, deltax est normal à la facette et d_N * deltaX = 0 {for (int js=1;js<=tab_taille-1;js++) {for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) (*raideur)(ism,jsm) += - phii(is) * penalisation * ( N(ix) * ( - phii(js) * N(iy)) + d_T(iddl_facette)(ix) * gap ); // dans le cas de l'axi il faut décaler jsm pour la suite, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus } }; // fin cas actif != 1 }; // fin du if then else sur : if (cas_collant) } // fin du cas if (d_T_pt != NULL) else {if (cas_collant) {for (int js=1;js<=tab_taille-1;js++) {for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) (*raideur)(ism,jsm) += - phii(is) * ( penalisation * (N(ix) * ( - phii(js) * N(iy))) + penalisation_tangentielle * phi_theta_0(js) *(-IdHB(ix,iy)+N(ix)*N(iy)) // partie collant ); // dans le cas de l'axi il faut décaler jsm pour la suite, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus }; } else // sinon cas non collant {for (int js=1;js<=tab_taille-1;js++) {for (int iy=1;iy<=nb_ddl_en_var;iy++,jsm++,iddl_facette++) (*raideur)(ism,jsm) += - phii(is) * penalisation * ( N(ix) * ( - phii(js) * N(iy)) ); // dans le cas de l'axi il faut décaler jsm pour la suite, ce qui correspond au z // if (axi) jsm++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus }; }; }; // fin du cas où d_T_pt == NULL }; // dans le cas de l'axi il faut décaler ism pour la suite, ce qui correspond au z // if (axi) ism++; // qui n'est pas pris en compte // non: 20 mars 2019: maintenant le SM local a juste les ddl de X1 et X2 donc ce qui correspond à nb_ddl_en_var // du coup on ne décale plus // on abonde au niveau de la réaction normale et tangentielle pour les noeuds maîtres noe->ModifDdl_etendu(ddl_reaction_normale).Valeur() += -intens_force*phii(is); noe->ModifDdl_etendu(ddl_reaction_tangentielle).Valeur() += intensite_tangentielle*phii(is); }; } else // *** modif 23 juin 2015 , a virer si c'est ok // on s'occupe quand même des réactions, ce qui nous permettra de les récupérer même sur un solide {for (int ir=1;ir<=tab_taille-1;ir++) { Noeud* noe = tabNoeud(ir+1); // pour simplifier tabForce_cont(ir) = -force_contact*phii(ir); for (int ix=1;ix<=dima;ix++) { double force_imp = -force_contact(ix)*phii(ir); // on abonde au niveau de la réaction au noeud noe->Ajout_val_tdt(Enum_ddl(ix+posi),force_imp); }; // on abonde au niveau de la réaction normale et tangentielle pour les noeuds maîtres noe->ModifDdl_etendu(ddl_reaction_normale).Valeur() += -intens_force*phii(ir); noe->ModifDdl_etendu(ddl_reaction_tangentielle).Valeur() += intensite_tangentielle*phii(ir); }; // // --- debug //{ for (int ir=1;ir<=tab_taille-1;ir++) // { Noeud* noe = tabNoeud(ir+1); // pour simplifier // if ((noe->Num_noeud() == 5) && (noe->Num_Mail() == 2)) // {cout << "\n debug : ElContact::SM_K_charge_contact() "; // cout << " R_X= "; // for (int ix=1;ix<=dima;ix++) // cout << noe->Valeur_tdt(Enum_ddl(ix+posi)) << " "; // } // } //} // // --- fin debug }; // retour // return residu; // SM est nul à ce niveau du programme //// --- debug //cout << "\n debug : ElContact::SM_K_charge_contact() "; //SM_res.Affiche(); //raideur->Affiche(); //// --- fin debug return el; }; //================================= methodes privées =========================== // calcul la normale en fonction de differente conditions Coordonnee ElContact::Calcul_Normale(int dim, Plan & pl, Droite & dr,int indic) { Coordonnee N(dim); // le vecteur normal initialisé à 0 ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite // --- le cas des angles morts est particulier, on commence par le traiter if ((elfront->Angle_mort()) && (elfro.Type_geom_front() == POINT_G)) { // s'il s'agit d'un point il n'y a pas d'existance de normale // on utilise la normale au noeud const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t); const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); N = gr.ConteneurCoordonnee_const(); } // ---- maintenant on regarde 1) si on veut une normale lissée // 2) qui représente aussi le cas d'angle mort avec une ligne normalement qu'en 3D, // pour laquelle il n'y a pas d'existance d'une normale classique // on utilise ici la normale interpolée aux noeuds else if ( (normale_lisser) || ((elfront->Angle_mort()) && (elfro.Type_geom_front() == LIGNE)) ) {// dans le cas d'une normale lissée on va utiliser une interpolation des normales aux noeuds ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite // recup et actualisation du dernier plan tangent, coordonnées locale etc. // dimensionnement des variables de tangence // //debug // cout << "\n debug ElContact::Calcul_Normale: "; // cout << "\n elfro.DernierTangent " << flush; // // fin debug Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence elfro.DernierTangent(dr,pl,indic,false); // récup des fonctions phi const Vecteur& phii = elfro.Phi(); // //debug // cout << "\n debug ElContact::Calcul_Normale: "; // cout << "\n phii= " << phii; // // fin debug // on parcours les noeuds de la frontière // retourne le tableau de noeud en lecture lecture/écriture Tableau & tNfront = elfro.TabNoeud(); int nbNfront = tNfront.Taille(); Coordonnee Nnoe(dim); // variable inter for (int inoe = 1;inoe <= nbNfront;inoe++) { const TypeQuelconque& tiq = tNfront(inoe)->Grandeur_quelconque(N_FRONT_t); const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); Nnoe = gr.ConteneurCoordonnee_const(); ////debug // cout << "\n debug ElContact::Calcul_Normale: "; // cout << "\n Nnoe= " << Nnoe; //// fin debug //{TypeQuelconque& tiq_t = tNfront(inoe)->ModifGrandeur_quelconque(N_FRONT_t); // Grandeur_coordonnee& gr_t= *((Grandeur_coordonnee*) (tiq_t.Grandeur_pointee())); // Coordonnee& normale_t = *gr.ConteneurCoordonnee(); // //} N = phii(inoe)*Nnoe; // cout << "\n Nnoe= " << Nnoe; }; N.Normer(); } // --- on continue avec les éléments qui ne sont pas d'angle mort et sans lissage else if ( (dim == 3) && (indic == 2)) // cas 3d avec une surface { N = pl.Vecplan(); } else if ( ((dim == 2) && (indic == 1)) || (ParaGlob::AxiSymetrie() && (dim == 3) && (indic == 1)) ) // cas 2D avec une ligne, (ou 3D bloqué avec une ligne ??) // ou 3D axisymétrique avec une ligne // else if((dim == 2) && (indic == 1)) // // cas 2D ou 3D bloqué avec une ligne { Coordonnee V = dr.VecDroite(); // le vecteur tangent // ici il faut faire attention: la normale est prise de manière à être sortante de l'élément // car normalement, le vecteur unitaire de la droite est déterminée avec deux points suivants la numérotation de l'élément // donc la normale à ce vecteur, en suivant le sens trigo, se trouve dirigé vers l'intérieur de l'élément N(1) = V(2); N(2) = -V(1); } else if ( (dim == 3) && (indic == 1)) // cas 3d avec une ligne, on va considérer que la normale est la normale dans la direction // du noeud esclave (si celui-ci est suffisemment externe à la ligne { Coordonnee V = dr.VecDroite(); // On récupère le vecteur tangent à la ligne Coordonnee A = noeud->Coord2() - tabNoeud(2)->Coord2(); // on construit un vecteur entre le premier point de ligne -> le noeud esclave Coordonnee B = Util:: ProdVec_coor(V,A); // produit vectoriel // on ne continue que si ce produit vectoriel n'est pas nul if (B.Norme() > ConstMath::petit) { N = Util:: ProdVec_coor(V,B);} else // sinon cela veut dire que le noeud est sur la ligne, dans ce cas particulier, d'une manière arbitraire // on prend la normale dans un des 3 plans de base, en commençant par le plan xy qui serait celui où par exemple // on travaille en 3D, mais avec le z bloqué { // on commence par voir si c'est possible en xy if (DabsMaX(V(1),V(2)) > ConstMath::petit) // là c'est ok la ligne n'est pas // à z // ici il faut faire attention: la normale est prise de manière à être sortante de l'élément // car normalement, le vecteur unitaire de la droite est déterminée avec deux points suivants la numérotation de l'élément // donc la normale à ce vecteur, en suivant le sens trigo, se trouve dirigé vers l'intérieur de l'élément { N(1) = V(2); N(2) = -V(1); N(3)=0.; // on pose z=0 d'où une normale dans le plan xy } else // sinon cela veut dire que l'on est suivant l'axe de z, on choisit x par défaut { N(1)=1; N(2)=N(3)=0.;}; }; } else if ( (dim == 1) && (indic == 0)) // cas 1D avec un point { // la normale sortante de l'élément c'est soit 1 ou -1 // pour le déterminer on commence par trouver le coté de l'élément frontière qui va rentrer en contact // on est obligé de considérer l'élément // pour cela on cherche le noeud de l'élément qui correspond au noeud frontière // qui normalement est le second noeud du global des noeuds de l'élément de contact Noeud * no_front= tabNoeud(2); int num_no_front = no_front->Num_noeud(); const Tableau < Noeud * > t_N_elem = elfront->PtEI()->Tab_noeud_const(); int nbN_elem = t_N_elem.Taille(); int num_N=0; for (int i=1;i<= nbN_elem;i++) if (t_N_elem(i)->Num_noeud() == num_no_front) {num_N = i;break;}; if (num_N == 0) { cout << "\n *** warning on n'a pas trouve de frontiere correct pour le calcul de la normale " << "\n cas du contact entre deux points ... le contact sera mal calcule " << endl; }; // maintenant on prend au autre noeud de l'élément et on regarde le sens for (int i=1;i<= nbN_elem;i++) if (i != num_N) {// il s'agit bien d'un autre noeud if (t_N_elem(i)->Coord2()(1) > no_front->Coord2()(1) ) {N(1) = -1.;} else {N(1) = 1.;} break; }; } else // autres cas non traite pour l'instant { cout << "\n erreur, desole mais le cas de contact en dimension = " << dim << ", avec "; if (indic == 1) cout << " une droite "; else cout << " un plan "; cout << " n\'est pas actuellement traite \n" ; Affiche(1); cout << "ElContact::Normal(.." << endl; Sortie(1); }; return N; }; // récupération d'informations des classes internes // N: le vecteur normal // M_impact: le point d'impact sur la surface (ou ligne ou point) // phii : les fonctions d'interpolation au point d'impact // si avec_var est vrai: il y a retour du tableau de variation de la normale Tableau * ElContact::RecupInfo(Coordonnee& N,Coordonnee& M_impact,Vecteur& phii,bool avec_var) { ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite // recup de la dimension int dim = noeud->Dimension(); // recup du dernier plan tangent (ou droite) // dimensionnement des variables de tangence Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence Tableau * d_T = NULL; // init par défaut = elfro.DernierTangent(dr,pl,indic,avec_var); // different cas N.Change_dim(dim); // le vecteur normal // on commence par traiter le cas des éléments particuliers qui décrivent les angles morts if ((elfront->Angle_mort()) && (elfro.Type_geom_front() == POINT_G)) { // s'il s'agit d'un point il n'y a pas d'existance de normale // on utilise la normale au noeud elfro.DernierTangent(dr,pl,indic,false); // on ne considère pas la variation de la normale d_T = NULL; const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t); const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); N = gr.ConteneurCoordonnee_const(); M_impact = elfro.Ref(); } // ---- maintenant on regarde si on veut une normale lissée // ---- maintenant on regarde 1) si on veut une normale lissée // 2) qui représente aussi le cas d'angle mort avec une ligne normalement qu'en 3D, // pour laquelle il n'y a pas d'existance d'une normale classique // on utilise ici la normale interpolée aux noeuds else if ( (normale_lisser) || ((elfront->Angle_mort()) && (elfro.Type_geom_front() == LIGNE)) ) {// dans le cas d'une normale lissée on va utiliser une interpolation des normales aux noeuds ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite // recup et actualisation du dernier plan tangent, coordonnées locale etc. // dimensionnement des variables de tangence Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence elfro.DernierTangent(dr,pl,indic,false); // le vecteur normal est fixe ici pendant les itérations d_T = NULL; // récup des fonctions phi const Vecteur& phii = elfro.Phi(); // on parcours les noeuds de la frontière // retourne le tableau de noeud en lecture lecture/écriture Tableau & tNfront = elfro.TabNoeud(); int nbNfront = tNfront.Taille(); Coordonnee Nnoe(dim); // variable inter N.Zero(); for (int inoe = 1;inoe <= nbNfront;inoe++) { const TypeQuelconque& tiq = tNfront(inoe)->Grandeur_quelconque(N_FRONT_t); const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); Nnoe = gr.ConteneurCoordonnee_const(); N += phii(inoe)*Nnoe; }; N.Normer(); ////debug //if (Permet_affichage() > 4) // {cout << "\n debug ElContact::RecupInfo N= "<Coord2(); break; default: {cout << "\n *** erreur indic =" << indic <<" c-a-d diff de 2 ou 1 ou 0 " << " on ne peut pas continuer !"; this->Affiche(1); } break; }; } else { // cas d'élément de contact autres que ceux qui représentent les angles morts // et des normales lissées Tableau * d_T = elfro.DernierTangent(dr,pl,indic,avec_var); if (ParaGlob::AxiSymetrie() && (dim == 3)) N(3) = 0.; // init pour le cas axi if ( (dim == 3) && (indic == 2)) // cas 3d avec une surface { N = pl.Vecplan(); M_impact = pl.PointPlan(); } else if ( ((dim == 2) && (indic == 1)) || (ParaGlob::AxiSymetrie() && (dim == 3) && (indic == 1)) ) // cas 2D avec une ligne // ou 3D axisymétrique avec une ligne { Coordonnee V = dr.VecDroite(); // le vecteur tangent // N(1) = -V(2); N(2) = V(1); // un vecteur normal // !!!! pour l'instant il y a une erreur de numérotation: on utilise en fait une numérotation // exactement inverse de la doc !!! N(1) = V(2); N(2) = -V(1); // un vecteur normal // en axi N(3) est déjà initialisé à 0 M_impact = dr.PointDroite(); } else if ( (dim == 3) && (indic == 1)) // cas 3d avec une ligne, on va considérer que la normale est la normale dans la direction // **** a revoir car c'est problématique quand le noeud esclave est très proche de la ligne // du coup on peut avoir une normale qui oscille d'un coté à l'autre de la ligne ce qui fait que l'on va // avoir ensuite une force de contact ou de collage !!! aléatoire !!! pas bon du tout // sans doute trouver autre chose : peut-être que c'est ok quand le point est franchement hors de la ligne // pas pas autrement // du noeud esclave (si celui-ci est suffisemment externe à la ligne { Coordonnee V = dr.VecDroite(); // On récupère le vecteur tangent à la ligne Coordonnee A = noeud->Coord2() - tabNoeud(2)->Coord2(); // on construit un vecteur entre le premier point de ligne -> le noeud esclave Coordonnee B = Util:: ProdVec_coor(V,A); // produit vectoriel // on ne continue que si ce produit vectoriel n'est pas nul if (B.Norme() > ConstMath::petit) { N = Util:: ProdVec_coor(V,B).Normer();} else // sinon cela veut dire que le noeud est sur la ligne, dans ce cas particulier, d'une manière arbitraire // on prend la normale dans un des 3 plans de base, en commençant par le plan xy qui serait celui où par exemple // on travaille en 3D, mais avec le z bloqué { // on commence par voir si c'est possible en xy if (DabsMaX(V(1),V(2)) > ConstMath::petit) // là c'est ok la ligne n'est pas // à z // ici il faut faire attention: la normale est prise de manière à être sortante de l'élément // car normalement, le vecteur unitaire de la droite est déterminée avec deux points suivants la numérotation de l'élément // donc la normale à ce vecteur, en suivant le sens trigo, se trouve dirigé vers l'intérieur de l'élément { N(1) = V(2); N(2) = -V(1); N(3)=0.; // on pose z=0 d'où une normale dans le plan xy } else // sinon cela veut dire que l'on est suivant l'axe de z, on choisit x par défaut { N(1)=1; N(2)=N(3)=0.;}; }; M_impact = dr.PointDroite(); } // else if ( (dim == 2) && (indic == 1)) // // cas 2D avec une ligne // { Coordonnee V = dr.VecDroite(); // le vecteur tangent //// N(1) = -V(2); N(2) = V(1); // un vecteur normal //// !!!! pour l'instant il y a une erreur de numérotation: on utilise en fait une numérotation //// exactement inverse de la doc !!! // N(1) = V(2); N(2) = -V(1); // un vecteur normal // M_impact = dr.PointDroite(); // } else if ( (dim == 1) && (indic == 0)) // cas 1D avec un point { // la normale sortante de l'élément c'est soit 1 ou -1 // pour le déterminer on commence par trouver le coté de l'élément frontière qui va rentrer en contact // on est obligé de considérer l'élément // pour cela on cherche le noeud de l'élément qui correspond au noeud frontière // qui normalement est le second noeud du global des noeuds de l'élément de contact Noeud * no_front= tabNoeud(2); int num_no_front = no_front->Num_noeud(); const Tableau < Noeud * > t_N_elem = elfront->PtEI()->Tab_noeud_const(); int nbN_elem = t_N_elem.Taille(); int num_N=0; for (int i=1;i<= nbN_elem;i++) if (t_N_elem(i)->Num_noeud() == num_no_front) {num_N = i;break;}; if (num_N == 0) { cout << "\n *** warning on n'a pas trouve de frontiere correct pour le calcul de la normale " << "\n cas du contact entre deux points ... le contact sera mal calcule " << endl; }; // maintenant on prend au autre noeud de l'élément et on regarde le sens for (int i=1;i<= nbN_elem;i++) if (i != num_N) {// il s'agit bien d'un autre noeud if (t_N_elem(i)->Coord2()(1) > no_front->Coord2()(1) ) {N(1) = -1.;} else {N(1) = 1.;} break; }; // non erreur !! M_impact = N; // sera modifié ensuite, la valeur n'a pas d'importance // a priori le point d'impact c'est le noeud frontière M_impact = no_front->Coord2(); } else if ( (dim == 1) && (indic == 1)) // cas 1D avec une droite { N(1) = 1.; // le vecteur normal, qui correspond au seul vecteur disponible M_impact = dr.PointDroite();; } else // autres cas non traite pour l'instant { cout << "\n erreur, desole mais le cas de contact en dimension = " << dim << ", avec "; if (indic == 1) cout << " une droite "; else cout << " un plan "; cout << " n\'est pas actuellement traite \n" ; Affiche(1); cout << "ElContact::RecupInfo(.."<< endl; Sortie(1); }; }; // récup des fonctions phi phii = elfro.Phi(); // retour return d_T; }; // mise à jour de cas_solide et donc de ddlElement en fonction de l'activité des ddl void ElContact::Mise_a_jour_ddlelement_cas_solide_assemblage() {cas_solide = 0; // init par défaut: déformable-déformable // on regarde l'activité des noeuds // 1-- cas du noeud esclave int dim = ParaGlob::Dimension(); bool esclave_solide = true; switch (dim) {case 3: esclave_solide = noeud->Ddl_fixe(X3) && esclave_solide ; case 2: esclave_solide = noeud->Ddl_fixe(X2) && esclave_solide ; case 1: esclave_solide = noeud->Ddl_fixe(X1) && esclave_solide ; }; // 2-- la frontière int tal = tabNoeud.Taille(); bool front_solide = true; for (int i=2;i<=tal;i++) { Noeud* noe = tabNoeud(i); switch (dim) {case 3: front_solide = noe->Ddl_fixe(X3) && front_solide ; case 2: front_solide = noe->Ddl_fixe(X2) && front_solide ; case 1: front_solide = noe->Ddl_fixe(X1) && front_solide ; }; }; // maintenant on traite les différents cas int dim_ddlelement = 0; if (esclave_solide && front_solide) {cas_solide = 3;dim_ddlelement=0;} else if (esclave_solide) {cas_solide = 2;dim_ddlelement=tal-1;} else if (front_solide) {cas_solide = 1;dim_ddlelement=1;} else {cas_solide = 0;dim_ddlelement=tal;}; // maintenant on dimensionne ddlElement éventuellement if (ddlElement_assemblage->NbNoeud() != dim_ddlelement) // il faut redimensionner { // en fait tous les éléments du tableau sont identiques et il suffit qu'ils existent // on commence par parcourir la liste pour trouver un bon candidat list ::iterator ili,ilifin=list_Ddl_global.end(); bool trouve = false; for (ili=list_Ddl_global.begin();ili !=ilifin; ili++) if ((*ili).NbNoeud() == dim_ddlelement) // on a trouvé un candidat {ddlElement_assemblage = &(*ili); trouve = true;}; // dans le cas où on n'a pas trouvé, on en cré un if (!trouve) { int dima = ParaGlob::Dimension(); // dans le cas où on est en axisymétrie il faut diminuer de 1 if (ParaGlob::AxiSymetrie()) dima--; DdlElement tab_ddl(dim_ddlelement,dima); int posi = Id_nom_ddl("X1") -1; for (int i =1; i<= dima; i++) for (int j=1; j<= dim_ddlelement; j++) tab_ddl.Change_Enum(j,i,Enum_ddl(i+posi)); list_Ddl_global.push_front(tab_ddl); ddlElement_assemblage = &(*list_Ddl_global.begin()); }; // ensuite il faut s'occuper du tableau d'assemblage tabNoeud_pour_assemblage.Change_taille(0); switch (cas_solide) {case 0: tabNoeud_pour_assemblage = tabNoeud; // c'est le cas le plus complet break; case 1: tabNoeud_pour_assemblage.Change_taille(1); // seulement le noeud libre tabNoeud_pour_assemblage(1)=noeud; break; case 2: tabNoeud_pour_assemblage.Change_taille(tal-1); // la facette seule est libre for (int i=2;i<=tal;i++) tabNoeud_pour_assemblage(i-1) = tabNoeud(i); break; case 3: tabNoeud_pour_assemblage.Change_taille(0); // tout est bloqué break; }; }; }; // récup d'une place pour le résidu local et mise à jour de list_SM éventuellement void ElContact::RecupPlaceResidu(int nbddl) { // tout d'abord on vérifie que la place actuelle ne convient pas if (residu != NULL) {if (residu->Taille() == nbddl) {residu->Zero(); return; } // cas courant }; // dans tous les autres cas il faut soit récupérer une place qui convient soit en créé une autre // on commence par parcourir la liste pour trouver un bon candidat list ::iterator ilv,ilvfin=list_SM.end(); bool trouve = false; for (ilv=list_SM.begin();ilv !=ilvfin; ilv++) if ((*ilv).Taille() == nbddl) // on a trouvé un candidat {residu = &(*ilv); residu->Zero();trouve = true;}; // dans le cas où on n'a pas trouvé, on en cré un if (!trouve) { Vecteur interv(nbddl,0.); // vecteur mis à 0 par défaut list_SM.push_front(interv); residu = &(*list_SM.begin()); }; }; // récup d'une place pour la raideur locale et mise à jour de list_raideur éventuellement void ElContact::RecupPlaceRaideur(int nbddl) { // tout d'abord on vérifie que la place actuelle ne convient pas if (raideur != NULL) {if (raideur->Nb_ligne() == nbddl) {raideur->Initialise(0.); return; } // cas courant }; // dans tous les autres cas il faut soit récupérer une place qui convient soit en créé une autre // on commence par parcourir la liste pour trouver un bon candidat list ::iterator ilv,ilvfin=list_raideur.end(); bool trouve = false; for (ilv=list_raideur.begin();ilv !=ilvfin; ilv++) if ((*ilv).Nb_ligne() == nbddl) // on a trouvé un candidat {raideur = &(*ilv); raideur->Initialise(0.);trouve = true;}; // dans le cas où on n'a pas trouvé, on en cré un if (!trouve) { Mat_pleine interv(nbddl,nbddl,0.); // init à 0 par défaut list_raideur.push_front(interv); raideur = &(*list_raideur.begin()); }; }; // calcul du facteur de pénalisation en pénétration, en fonction de la géométrie // du module de compressibilité et des différents cas possibles // sauvegarde du gap // éventuellement, calcul de la dérivée: d_beta_gapdu, du facteur par rapport au gap // la sensibilité dépend du type de calcul du facteur de pénalisation double ElContact::CalFactPenal(const double& gap,double & d_beta_gap,int contact_type) { // récup du facteur donné par l'utilisateur double fact_penal=0.; // init par défaut if (fct_nD_contact.fct_nD_penalisationPenetration != NULL) {fact_penal = Valeur_fct_nD(fct_nD_contact.fct_nD_penalisationPenetration ,ElContact::Fct_nD_contact::tqi_fct_nD_penalisationPenetration ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penalisationPenetration ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penalisationPenetration); } else {fact_penal = ParaGlob::param->ParaAlgoControleActifs().PenalisationPenetrationContact();} double facteur=1.; // un facteur // choix en fonction du type de calcul int typ_cal= ParaGlob::param->ParaAlgoControleActifs().TypeCalculPenalisationPenetration(); // --- premier choix en fonction de typ_cal --- if (typ_cal != 0) {if (contact_type != 41) {switch (typ_cal) { case 1: // cas le plus simple : pas de calcul particulier break; // on ne fait rien case 2: case 3: case 4: case 5: case 6: case 7: case 8:// calcul du type ls-dyna avec la compressibilité du maître { // on choisit en fonction de la géométrie si il existe une compressibilité ElFrontiere* elfrontiere = elfront->Eleme(); // pour simplifier l'écriture const Element * elem = elfront->PtEI(); // idem if (((const ElemMeca*) elem)->CompressibiliteMoyenne() > 0.) { switch (elfrontiere->Type_geom_front()) { case SURFACE: if (elem->ElementGeometrie(X1).Dimension()==3) { // cas d'une surface frontière d'un volume double surf = elfrontiere->SurfaceApprox(); double volume = elem->Volume(); if (Dabs(volume) <= ConstMath::pasmalpetit) { if (Permet_affichage() > 1) cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; if (Permet_affichage() > 4) cout << "\n ElContact::CalFactPenal()"; facteur = elfrontiere->MaxDiagonale_tdt(); } else // sinon c'est bon {facteur = surf*surf/volume;} } else if (elem->ElementGeometrie(X1).Dimension()==2) { // cas d'une surface frontière d'un élément coque ou plaque double surf = elfrontiere->SurfaceApprox(); double maxdiag = elfrontiere->MaxDiagonale_tdt(); if (maxdiag <= ConstMath::pasmalpetit) { if (Permet_affichage() > 1) cout << "\n *** attention, la surface de l'element " << elem->Num_elt_const() << " du maillage " << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; if (Permet_affichage() > 4) cout << "\n ElContact::CalFactPenal()"; facteur = maxdiag; } else // sinon c'est bon {facteur = surf/maxdiag;} }; break; // pour mémoire pour la suite POINT_G = 1 , LIGNE , SURFACE, VOLUME, RIEN_TYPE_GEOM //******** partie en construction case LIGNE: // --- on traite les cas particuliers des éléments d'angle mort if (elfront->Angle_mort()) { // il n'y a pas de notion de surface de contact et l'idée est de récupérer // une grandeur du type module de compressibilité * une longueur caractéristique de l'élément // on retient une expression simple double longueur = elem->LongueurGeometrique_mini(1); facteur = longueur; //%%%% essai %%%% longueur = elem->LongueurGeometrique_mini(1); double volume = elem->Volume(); facteur = volume / (longueur * longueur); //%%%% fin essai %%%% } else {// --- cas d'un élément de contact autre que d'angle mort if (ParaGlob::AxiSymetrie()) // cas où l'on est en axisymétrique, { // on regarde la dimension de l'élément associé if (elem->ElementGeometrie(X1).Dimension()==2) // élément 2D en axi donc 3D en réalité { // récup de la longueur de la frontière double longueur = elfrontiere->LongueurApprox(); // on calcul la surface associée // M_impact(1) = x donc = le rayon car on est axi autour de y //// double surf = longueur * M_impact(1) * ConstMath::Pi * 2.; // 9 nov 2016: je crois qu'ici il y a une incompréhension. en fait il s'agit de la surface du maître // donc il ne faut pas tenir compte du point d'impact ?? sinon par exemple si le point est au centre // M_impact(1) = 0 et on n'a pas de surf !! // mais plus généralement, le pt d'impact ne doit pas interférer à mon avis double surf = longueur * ConstMath::Pi * 2.; // on récupère la surface de l'élément qui a créé la frontière double volume = elem->Volume(); if (Dabs(volume) <= ConstMath::pasmalpetit) { if (Permet_affichage() > 1) cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; if (Permet_affichage() > 4) cout << "\n ElContact::CalFactPenal()"; facteur = elfrontiere->MaxDiagonale_tdt(); } else // sinon c'est bon {facteur = surf*surf/volume;} } else // sinon ce n'est pas pris en charge {cout << "\n *** erreur: cas non pris en charge pour l'instant: " << " type de frontiere : " << Nom_type_geom(elfrontiere->Type_geom_front()) << " contact sur une ligne axisyme trique" << " \n ElContact::CalFactPenal() "; Sortie(1); }; } else if (elem->ElementGeometrie(X1).Dimension()==3) // cas d'une ligne frontière d'un élément coque ou plaque { // il faut finir la fonction ElFrontiere::LongueurApprox() // cf le laïus qui est écrit dans la méthode double surf = elfrontiere->SurfaceApprox(); double volume = elem->Volume(); // virtual double EpaisseurMoyenne(Enum_dure ) if (Dabs(volume) <= ConstMath::pasmalpetit) { if (Permet_affichage() > 1) cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; if (Permet_affichage() > 4) cout << "\n ElContact::CalFactPenal()"; facteur = elfrontiere->MaxDiagonale_tdt(); } else // sinon c'est bon {facteur = surf*surf/volume;} } else if (elem->ElementGeometrie(X1).Dimension()==2) // cas d'une ligne frontière d'une plaque { double longueur = elfrontiere->LongueurApprox(); double epais = 0.; // init if (!ParaGlob::AxiSymetrie()) // si on n'est pas en axi {epais = ((ElemMeca*) elem)->EpaisseurMoyenne(TEMPS_tdt );} else // si on est en axi {epais = 1.;}; double surf = elem->Volume() / epais; facteur = surf/longueur; }; }; break; //******** fin partie en construction case POINT_G: // --- on traite les cas particuliers des éléments d'angle mort if (elfront->Angle_mort()) { // il n'y a pas de notion de surface de contact et l'idée est de récupérer // une grandeur du type module de compressibilité * une longueur caractéristique de l'élément // on retient une expression simple double longueur = elem->LongueurGeometrique_mini(1); facteur = longueur; //%%%% essai %%%% longueur = elem->LongueurGeometrique_mini(1); double volume = elem->Volume(); facteur = volume / (longueur * longueur); //%%%% fin essai %%%% } else { // --- cas d'un élément de contact autre que d'angle mort // on ne considère que le cas d'une dimension 1D non axi, qui est un cas d'école // pour les autres cas, on considère que l'on ne doit pas prendre en compte des frontières // de type point, du coup on génèrera une erreur if (ParaGlob::Dimension() == 1 ) {// on veut calculer Ke * Ae * Ae / vol , comme on est dans le cas d'une barre // vol = long * section, et Ae = section , du coup: // Ke * Ae * Ae / vol = Ke * vol / (long * long) avec long = la longueur de l'élément double longueur = elem->LongueurGeometrique_mini(1); double volume = elem->Volume(); facteur = volume / (longueur * longueur); } else {cout << "\n *** erreur: cas non pris en charge pour l'instant: " << " type de frontiere : " << Nom_type_geom(elfrontiere->Type_geom_front()) << "\n contact en dimension " << ParaGlob::Dimension() << " \n ElContact::CalFactPenal() "; Sortie(1); }; } break; default: cout << "\n *** erreur: cas non pris en charge pour l'instant: " << " type de frontiere : " << Nom_type_geom(elfrontiere->Type_geom_front()) << " \n ElContact::CalFactPenal() "; Sortie(1); }; // on tiens compte maintenant de la compressibilité (moyenne) de l'élément facteur *= ((const ElemMeca*) elem)->CompressibiliteMoyenne(); }// fin du test d'existance d'une compressibilité else // on n'interdit pas car on pourrait avoir plusieurs maîtres, certains bien définit et // d'autres sans comportement donc sans compressibilité: {if (Permet_affichage() >5) cout << "\n contact: *** attention le module de compressibilite n'est (encore?) pas defini " << " on n'en tient pas compte (pour l'instant) pour le calcul du facteur de penalite normal ! "; }; break; } default: cout << "\n **** erreur cas de calcul du facteur de penalisation en penetration non existant " << " typ_cal= " << typ_cal << "\n ElContact::CalFactPenal()"; Sortie(1); }; // --- deuxième choix en fonction de typ_cal --- // on calcul la valeur final du facteur switch (typ_cal) { case 1: case 2: { fact_penal *= facteur;d_beta_gap = 0.; if (Permet_affichage() >= 5) cout << "\n fact_penal= " << fact_penal ; break; } case 3: // on fait varier le coef de pénalisation en fonction de l'enfoncement { if ((gap < gap_t) && (gap < 0.)) // la première fois gap_t = 0. donc on a forcément nb_pene_tdt >= 1 après le test { nb_pene_tdt++;} else if ((gap > gap_t ) && (gap < 0.)) // on diminue que si gap diminue, { if (nb_pene_tdt > 1) // on limite à 1 pour le calcul correcte de la puissance { nb_pene_tdt--; // cout << "\n ************* on diminue np_pene_tdt *************** "; }; }; // l'évolution du facteur prend en compte le nombre de fois où la pénétration augmente // 2**10 = 10^2, 2**100 = 10^30 !! const double a = 1.2; const double b=2.; fact_penal *= pow(a,nb_pene_tdt*b); //PUISSN(2.,nb_pene_tdt); // on modifie en conséquence le facteur d_beta_gap = 0.; // a priori on n'a pas vraiment de dérivée calculable ... if (Permet_affichage() >= 5) cout << "\n fact_penal= " << fact_penal << " nb_pene_tdt= " << nb_pene_tdt << " pow("<ParaAlgoControleActifs().Penetration_borne_regularisation();} double max_pene; if (fct_nD_contact.fct_nD_penetration_contact_maxi != NULL) {max_pene = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_contact_maxi); } else {max_pene = ParaGlob::param->ParaAlgoControleActifs().Penetration_contact_maxi();} // on calcul un facteur de majoration en fonction de la pénétration du pas précédent // double mult_pene = MaX(1., -gap_t/max_pene); double mult_pene = MaX(1., (-gap_t/max_pene)*mult_pene_t); // il faut appliquer au facteur précédent ! mult_pene_tdt = 0.5*(mult_pene_t + mult_pene); // on fait une moyenne glissante de 2 //////debug //if (mult_pene_tdt == 1) // {cout << "\n debug ElContact::CalFactPenal() temps " << ParaGlob::Variables_de_temps().TempsCourant(); // }; ////fin debug if (Permet_affichage() > 5) cout << "\n mult_pene_tdt= " << mult_pene_tdt ; // on calcul un facteur en fonction de la distance: thèse dominique chamoret, page 94 formule (3.3) if (gap <= (-borne_regularisation)) { fact_penal *= facteur*mult_pene_tdt;d_beta_gap = 0.;} else if (Dabs(gap) < borne_regularisation) // borne_regularisation est positif, donc on continue à avoir une pénalisation pour une pénétration positive !! // { fact_penal *= facteur * ((gap*0.5/borne_regularisation + 1.)*gap*0.5+borne_regularisation*0.25); // formule fausse de marmoret { fact_penal *= facteur *mult_pene_tdt * Sqr((gap-borne_regularisation)/borne_regularisation) * 0.25; // formule du même type d_beta_gap = facteur *mult_pene_tdt * (gap-borne_regularisation)/borne_regularisation/borne_regularisation * 0.5;; } else { fact_penal = 0;d_beta_gap = 0.;}; break; } case 5:case 7: // on fait varier le coef de pénalisation lorsque l'enfoncement est plus petit qu'un maxi { // ON récupe la pénétration maxi double borne_regularisation; if (fct_nD_contact.fct_nD_penetration_borne_regularisation != NULL) {borne_regularisation = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_borne_regularisation); } else {borne_regularisation = ParaGlob::param->ParaAlgoControleActifs().Penetration_borne_regularisation();} // on calcul un facteur en fonction de la distance: thèse dominique chamoret, page 94 formule (3.3) if (gap <= (-borne_regularisation)) {fact_penal *= facteur;d_beta_gap = 0.;} else if (Dabs(gap) < borne_regularisation) // borne_regularisation est positif, donc on continue à avoir une pénalisation pour une pénétration positive !! // { fact_penal *= facteur * ((gap*0.5/borne_regularisation + 1.)*gap*0.5+borne_regularisation*0.25); // formule fausse de marmoret { fact_penal *= facteur * Sqr((gap-borne_regularisation)/borne_regularisation) * 0.25; // formule du même type d_beta_gap = facteur * (gap-borne_regularisation)/borne_regularisation/borne_regularisation * 0.5;; } else { fact_penal = 0;d_beta_gap = 0.;}; break; } case 6: // idem cas 5 mais avec une fonction différente { // on récupère la pénétration maxi double borne_regularisation; if (fct_nD_contact.fct_nD_penetration_borne_regularisation != NULL) {borne_regularisation = Dabs(Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_borne_regularisation)); } else {borne_regularisation = Dabs(ParaGlob::param->ParaAlgoControleActifs().Penetration_borne_regularisation());} // on calcul un facteur en fonction de la distance: = (0.5*(tanh(4.*gap/borne)-1.)) // 0.25*exp(-gap/(0.25*borne)) if (gap <= (-borne_regularisation)) {fact_penal *= facteur;d_beta_gap = 0.;} else if (Dabs(gap) < borne_regularisation) { fact_penal *= -facteur * 0.5 * (tanh(gap*4./borne_regularisation)-1.); // formule du même type d_beta_gap = 4.*(1.-fact_penal*fact_penal)/borne_regularisation; // } else { fact_penal = 0;d_beta_gap = 0.;}; break; } case 8: // idem 4 mais pour un contact collant: quelque soit le sens de gap // on fait varier le coef de pénalisation lorsque l'enfoncement est plus petit qu'un maxi { // //------ essai à virer // fact_penal *= facteur;d_beta_gap = 0.; // if (Permet_affichage() >= 5) // cout << "\n fact_penal= " << fact_penal ; // break; // //------ fin essai // ON récupe la pénétration maxi double borne_regularisation; if (fct_nD_contact.fct_nD_penetration_borne_regularisation != NULL) {borne_regularisation = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_borne_regularisation); } else {borne_regularisation = ParaGlob::param->ParaAlgoControleActifs().Penetration_borne_regularisation();} double max_pene; if (fct_nD_contact.fct_nD_penetration_contact_maxi != NULL) {max_pene = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_contact_maxi); } else {max_pene = ParaGlob::param->ParaAlgoControleActifs().Penetration_contact_maxi();} // on calcul un facteur de majoration en fonction de la pénétration du pas précédent // double mult_pene = MaX(1., -gap_t/max_pene); double mult_pene = MaX(1., (Dabs(gap_t)/max_pene)*mult_pene_t); // il faut appliquer au facteur précédent ! mult_pene_tdt = 0.5*(mult_pene_t + mult_pene); // on fait une moyenne glissante de 2 //////debug //if (mult_pene_tdt == 1) // {cout << "\n debug ElContact::CalFactPenal() temps " << ParaGlob::Variables_de_temps().TempsCourant(); // }; ////fin debug if (Permet_affichage() > 5) cout << "\n mult_pene_tdt= " << mult_pene_tdt ; // on calcul un facteur en fonction de la distance: thèse dominique chamoret, page 94 formule (3.3) if (Dabs(gap) > borne_regularisation) { fact_penal *= facteur*mult_pene_tdt;d_beta_gap = 0.;} else // borne_regularisation est positif, donc on continue à avoir une pénalisation pour une pénétration positive !! // { fact_penal *= facteur * ((gap*0.5/borne_regularisation + 1.)*gap*0.5+borne_regularisation*0.25); // formule fausse de marmoret { fact_penal *= facteur *mult_pene_tdt * Sqr((-Dabs(gap)-borne_regularisation)/borne_regularisation) * 0.25; // formule du même type d_beta_gap = facteur *mult_pene_tdt * (-Dabs(gap)-borne_regularisation)/borne_regularisation/borne_regularisation * 0.5;; } break; } default: cout << "\n **** erreur 2 cas de calcul du facteur de penalisation en penetration non existant " << " typ_cal= " << typ_cal << "\n ElContact::CalFactPenal()"; Sortie(1); }; } else if (contact_type == 41)// si == 41, on utilise une mise à jour de la pénalisation { // la méthode est analogue au cas: typ_cal == 4 // ON récupe la pénétration maxi double borne_regularisation; if (fct_nD_contact.fct_nD_penetration_borne_regularisation != NULL) {borne_regularisation = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_borne_regularisation ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_borne_regularisation); } else {borne_regularisation = ParaGlob::param->ParaAlgoControleActifs().Penetration_borne_regularisation();} double max_pene; if (fct_nD_contact.fct_nD_penetration_contact_maxi != NULL) {max_pene = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_contact_maxi ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_contact_maxi); } else {max_pene = ParaGlob::param->ParaAlgoControleActifs().Penetration_contact_maxi();} // on calcul un facteur de majoration en fonction de la pénétration du pas précédent double mult_pene = MaX(1., (-gap_t/max_pene)); // il faut appliquer au facteur précédent ! if (Permet_affichage() > 5) cout << "\n mult_pene_tdt= " << mult_pene ; // on calcul un facteur en fonction de la distance: // if (gap <= (-borne_regularisation)) { fact_penal *= penalisation * mult_pene;d_beta_gap = 0.;} // else if (Dabs(gap) < borne_regularisation) // borne_regularisation est positif, donc on continue à avoir une pénalisation pour une pénétration positive !! // { fact_penal *= penalisation * mult_pene * Sqr((gap-borne_regularisation)/borne_regularisation) * 0.25; // formule du même type // d_beta_gap = penalisation * mult_pene * (gap-borne_regularisation)/borne_regularisation/borne_regularisation * 0.5;; // } // else // sinon on annule tout // { fact_penal = 0;d_beta_gap = 0.;}; } else if (contact_type == 42)// si == 41, on utilise une mise à jour de la pénalisation { }; }; // retour return fact_penal; }; // calcul du facteur de pénalisation en tangentiel, en fonction de la géométrie // du module de compressibilité et des différents possibles // sauvegarde du dep_T = déplacement tangentiel // éventuellement, calcul de la dérivée: d_beta_dep_Tdu, du facteur par rapport au dep_T // la sensibilité dépend du type de calcul du facteur de pénalisation // suis la même logique que pour la pénétration double ElContact::CalFactPenalTangentiel(const double& dep_T,double & d_beta_dep_T) { // récup du facteur donné par l'utilisateur double fact_penal = 0.; // init par défaut if (fct_nD_contact.fct_nD_penalisationTangentielle != NULL) {fact_penal = Valeur_fct_nD(fct_nD_contact.fct_nD_penalisationTangentielle ,ElContact::Fct_nD_contact::tqi_fct_nD_penalisationTangentielle ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penalisationTangentielle ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penalisationTangentielle); } else {fact_penal = ParaGlob::param->ParaAlgoControleActifs().PenalisationTangentielleContact();} double facteur=1.; // un facteur // choix en fonction du type de calcul int typ_cal= ParaGlob::param->ParaAlgoControleActifs().TypeCalculPenalisationTangentielle(); // --- premier choix en fonction de typ_cal --- if (typ_cal != 0) {switch (typ_cal) { case 1: // cas le plus simple : pas de calcul particulier break; // on ne fait rien case 2: case 3: case 4: case 5: case 6: case 7: case 8:// calcul du type ls-dyna avec la compressibilité du maître { // on choisit en fonction de la géométrie ElFrontiere* elfrontiere = elfront->Eleme(); // pour simplifier l'écriture const Element * elem = elfront->PtEI(); // idem if (((const ElemMeca*) elem)->CompressibiliteMoyenne() > 0.) { switch (elfrontiere->Type_geom_front()) { case SURFACE: if (elem->ElementGeometrie(X1).Dimension()==3) { // cas d'une surface frontière d'un volume double surf = elfrontiere->SurfaceApprox(); double volume = elem->Volume(); if (Dabs(volume) <= ConstMath::pasmalpetit) { if (Permet_affichage() > 1) cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; if (Permet_affichage() > 4) cout << "\n ElContact::CalFactPenalTangentiel()"; facteur = elfrontiere->MaxDiagonale_tdt(); } else // sinon c'est bon {facteur = surf*surf/volume;} } else if (elem->ElementGeometrie(X1).Dimension()==2) { // cas d'une surface frontière d'un élément coque ou plaque double surf = elfrontiere->SurfaceApprox(); double maxdiag = elfrontiere->MaxDiagonale_tdt(); if (maxdiag <= ConstMath::pasmalpetit) { if (Permet_affichage() > 1) cout << "\n *** attention, la surface de l'element " << elem->Num_elt_const() << " du maillage " << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; if (Permet_affichage() > 4) cout << "\n ElContact::CalFactPenalTangentiel()"; facteur = maxdiag; } else // sinon c'est bon {facteur = surf/maxdiag;} }; break; // pour mémoire pour la suite POINT_G = 1 , LIGNE , SURFACE, VOLUME, RIEN_TYPE_GEOM //******** partie en construction case LIGNE: // --- on traite les cas particuliers des éléments d'angle mort if (elfront->Angle_mort()) { // il n'y a pas de notion de surface de contact et l'idée est de récupérer // une grandeur du type module de compressibilité * une longueur caractéristique de l'élément // on retient une expression simple double longueur = elem->LongueurGeometrique_mini(1); facteur = longueur; //%%%% essai %%%% longueur = elem->LongueurGeometrique_mini(1); double volume = elem->Volume(); facteur = volume / (longueur * longueur); //%%%% fin essai %%%% } else {// --- cas d'un élément de contact autre que d'angle mort if (ParaGlob::AxiSymetrie()) // cas où l'on est en axisymétrique, { // on regarde la dimension de l'élément associé if (elem->ElementGeometrie(X1).Dimension()==2) // élément 2D en axi donc 3D en réalité { // récup de la longueur de la frontière double longueur = elfrontiere->LongueurApprox(); // on calcul la surface associée double surf = longueur * ConstMath::Pi * 2.; // on récupère la surface de l'élément qui a créé la frontière double volume = elem->Volume(); if (Dabs(volume) <= ConstMath::pasmalpetit) { if (Permet_affichage() > 1) cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; if (Permet_affichage() > 4) cout << "\n ElContact::CalFactPenalTangentiel()"; facteur = elfrontiere->MaxDiagonale_tdt(); } else // sinon c'est bon {facteur = surf*surf/volume;} } else // sinon ce n'est pas pris en charge {cout << "\n *** erreur: cas non pris en charge pour l'instant: " << " type de frontiere : " << Nom_type_geom(elfrontiere->Type_geom_front()) << " contact sur une ligne axisyme trique" << " \n ElContact::CalFactPenalTangentiel() "; Sortie(1); }; } else if (elem->ElementGeometrie(X1).Dimension()==3) // cas d'une ligne frontière d'un élément coque ou plaque { // il faut finir la fonction ElFrontiere::LongueurApprox() // cf le laïus qui est écrit dans la méthode double surf = elfrontiere->SurfaceApprox(); double volume = elem->Volume(); // virtual double EpaisseurMoyenne(Enum_dure ) if (Dabs(volume) <= ConstMath::pasmalpetit) { if (Permet_affichage() > 1) cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; if (Permet_affichage() > 4) cout << "\n ElContact::CalFactPenalTangentiel()"; facteur = elfrontiere->MaxDiagonale_tdt(); } else // sinon c'est bon {facteur = surf*surf/volume;} } else if (elem->ElementGeometrie(X1).Dimension()==2) // cas d'une ligne frontière d'une plaque { double longueur = elfrontiere->LongueurApprox(); double epais = 0.; // init if (!ParaGlob::AxiSymetrie()) // si on n'est pas en axi {epais = ((ElemMeca*) elem)->EpaisseurMoyenne(TEMPS_tdt );} else // si on est en axi {epais = 1.;}; double surf = elem->Volume() / epais; facteur = surf/longueur; }; } break; //******** fin partie en construction case POINT_G: // --- on traite les cas particuliers des éléments d'angle mort if (elfront->Angle_mort()) { // il n'y a pas de notion de surface de contact et l'idée est de récupérer // une grandeur du type module de compressibilité * une longueur caractéristique de l'élément // on retient une expression simple double longueur = elem->LongueurGeometrique_mini(1); facteur = longueur; } else { // --- cas d'un élément de contact autre que d'angle mort // on ne considère que le cas d'une dimension 1D non axi, qui est un cas d'école // pour les autres cas, on considère que l'on ne doit pas prendre en compte des frontières // de type point, du coup on génèrera une erreur if (ParaGlob::Dimension() == 1 ) {// on veut calculer Ke * Ae * Ae / vol , comme on est dans le cas d'une barre // vol = long * section, et Ae = section , du coup: // Ke * Ae * Ae / vol = Ke * vol / (long * long) avec long = la longueur de l'élément double longueur = elem->LongueurGeometrique_mini(1); double volume = elem->Volume(); facteur = volume / (longueur * longueur); } else {cout << "\n *** erreur: cas non pris en charge pour l'instant: " << " type de frontiere : " << Nom_type_geom(elfrontiere->Type_geom_front()) << "\n contact en dimension " << ParaGlob::Dimension() << " \n ElContact::CalFactPenalTangentiel() "; Sortie(1); }; }; break; default: cout << "\n *** erreur: cas non pris en charge pour l'instant: " << " type de frontiere : " << Nom_type_geom(elfrontiere->Type_geom_front()) << " \n ElContact::CalFactPenalTangentiel() "; Sortie(1); }; // on tiens compte maintenant de la compressibilité (moyenne) de l'élément facteur *= ((const ElemMeca*) elem)->CompressibiliteMoyenne(); }// fin du test d'existance d'une compressibilité else // on n'interdit pas car on pourrait avoir plusieurs maîtres, certains bien définit et // d'autres sans comportement donc sans compressibilité: {if (Permet_affichage() >5) cout << "\n contact: *** attention le module de compressibilite n'est (encore?) pas defini " << " on n'en tient pas compte (pour l'instant) pour le calcul du facteur de penalite tangentiel ! "; }; break; } default: cout << "\n **** erreur cas de calcul du facteur de penalisation tangentielle non existant " << " typ_cal= " << typ_cal << "\n ElContact::CalFactPenalTangentiel()"; Sortie(1); }; // --- deuxième choix en fonction de typ_cal --- // on calcul la valeur final du facteur switch (typ_cal) { case 1: case 2: { fact_penal *= facteur;d_beta_dep_T = 0.; if (Permet_affichage() >= 5) cout << "\n fact_penal= " << fact_penal ; break; } case 5:case 7: // on fait varier le coef de pénalisation lorsque le déplacement est plus petit qu'une borne { // de régularisation double borne_regularisation; if (fct_nD_contact.fct_nD_tangentielle_borne_regularisation != NULL) {borne_regularisation = Valeur_fct_nD(fct_nD_contact.fct_nD_tangentielle_borne_regularisation ,ElContact::Fct_nD_contact::tqi_fct_nD_tangentielle_borne_regularisation ,ElContact::Fct_nD_contact::tqi_const_fct_nD_tangentielle_borne_regularisation ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_tangentielle_borne_regularisation); } else {borne_regularisation = ParaGlob::param->ParaAlgoControleActifs().Tangentielle_borne_regularisation();} // on calcul un facteur en fonction de la distance if (Dabs(dep_T) > borne_regularisation) {fact_penal *= facteur;d_beta_dep_T = 0.;} else { fact_penal *= facteur * Sqr((dep_T-borne_regularisation)/borne_regularisation) * 0.25; d_beta_dep_T = facteur * (dep_T-borne_regularisation)/borne_regularisation/borne_regularisation * 0.5;; }; break; } case 6: // idem cas 5 mais avec une fonction différente { // on récupère le déplacement maxi double borne_regularisation; if (fct_nD_contact.fct_nD_tangentielle_borne_regularisation != NULL) {borne_regularisation = Dabs(Valeur_fct_nD(fct_nD_contact.fct_nD_tangentielle_borne_regularisation ,ElContact::Fct_nD_contact::tqi_fct_nD_tangentielle_borne_regularisation ,ElContact::Fct_nD_contact::tqi_const_fct_nD_tangentielle_borne_regularisation ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_tangentielle_borne_regularisation)); } else {borne_regularisation = Dabs(ParaGlob::param->ParaAlgoControleActifs().Tangentielle_borne_regularisation());} // on calcul un facteur en fonction de la distance: = (0.5*(tanh(4.*dep_T/borne)-1.)) // 0.25*exp(-dep_T/(0.25*borne)) if (Dabs(dep_T) > borne_regularisation) {fact_penal *= facteur;d_beta_dep_T = 0.;} else { fact_penal *= -facteur * 0.5 * (tanh(dep_T*4./borne_regularisation)-1.); // formule du même type d_beta_dep_T = 4.*(1.-fact_penal*fact_penal)/borne_regularisation; // }; break; } case 4: case 8: // quelque soit le sens de dep_T // on fait varier le coef de pénalisation lorsque le déplacement est plus petit qu'un maxi { ////------ essai à virer //fact_penal *= facteur;d_beta_dep_T = 0.; // if (Permet_affichage() >= 5) // cout << "\n fact_penal= " << fact_penal ; // break; ////------ fin essai double borne_regularisation; if (fct_nD_contact.fct_nD_tangentielle_borne_regularisation != NULL) {borne_regularisation = Valeur_fct_nD(fct_nD_contact.fct_nD_tangentielle_borne_regularisation ,ElContact::Fct_nD_contact::tqi_fct_nD_tangentielle_borne_regularisation ,ElContact::Fct_nD_contact::tqi_const_fct_nD_tangentielle_borne_regularisation ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_tangentielle_borne_regularisation); } else {borne_regularisation = ParaGlob::param->ParaAlgoControleActifs().Tangentielle_borne_regularisation();} double max_dep_T; if (fct_nD_contact.fct_nD_tangentielle_contact_maxi != NULL) {max_dep_T = Valeur_fct_nD(fct_nD_contact.fct_nD_tangentielle_contact_maxi ,ElContact::Fct_nD_contact::tqi_fct_nD_tangentielle_contact_maxi ,ElContact::Fct_nD_contact::tqi_const_fct_nD_tangentielle_contact_maxi ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_tangentielle_contact_maxi); } else {max_dep_T = ParaGlob::param->ParaAlgoControleActifs().Tangentielle_contact_maxi();} // on calcul un facteur de majoration en fonction du déplacement du pas précédent // double mult_tang = MaX(1., -dep_T_t/max_pene); double mult_tang = MaX(1., (Dabs(dep_T_t)/max_dep_T)*mult_tang_t); // il faut appliquer au facteur précédent ! mult_tang_tdt = 0.5*(mult_tang_t + mult_tang); // on fait une moyenne glissante de 2 //////debug //if (mult_tang_tdt == 1) // {cout << "\n debug ElContact::CalFactPenalTangentiel() temps " << ParaGlob::Variables_de_temps().TempsCourant(); // }; ////fin debug if (Permet_affichage() >= 5) cout << "\n mult_tang_tdt= " << mult_tang_tdt ; // on calcul un facteur en fonction de la distance: du même genre que pour la pénétration // sauf qu'ici c'est vrai quelque soit la direction du déplacement if (Dabs(dep_T) > borne_regularisation) { fact_penal *= facteur * mult_tang_tdt; d_beta_dep_T = 0.;} else { fact_penal *= facteur * mult_tang_tdt * Sqr((-Dabs(dep_T)-borne_regularisation)/borne_regularisation) * 0.25; d_beta_dep_T = facteur * mult_tang_tdt * (-Dabs(dep_T)-borne_regularisation)/borne_regularisation/borne_regularisation * 0.5;; }; break; } default: cout << "\n **** erreur 2 cas d'un calcul du facteur de penalisation tangentiel non existant " << " typ_cal= " << typ_cal << "\n ElContact::CalFactPenalTangentiel()"; Sortie(1); }; }; // dep_T_tdt = dep_T; // sauvegarde // retour return fact_penal; }; // calcul éventuel de la moyenne glissante des positions successive du noeud esclave // et changement des coordonnées du point passé en paramètre void ElContact::ChangeEnMoyGlissante(Coordonnee& Noe_atdt) { // pour simplifier int taille_moyenne_glissante = ParaGlob::param->ParaAlgoControleActifs().Nb_moy_glissant(); // on dimensionne éventuellement le tableau des positions successives if (tab_posi_esclave.Taille() != taille_moyenne_glissante) tab_posi_esclave.Change_taille(taille_moyenne_glissante,Coordonnee(ParaGlob::Dimension())); // on calcul éventuellement la moyenne glissante if (taille_moyenne_glissante > 1) { if (nb_posi_esclave_stocker_t < taille_moyenne_glissante) { // c'est le cas où la moyenne n'est pas finalisée nb_posi_esclave_stocker_tdt = nb_posi_esclave_stocker_t+1; indice_stockage_glissant_tdt=1; tab_posi_esclave(nb_posi_esclave_stocker_tdt) = noeud->Coord2(); // calcul de la moyenne Noe_atdt=tab_posi_esclave(1); // init for (int i=2;i<=nb_posi_esclave_stocker_tdt;i++) Noe_atdt += tab_posi_esclave(i); Noe_atdt /= nb_posi_esclave_stocker_tdt; } else // cas ou la moyenne est finalisée { tab_posi_esclave(indice_stockage_glissant_t) = noeud->Coord2(); indice_stockage_glissant_tdt=indice_stockage_glissant_t+1; // mise à jour de l'indice pour le prochain stockage // si l'indice dépasse la taille du taille, on remet au début if (indice_stockage_glissant_tdt > taille_moyenne_glissante) indice_stockage_glissant_tdt = 1; // calcul de la moyenne Noe_atdt=tab_posi_esclave(1); // init for (int i=2;i<=taille_moyenne_glissante;i++) Noe_atdt += tab_posi_esclave(i); Noe_atdt /= taille_moyenne_glissante; }; }; }; // calcul de la moyenne glissante de la pénalisation void ElContact::Moyenne_glissante_penalisation(double& penalisation, double& essai_penalisation) { // penalisation = la somme pondérée int tai_val_penal = val_penal.Taille(); // on vérifie que la taille n'a pas changée if (fct_pilotage_contact4 != NULL) {if (fct_pilotage_contact4->NbComposante() == 2) {Tableau & tava = fct_pilotage_contact4->Valeur_pour_variables_globales(); int nouvelle_taille = tava(2); if (nouvelle_taille != tai_val_penal) {if (nouvelle_taille > tai_val_penal) // on se contente d'augmenter la taille de stockage {val_penal.Change_taille(nouvelle_taille);} else // on va changer la taille en gardant les dernières valeurs (= les plus récentes) {int diff = tai_val_penal-nouvelle_taille; for (int i=1;i<= nouvelle_taille;i++) val_penal(i) = val_penal(i+diff); val_penal.Change_taille(nouvelle_taille); if (pt_dans_val_penal > nouvelle_taille) pt_dans_val_penal = 1; tai_val_penal = nouvelle_taille; }; if (Permet_affichage() > 2) cout << "\n -> Moyenne_glissante_penalisation: change taille " << tai_val_penal << " -> " << nouvelle_taille; }; }; }; if (val_penal(tai_val_penal) == 0.) // cas où le tableau n'est pas complètement rempli { penalisation *= (pt_dans_val_penal-1); // correspond à la somme penalisation = (penalisation+essai_penalisation)/pt_dans_val_penal; val_penal(pt_dans_val_penal)=essai_penalisation; // pour pt_dans_val_penal == 1 cela donne directement essai_penalisation } else // sinon le tableau est rempli, on effectue des remplacements { penalisation *= tai_val_penal; // la somme penalisation = (penalisation - val_penal(pt_dans_val_penal) + essai_penalisation )/tai_val_penal; val_penal(pt_dans_val_penal) = essai_penalisation; } pt_dans_val_penal++; // mise en place du retour périodique if (pt_dans_val_penal > tai_val_penal) pt_dans_val_penal = 1; }; // calcul d'une fonction nD relative à des données de contact double ElContact::Valeur_fct_nD(Fonction_nD * pt_fonct,Tableau < TypeQuelconque * >& tqi ,Tableau < const TypeQuelconque * >& tqii,Tableau & t_num_ordre ) const { // on commence par récupérer les conteneurs des grandeurs à fournir List_io & li_enu_scal = pt_fonct->Li_enu_etendu_scalaire(); List_io & li_quelc = pt_fonct->Li_equi_Quel_evolue(); const Tableau & tab_queconque = pt_fonct->Tab_enu_quelconque(); // on passe en revue les grandeurs pour les renseigner // --- on balaie maintenant la liste des grandeurs à sortir // 1) tout d'abord les ddl étendues Tableau val_ddl_enum (li_enu_scal.size()); List_io::const_iterator ie,iefin=li_enu_scal.end(); int it; // it est l'indice dans le tableau de retour for (it=1,ie=li_enu_scal.begin(); ie!=iefin;ie++,it++) { // si c'est un ddl pure, on regarde si elle est définie au noeud const Ddl_enum_etendu& enu = (*ie); // pour simplifier bool trouver = false; // on regarde d'abord s'il s'agit d'une info spécifique au contact int posi = (*ie).Position()-NbEnum_ddl(); switch (posi) { case 90: // reaction_normale -> en fait pas nécessaire, mais je laisse // car existe au noeud en contact {val_ddl_enum(it) = force_contact * N; trouver=true; break;} // case 91: // reaction_tangentielle déjà stockée au noeud // donc sera traité par le noeud en contact ///**** en fait les cas X1, X1_t et X1_t0 ne vont pas être traités ici /// mais vont être traités en grandeurs quelconques car il s'agit de vecteur // Je laisse quand même mais on ne passera sans doute jamais ici case 123: // "X1_t" { val_ddl_enum(it) = noeud->Coord1()(1);trouver=true; break;} case 124: // "X2_t" { val_ddl_enum(it) = noeud->Coord1()(2);trouver=true; break;} case 125: // X3_t { val_ddl_enum(it) = noeud->Coord1()(3);trouver=true; break;} case 126: // X1_t0 { val_ddl_enum(it) = noeud->Coord0()(1);trouver=true; break;} case 127: // X2_t0 { val_ddl_enum(it) = noeud->Coord0()(2);trouver=true; break;} case 128: // X3_t0 { val_ddl_enum(it) = noeud->Coord0()(3);trouver=true; break;} //*** fin default: // sinon on le stock en ddl_étendu trouver=false; break; };// fin du switch // si l'on n'a pas trouver, on regarde si l'info est dispo au noeud en contact if (!trouver) {if (enu.Position() <= NbEnum_ddl()) // cas d'un ddl pur {if (noeud->Existe_ici(enu.Enum())) {val_ddl_enum(it) = noeud->Ddl_noeud_tdt(enu.Enum()).Valeur(); trouver = true; }; } else // cas d'un ddl étendu { if (noeud->Existe_ici_ddlEtendu(enu)) {val_ddl_enum(it) = noeud->DdlEtendue(enu).ConstValeur(); trouver = true; }; }; }; // si on n'a rien trouvé if (!trouver) {cout << "\n erreur***: le ddl " << enu.Nom_plein() << " n'existe pas pour le noeud en contact " << " la fonction nD " << pt_fonct->NomFonction() << " ne peut pas etre renseignee " << flush; if (ParaGlob::NiveauImpression() > 2) {this->Affiche(); pt_fonct->Affiche(); } Sortie(1); }; };// -- fin de la boucle sur la liste de Ddl_enum_etendu // 2) puis les grandeurs quelconques {List_io::iterator ipq,ipqfin=li_quelc.end(); for (ipq=li_quelc.begin();ipq!=ipqfin;ipq++) { bool trouver = false; TypeQuelconque& enuTQ = (*ipq); // pour simplifier const TypeQuelconque_enum_etendu& en_ET_TQ = enuTQ.EnuTypeQuelconque().EnumTQ(); // on regarde tout d'abord les grandeurs spécifiques à l'élément contact switch (enuTQ.EnuTypeQuelconque().EnumTQ()) { // il y a des grandeurs vectorielles qui pour l'instant ne sont pas prises en compte // cf. fct nD // NORMALE_CONTACT, GLISSEMENT_CONTACT ,PENETRATION_CONTACT,FORCE_CONTACT, case CONTACT_NB_DECOL: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurEntier()) = nb_decol_tdt; trouver = true; break; } case CONTACT_PENALISATION_N: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurDouble()) = penalisation; trouver = true; break; } case CONTACT_PENALISATION_T: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurDouble()) = penalisation_tangentielle; trouver = true; break; } case CONTACT_NB_PENET: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurEntier()) = nb_pene_tdt; trouver = true; break; } case CONTACT_CAS_SOLIDE: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurEntier()) = cas_solide; trouver = true; break; } case CONTACT_ENERG_GLISSE_ELAS: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurDouble()) = energie_frottement.EnergieElastique(); trouver = true; break; } case CONTACT_ENERG_GLISSE_PLAS: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurDouble()) = energie_frottement.DissipationPlastique(); trouver = true; break; } case CONTACT_ENERG_GLISSE_VISQ: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurDouble()) = energie_frottement.DissipationVisqueuse(); trouver = true; break; } case CONTACT_ENERG_PENAL: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurDouble()) = energie_penalisation; trouver = true; break; } case NOEUD_PROJECTILE_EN_CONTACT: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); if (actif) {*(gr.ConteneurDouble()) = 100.;} else {*(gr.ConteneurDouble()) = -0.1; }; trouver = true; break; } case NOEUD_FACETTE_EN_CONTACT: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee())); int NBNF = tabNoeud.Taille(); if (actif) { for (int i=2;i<=NBNF;i++) *(gr.ConteneurDouble()) += 100.; } else { for (int i=2;i<=NBNF;i++) *(gr.ConteneurDouble()) -= 0.1; }; trouver = true; break; } case NUM_NOEUD: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurEntier()) = noeud->Num_noeud(); trouver = true; break; } case NUM_MAIL_NOEUD: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurEntier()) = noeud->Num_Mail(); trouver = true; break; } case POSITION_GEOMETRIQUE: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= noeud->Coord2(); trouver = true; break; } case POSITION_GEOMETRIQUE_t: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= noeud->Coord1(); trouver = true; break; } case POSITION_GEOMETRIQUE_t0: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= noeud->Coord0(); trouver = true; break; } // *** pour l'instant la suite n'est pas opérationelle car il s'agit de vecteur // dont les composantes n'ont pas d'équivalent en ddl_enum_etendu // il faudrait les définir si on veut pouvoir s'en servir case PENETRATION_CONTACT: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= Mtdt-noeud->Coord2(); trouver = true; break; } case PENETRATION_CONTACT_T: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= Mt-noeud->Coord1(); trouver = true; break; } case GLISSEMENT_CONTACT: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= dep_tangentiel; trouver = true; break; } case GLISSEMENT_CONTACT_T: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= dep_T_t; trouver = true; break; } case NORMALE_CONTACT: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= N; trouver = true; break; } case FORCE_CONTACT: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= force_contact; trouver = true; break; } case FORCE_CONTACT_T: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= force_contact_t; trouver = true; break; } case CONTACT_COLLANT: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurEntier()) = cas_collant; trouver = true; break; } case NUM_ZONE_CONTACT: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurEntier()) = num_zone_contact; trouver = true; break; } case NUM_ELEMENT: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurEntier()) = Elfront()->PtEI()->Num_elt(); trouver = true; break; } case NUM_MAIL_ELEM: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee())); *(gr.ConteneurEntier()) = Elfront()->PtEI()->Num_maillage(); trouver = true; break; } //*** fin default: trouver = false; break; }; // si l'on n'a pas trouver, on regarde si l'info est dispo au noeud en contact if (!trouver) {if (noeud->Existe_ici(en_ET_TQ)) { (*(enuTQ.Grandeur_pointee())) = (*noeud->Grandeur_quelconque(en_ET_TQ).Const_Grandeur_pointee()); trouver = true; }; }; // si on n'a rien trouvé if (!trouver) {cout << "\n erreur***: la grandeur " << en_ET_TQ.NomPlein() << " n'existe pas pour le noeud en contact " << " la fonction nD " << pt_fonct->NomFonction() << " ne peut pas etre renseignee " << flush; if (ParaGlob::NiveauImpression() > 2) {this->Affiche(); pt_fonct->Affiche(); } Sortie(1); }; }; }; // 3) et enfin les grandeurs quelconques: // les conteneurs sont normalement déjà bien dimensionné int tail_tab_queconque = tab_queconque.Taille(); // on va balayer les grandeurs quelconques for (int i=1;i<= tail_tab_queconque;i++) { bool trouver = false; EnumTypeQuelconque enu = tab_queconque(i); // on regarde tout d'abord les grandeurs spécifiques à l'élément contact switch (enu) { // il y a des grandeurs vectorielles qui pour l'instant ne sont pas prises en compte // cf. fct nD // NORMALE_CONTACT, GLISSEMENT_CONTACT ,PENETRATION_CONTACT,FORCE_CONTACT, case CONTACT_NB_DECOL: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurEntier()) = nb_decol_tdt; trouver = true; break; } case CONTACT_PENALISATION_N: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurDouble()) = penalisation; trouver = true; break; } case CONTACT_PENALISATION_T: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurDouble()) = penalisation_tangentielle; trouver = true; break; } case CONTACT_NB_PENET: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurEntier()) = nb_pene_tdt; trouver = true; break; } case CONTACT_CAS_SOLIDE: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurEntier()) = cas_solide; trouver = true; break; } case CONTACT_ENERG_GLISSE_ELAS: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurDouble()) = energie_frottement.EnergieElastique(); trouver = true; break; } case CONTACT_ENERG_GLISSE_PLAS: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurDouble()) = energie_frottement.DissipationPlastique(); trouver = true; break; } case CONTACT_ENERG_GLISSE_VISQ: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurDouble()) = energie_frottement.DissipationVisqueuse(); trouver = true; break; } case CONTACT_ENERG_PENAL: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurDouble()) = energie_penalisation; trouver = true; break; } case NOEUD_PROJECTILE_EN_CONTACT: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*tqi(i)).Grandeur_pointee())); if (actif) {*(gr.ConteneurDouble()) = 100.;} else {*(gr.ConteneurDouble()) = -0.1; }; trouver = true; break; } case NOEUD_FACETTE_EN_CONTACT: { Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*tqi(i)).Grandeur_pointee())); int NBNF = tabNoeud.Taille(); if (actif) { for (int i=2;i<=NBNF;i++) *(gr.ConteneurDouble()) += 100.; } else { for (int i=2;i<=NBNF;i++) *(gr.ConteneurDouble()) -= 0.1; }; trouver = true; break; } case NUM_NOEUD: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurEntier()) = noeud->Num_noeud(); trouver = true; break; } case NUM_MAIL_NOEUD: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurEntier()) = noeud->Num_Mail(); trouver = true; break; } case POSITION_GEOMETRIQUE: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*tqi(i)).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= noeud->Coord2(); trouver = true; break; } case POSITION_GEOMETRIQUE_t: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*tqi(i)).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= noeud->Coord1(); trouver = true; break; } case POSITION_GEOMETRIQUE_t0: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*tqi(i)).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= noeud->Coord0(); trouver = true; break; } // *** pour l'instant la suite n'est pas opérationelle car il s'agit de vecteur // dont les composantes n'ont pas d'équivalent en ddl_enum_etendu // il faudrait les définir si on veut pouvoir s'en servir case PENETRATION_CONTACT: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*tqi(i)).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= Mtdt-noeud->Coord2(); trouver = true; break; } case PENETRATION_CONTACT_T: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*tqi(i)).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= Mt-noeud->Coord1(); trouver = true; break; } case GLISSEMENT_CONTACT: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*tqi(i)).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= dep_tangentiel; trouver = true; break; } case GLISSEMENT_CONTACT_T: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*tqi(i)).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= dep_T_t; trouver = true; break; } case NORMALE_CONTACT: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*tqi(i)).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= N; trouver = true; break; } case FORCE_CONTACT: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*tqi(i)).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= force_contact; trouver = true; break; } case FORCE_CONTACT_T: {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*tqi(i)).Grandeur_pointee())); (*gr.ConteneurCoordonnee())= force_contact_t; trouver = true; break; } case CONTACT_COLLANT: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurEntier()) = cas_collant; trouver = true; break; } case NUM_ZONE_CONTACT: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurEntier()) = num_zone_contact; trouver = true; break; } case NUM_ELEMENT: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurEntier()) = Elfront()->PtEI()->Num_elt(); trouver = true; break; } case NUM_MAIL_ELEM: { Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*tqi(i)).Grandeur_pointee())); *(gr.ConteneurEntier()) = Elfront()->PtEI()->Num_maillage(); trouver = true; break; } //*** fin default: trouver = false; break; }; // si on n'a rien trouvé if (!trouver) {cout << "\n erreur***: la grandeur " << NomTypeQuelconque(enu) << " n'existe pas pour le noeud en contact " << " la fonction nD " << pt_fonct->NomFonction() << " ne peut pas etre renseignee " << flush; if (ParaGlob::NiveauImpression() > 2) {this->Affiche(); pt_fonct->Affiche(); } Sortie(1); }; }; // tqii pointe aux mêmes endroit que tqi mais est déclaré en const ... // calcul de la valeur et retour dans tab_ret Tableau & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,&(tqii),&t_num_ordre); #ifdef MISE_AU_POINT if (tab_val.Taille() != 1) { cout << "\nErreur : la fonction nD " << pt_fonct->NomFonction() << " doit calculer un scalaire or le tableau de retour est de taille " << tab_val.Taille() << " ce n'est pas normal !\n" << flush; if (Permet_affichage() > 2 ) cout << "\n ElContact::Valeur(... \n"; Sortie(1); }; #endif // on récupère le premier élément du tableau uniquement return tab_val(1); };