rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/herezh_pp/contact/ElContact_2.cc

2656 lines
138 KiB
C++
Executable file
Raw Permalink Blame History

// 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) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2021 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 <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
#include "ElContact.h"
#include "Droite.h"
#include "ConstMath.h"
#include "MathUtil.h"
#include "ElemMeca.h"
#include "Util.h"
#include <cmath>
#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 (niveau_commentaire > 5)
cout << "\n -- ElContact::SM_charge_contact: ";
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
// 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:"<<noeud->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()"<<endl;
//// fin debug
// limitation éventuelle au niveau de la pénétration maxi
// double max_pene = ParaGlob::param->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);}
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);}
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;
}
// 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
// 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));}
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 (niveau_commentaire > 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 (niveau_commentaire > 5) cout << "\n (penalisation limite) ";
// };
if (niveau_commentaire > 6) cout << " F_N(pur)= " << F_N << " ";
break;
}
default: break; // déjà vu au dessus
};
F_N_max = Dabs(intens_force); // sauvegarde
//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 (niveau_commentaire > 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 (niveau_commentaire >= 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 "<<force_limiter <<endl
// ;
// };
//// fin debug
// -- debug
//cout << ", F_N= " << (F_N * N) << " (ElContact::SM_charge_contact())" << endl;
// -- fin debug */
// énergie élastique de pénalisation: énergie élastique = 1/2 k x^2 = 1/2 F x
energie_penalisation = 0.5 * intens_force * 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
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);}
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;
}
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);}
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;
};
F_T_max = F_T.Norme();
// --- 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
force_contact = F_T + F_N;
if (niveau_commentaire >= 5)
{ cout << "\n noeud: " << noeud->Num_noeud() << ": force contact " << force_contact ;
if (force_limiter) cout << " force normale en limite " ;
if (force_T_limiter) cout << " force tangentielle en limite " ;
};
/* ////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 (niveau_commentaire >= 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 <<endl;
// fin debug */
return residu;
};
// -> implicite,
Element::ResRaid ElContact::SM_K_charge_contact()
{ if (niveau_commentaire > 5)
cout << "\n -- SM_K_charge_contact: ";
// 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()==56)
// {
// cout << "\n debug : ElContact::SM_K_charge_contact(): noeud 409"
// << endl;
// };
//// fin debug
Tableau <Coordonnee >* 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 ((niveau_commentaire > 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);}
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);}
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 (niveau_commentaire > 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 (niveau_commentaire > 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));}
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 (niveau_commentaire > 5) cout << "\n (penalisation limite) ";
// };
if (niveau_commentaire > 6) cout << " F_N(pur)= " << F_N << " ";
break;
}
default: break; // déjà vu au dessus
};
F_N_max = Dabs(intens_force); // sauvegarde
if (niveau_commentaire > 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 > 80.)
// {
// cout << "\n debug : ElContact::SM_K_charge_contact(): noeud " << noeud->Num_noeud()
// << " zone: " << num_zone_contact
// << "\n force " << F_N << " force_limiter "<<force_limiter <<endl
// <<"\n Noe_atdt= " << Noe_atdt << " M_impact= " << M_impact
// <<"\n gap= " << gap << " N= " << N
// << " penalisation= " << penalisation << endl
// ;
// };
// fin debug
// ici on a un traitement différent suivant les types de contact
// 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
if (niveau_commentaire > 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
if (niveau_commentaire >= 7)
{ 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: là il s'agit du déplacement relatif entre l'ancien point projeté et le nouveau
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);}
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 (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));}
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 ((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 (niveau_commentaire >= 6)
{cout << " ** limitation penalisation tangentielle= " << penalisation_tangentielle
<< " intensite force = " << max_force_T_noeud << flush;};
};
F_T_max = F_T.Norme(); // sauvegarde
////---- 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);
};
// --- calcul des puissances virtuelles
force_contact = F_T + F_N;
if (niveau_commentaire > 5)
{ cout << " noeud: " << noeud->Num_noeud() << ": force contact " << force_contact ;
if (force_limiter) cout << " force normale en limite " ;
if (force_T_limiter) cout << " force tangentielle en limite " ;
};
if (niveau_commentaire > 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 <Coordonnee >& 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 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 bloaqué
// 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 (niveau_commentaire >= 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 (niveau_commentaire > 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 <Coordonnee >& 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);
};
};
// 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
if ( (dim == 3) && (indic == 2))
// cas 3d avec une surface
{ N = pl.Vecplan();
return N;
}
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 "
<< "ElContact::Normal(int dim, Plan & pl, Droite & dr, indic)" << 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 <Coordonnee >* 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 <Coordonnee >* d_T = elfro.DernierTangent(dr,pl,indic,avec_var);
// different cas
N.Change_dim(dim); // le vecteur normal
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 " << "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 <DdlElement>::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 <Vecteur>::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 <Mat_pleine>::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 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);}
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 (niveau_commentaire > 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 (niveau_commentaire > 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 (niveau_commentaire > 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 (niveau_commentaire > 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:
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 (niveau_commentaire > 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 (niveau_commentaire > 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 (niveau_commentaire > 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 (niveau_commentaire > 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 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(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 (niveau_commentaire >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 (niveau_commentaire >= 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 (niveau_commentaire >= 5)
cout << "\n fact_penal= " << fact_penal << " nb_pene_tdt= "
<< nb_pene_tdt << " pow("<<a<<",nb_pene_tdt*"<<b<<") " <<pow(a,nb_pene_tdt*b);
break;
}
case 4: // 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);}
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);}
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 (niveau_commentaire > 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);}
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));}
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 (niveau_commentaire >= 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);}
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);}
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 (niveau_commentaire > 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);}
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);}
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 (niveau_commentaire > 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);}
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 (niveau_commentaire > 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 (niveau_commentaire > 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 (niveau_commentaire > 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 (niveau_commentaire > 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:
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 (niveau_commentaire > 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 (niveau_commentaire > 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 (niveau_commentaire > 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 (niveau_commentaire > 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 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(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 (niveau_commentaire >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 (niveau_commentaire >= 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);}
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));}
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 (niveau_commentaire >= 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);}
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);}
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 (niveau_commentaire >= 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 <double> & 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 (niveau_commentaire > 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) const
{ // on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
// 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 <double> val_ddl_enum (li_enu_scal.size());
List_io<Ddl_enum_etendu>::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 (niveau_commentaire > 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<TypeQuelconque>::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;
}
//*** 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 (niveau_commentaire > 2)
{this->Affiche();
pt_fonct->Affiche();
}
Sortie(1);
};
};
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
#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 (niveau_commentaire > 2 )
cout << "\n ElContact::Valeur(... \n";
Sortie(1);
};
#endif
// on récupère le premier élément du tableau uniquement
return tab_val(1);
};
Reference in a new issue View git blame Copy permalink