2021-09-27 12:42:13 +02:00
// 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.
//
2023-05-03 17:23:49 +02:00
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
2021-09-27 12:42:13 +02:00
// 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 "Debug.h"
# include "ElemMeca.h"
# include <iomanip>
# include "ConstMath.h"
# include "Util.h"
# include "Coordonnee2.h"
# include "Coordonnee3.h"
# include "TypeQuelconqueParticulier.h"
# include "TypeConsTens.h"
# include "Loi_iso_elas3D.h"
# include "Loi_iso_elas1D.h"
# include "Loi_iso_elas2D_C.h"
# include "FrontPointF.h"
# include "MathUtil2.h"
# include "Enum_StabMembrane.h"
// -------------------- stabilisation d'hourglass -------------------
// calcul d'élément de contrôle d'hourglass associée à un comportement
// str_precision : donne le type particulier d'élément à construire
void ElemMeca : : Init_hourglass_comp ( const ElemGeomC0 & elgeHour , const string & str_precision
, LoiAbstraiteGeneral * loiHourglass , const BlocGen & bloc )
{
//!!!! en fait pour l'instant on n'utilise pas elgehour, on considère que le nombre de pti complet est celui de l'élément sans info annexe !!!
// sauf pour l'hexaèdre quadratique pour lequel on a définit un str_precision particulier
// on met à jour les indicateurs
if ( Type_Enum_StabHourglass_existe ( bloc . Nom ( 1 ) ) )
{ type_stabHourglass = Id_Nom_StabHourglass ( bloc . Nom ( 1 ) . c_str ( ) ) ;
coefStabHourglass = bloc . Val ( 1 ) ;
}
else
{ cout < < " \n erreur, le type " < < bloc . Nom ( 2 ) < < " de gestion d'hourglass n'existe pas !! " ;
if ( ParaGlob : : NiveauImpression ( ) > 5 )
cout < < " \n ElemMeca::Init_hourglass_comp(... " ;
cout < < endl ;
Sortie ( 1 ) ;
} ;
switch ( type_stabHourglass )
{ case STABHOURGLASS_PAR_COMPORTEMENT :
{ // -- choix de l'element pour le calcul de la raideur -----
// on recherche un élément de même type,
// par défaut : numéro de maillage = le numéro de this, car c'est bien rattaché à ce maillage, mais l'élément est caché, donc non accessible par le maillage
// numéro d'élément = 1000000 + numéro de this : a priori c'est deux numéros n'ont pas d'importance, car ils ne sont pas
// rattachés à un maillage existant en tant que tel
int nUmMail = this - > Num_elt_const ( ) ;
int nUmElem = 1000000 + this - > Num_elt ( ) ;
tab_elHourglass . Change_taille ( 1 ) ;
tab_elHourglass ( 1 ) = ( ElemMeca * ) Element : : Choix_element ( nUmMail , nUmElem , Id_geometrie ( ) , Id_interpolation ( ) , Id_TypeProblem ( ) , str_precision ) ;
if ( tab_elHourglass ( 1 ) = = NULL )
{ // l'operation a echouee
cout < < " \n Erreur dans le choix d'un element interne utilise pour le blocage d'hourglass ****** " ;
cout < < " \n l \' element : " < < Nom_interpol ( id_interpol )
< < " " < < Nom_geom ( id_geom ) < < " " < < str_precision
< < " n \' est pas present dans le programme ! " < < endl ;
if ( ParaGlob : : NiveauImpression ( ) > 5 )
cout < < " \n ElemMeca::Init_hourglass_comp(... " < < endl ;
Sortie ( 1 ) ;
} ;
// --- définition des noeuds de l'élément -----
// on construit à partir des mêmes noeuds que ceux de l'élément père
// affectation des noeuds au nouvel élément
tab_elHourglass ( 1 ) - > Tab_noeud ( ) = this - > Tab_noeud ( ) ;
//--- def de la loi de comportement ---
// affectation de la loi
tab_elHourglass ( 1 ) - > DefLoi ( loiHourglass ) ;
break ;
}
case STABHOURGLASS_PAR_COMPORTEMENT_REDUIT :
{ // ici on a besoin de deux éléments. Le premier sert à faire le calcul complet
// le second sert à faire un calcul réduit, comme l'élément réel
tab_elHourglass . Change_taille ( 2 ) ;
// -- choix pour le calcul de la raideur -----
// on recherche un élément de même type,
// par défaut : numéro de maillage = le numéro de this, car c'est bien rattaché à ce maillage, mais l'élément est caché, donc non accessible par le maillage
// numéro d'élément = 1000000 + numéro de this : a priori c'est deux numéros n'ont pas d'importance, car ils ne sont pas
// rattachés à un maillage existant en tant que tel
//-- le premier élément:
int nUmMail = this - > Num_elt_const ( ) ;
int nUmElem = 1000000 + this - > Num_elt ( ) ;
tab_elHourglass ( 1 ) = ( ElemMeca * ) Element : : Choix_element ( nUmMail , nUmElem , Id_geometrie ( ) , Id_interpolation ( ) , Id_TypeProblem ( ) , str_precision ) ;
if ( tab_elHourglass ( 1 ) = = NULL )
{ // l'operation a echouee
cout < < " \n Erreur dans le choix d'un element interne utilise pour le blocage d'hourglass ****** " ;
cout < < " \n l \' element : " < < Nom_interpol ( id_interpol )
< < " " < < Nom_geom ( id_geom ) < < " " < < str_precision
< < " n \' est pas present dans le programme ! " < < endl ;
if ( ParaGlob : : NiveauImpression ( ) > 5 )
cout < < " \n ElemMeca::Init_hourglass_comp(... " < < endl ;
Sortie ( 1 ) ;
} ;
//-- le second élément identique à this: on utilise le même numéro de maillage
// même infos annexes
nUmElem = 20000000 + this - > Num_elt ( ) ; // on ajoute une unité pour le numéro
tab_elHourglass ( 2 ) = ( ElemMeca * ) Element : : Choix_element ( nUmMail , nUmElem , Id_geometrie ( ) , Id_interpolation ( ) , Id_TypeProblem ( ) , this - > infos_annexes ) ;
if ( tab_elHourglass ( 2 ) = = NULL )
{ // l'operation a echouee
cout < < " \n Erreur dans le choix d'un element interne utilise pour le blocage d'hourglass ****** " ;
cout < < " \n l \' element : " < < Nom_interpol ( id_interpol )
< < " " < < Nom_geom ( id_geom ) < < " " < < this - > infos_annexes
< < " n \' est pas present dans le programme ! " < < endl ;
if ( ParaGlob : : NiveauImpression ( ) > 5 )
cout < < " \n ElemMeca::Init_hourglass_comp(... " < < endl ;
Sortie ( 1 ) ;
} ;
// --- définition des noeuds de l'élément -----
// on construit à partir des mêmes noeuds que ceux de l'élément père
// affectation des noeuds au nouvel élément
tab_elHourglass ( 1 ) - > Tab_noeud ( ) = this - > Tab_noeud ( ) ;
tab_elHourglass ( 2 ) - > Tab_noeud ( ) = this - > Tab_noeud ( ) ;
//--- def de la loi de comportement ---
// affectation de la loi
tab_elHourglass ( 1 ) - > DefLoi ( loiHourglass ) ;
tab_elHourglass ( 2 ) - > DefLoi ( loiHourglass ) ;
// def de trois vecteurs de travail s'ils ne sont pas encore définit
if ( vec_hourglass . Taille ( ) < 3 )
vec_hourglass . Change_taille ( 3 ) ;
break ;
}
default :
cout < < " \n *** Erreur : cas de gestop, d'hourglass non defini ! \n " ;
cout < < " \n ElemMeca::Cal_implicit_hourglass() \n " ;
Sortie ( 1 ) ;
} ;
} ;
// stabilisation pour un calcul implicit
void ElemMeca : : Cal_implicit_hourglass ( )
{ E_Hourglass = 0. ; // init par défaut
switch ( type_stabHourglass )
{ case STABHOURGLASS_PAR_COMPORTEMENT :
{ // --- calcul de la raideur de l'élément, dans le cas implicite ---
Element : : ResRaid resraid = tab_elHourglass ( 1 ) - > Calcul_implicit ( ParaGlob : : param - > ParaAlgoControleActifs ( ) ) ;
// on tiend compte du facteur de modération
( * ( resraid . raid ) ) * = coefStabHourglass ;
( * ( resraid . res ) ) * = coefStabHourglass ;
// on met à jour la raideur et le résidu de l'élément principal
( * raideur ) + = ( * ( resraid . raid ) ) ;
( * residu ) + = ( * ( resraid . res ) ) ;
// on récupère les énergies stockées à l'élément
const EnergieMeca & energieTotale = tab_elHourglass ( 1 ) - > EnergieTotaleElement ( ) ;
E_Hourglass = coefStabHourglass * ( energieTotale . EnergieElastique ( )
+ energieTotale . DissipationPlastique ( ) + energieTotale . DissipationVisqueuse ( ) ) ;
//---- debug
//cout << "\n raideur locale après stabilisation d'hourglass ";
//matRaideur.Affiche();
//---- fin debug
} ;
break ;
case STABHOURGLASS_PAR_COMPORTEMENT_REDUIT :
{ // --- calcul de la raideur de l'élément, dans le cas implicite ---
// ici on calcul la partie supposée supplémentaire à celle déjà présente dans le calcul réduit
// au premier passage ou bien si on demande une grande précision de calcul, on calcul la raideur et le second membre,
// ensuite on ne calcul que le second membre
// if ((raid_hourglass_transitoire == NULL)||(prepa_niveau_precision > 8))
// en l'expérience prouve que c'est la forme des éléments du début qui est la moins tordu donc vaut mieux
// ne pas recalculer la raideur sur des éléments tordus, c'est dans ce cas que ça marche le mieux
if ( raid_hourglass_transitoire = = NULL )
{ // 1) calcul de la partie complète
Element : : ResRaid resraid1 = tab_elHourglass ( 1 ) - > Calcul_implicit ( ParaGlob : : param - > ParaAlgoControleActifs ( ) ) ;
// 2) calcul de la partie réduite
Element : : ResRaid resraid2 = tab_elHourglass ( 2 ) - > Calcul_implicit ( ParaGlob : : param - > ParaAlgoControleActifs ( ) ) ;
// on soustrait en gardant le résultat dans resraid1
// **** non ce n'est pas une bonne idée, lorsque l'on fait certain test de stabilité: lorsque le facteur est très grand
// si on ne fait pas la soustraction cela converge, avec la soustraction cela ne converge pas donc il ne faut pas soustraire
// (*(resraid1.raid)) -= (*(resraid2.raid));
// (*(resraid1.res)) -= (*(resraid2.res));
// ---nouvelle idée:
// on commence par récupérer dans le vecteur résidu les incréments de ddl correspondant à l'élément
int NBN = tab_noeud . Taille ( ) ;
int dim = ParaGlob : : Dimension ( ) ;
if ( ParaGlob : : AxiSymetrie ( ) )
dim - - ; // on passe à une dimension 2 ! et il y a uniquement les x et les y qui sont concernés
Vecteur & D_X = vec_hourglass ( 1 ) ; // pour simplier
Vecteur & X = vec_hourglass ( 2 ) ; // pour simplier
D_X . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
X . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
for ( int ine = 1 ; ine < = NBN ; ine + + )
{ Noeud & noe = * ( tab_noeud ( ine ) ) ; // pour simplifier
Coordonnee delta_X = noe . Coord2 ( ) - noe . Coord1 ( ) ;
Coordonnee Xnoeud = noe . Coord2 ( ) ;
int ipos = ( ine - 1 ) * dim ;
switch ( dim )
{ case 3 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; D_X ( ipos + 2 ) = delta_X ( 2 ) ; D_X ( ipos + 3 ) = delta_X ( 3 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; X ( ipos + 2 ) = Xnoeud ( 2 ) ; X ( ipos + 3 ) = Xnoeud ( 3 ) ; break ; }
case 2 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; D_X ( ipos + 2 ) = delta_X ( 2 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; X ( ipos + 2 ) = Xnoeud ( 2 ) ; break ; }
case 1 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; break ; }
} ;
} ;
// on multiplie la matrice par le vecteur delta_X
( * ( resraid1 . res ) ) = resraid1 . raid - > Prod_mat_vec ( D_X ) ;
// on sauvegarde la raideur initiale
if ( raid_hourglass_transitoire = = NULL )
raid_hourglass_transitoire = new Mat_pleine ( ( * ( resraid1 . raid ) ) ) ;
// premier calcul de l'accroissement de l'énergie
E_Hourglass = 0.5 * ( D_X * ( * ( resraid1 . res ) ) ) ;
// on tiend compte du facteur de modération
double coefMultiplicatif = coefStabHourglass ;
if ( E_Hourglass < 0. ) // on veut un accroissement d'énergie positive (1/2 k deltaX ^2)
{ coefMultiplicatif = - coefStabHourglass ; } ;
// grandeurs finales avec le bon signe
( * ( resraid1 . raid ) ) * = coefMultiplicatif ;
( * ( resraid1 . res ) ) * = coefMultiplicatif ;
// calcul définitif: correspond à une grandeur qui peut être comparée avec l'énergie élastique totale
E_Hourglass = 0.5 * ( X * ( * ( resraid1 . res ) ) ) ;
//debug
//cout << "\n premier passage, ";
//res.Affiche(); (resraid1.raid)->Affiche();
//fin debug
// on met à jour la raideur et le résidu de l'élément principal
( * raideur ) + = ( * ( resraid1 . raid ) ) ;
( * residu ) - = ( * ( resraid1 . res ) ) ;
}
else
{ // sinon on utilise la précédente raideur sauvegardée et on ne calcul que la partie résidue
// dans la nouvelle idée le vecteur doit exister à la bonne dimension, mais CalculResidu_tdt n'a pas besoin
// d'être appeler **** a améliorer
// Vecteur* resHourglass1 = NULL;
// resHourglass1 = tab_elHourglass(1)->CalculResidu_tdt (ParaGlob::param->ParaAlgoControleActifs());
/* Vecteur* resHourglass2 = NULL;
resHourglass2 = tab_elHourglass ( 2 ) - > CalculResidu_tdt ( ParaGlob : : param - > ParaAlgoControleActifs ( ) ) ;
// on soustrait les résidus en gardant le résultat dans resraid1
( * ( resHourglass1 ) ) - = ( * ( resHourglass2 ) ) ; */
// ---nouvelle idée:
// on commence par récupérer dans le vecteur résidu les incréments de ddl correspondant à l'élément
int NBN = tab_noeud . Taille ( ) ;
int dim = ParaGlob : : Dimension ( ) ;
if ( ParaGlob : : AxiSymetrie ( ) )
dim - - ; // on passe à une dimension 2 ! et il y a uniquement les x et les y qui sont concernés
Vecteur & D_X = vec_hourglass ( 1 ) ; // pour simplier
Vecteur & X = vec_hourglass ( 2 ) ; // pour simplier
D_X . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
X . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
for ( int ine = 1 ; ine < = NBN ; ine + + )
{ Noeud & noe = * ( tab_noeud ( ine ) ) ; // pour simplifier
Coordonnee delta_X = noe . Coord2 ( ) - noe . Coord1 ( ) ;
Coordonnee Xnoeud = noe . Coord2 ( ) ;
int ipos = ( ine - 1 ) * dim ;
switch ( dim )
{ case 3 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; D_X ( ipos + 2 ) = delta_X ( 2 ) ; D_X ( ipos + 3 ) = delta_X ( 3 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; X ( ipos + 2 ) = Xnoeud ( 2 ) ; X ( ipos + 3 ) = Xnoeud ( 3 ) ; break ; }
case 2 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; D_X ( ipos + 2 ) = delta_X ( 2 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; X ( ipos + 2 ) = Xnoeud ( 2 ) ; break ; }
case 1 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; break ; }
} ;
} ;
//debug
//cout << " res ddl, ";
//res.Affiche();
//fin debug
Vecteur & res = vec_hourglass ( 2 ) ; // pour simplier
res . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
// on multiplie le vecteur delta_X par la matrice
res = raid_hourglass_transitoire - > Prod_mat_vec ( D_X ) ;
// on définit un facteur de modération qui est fondé sur le ratio des énergies
// --- essai d'adaptation du coef, mais ça ne marche pas du tout, il n'y a que quand c'est fixe que c'est ok
/*
double energieHourglassDeBase = 0.5 * ( D_X * res ) ;
const EnergieMeca & energieActuelle = this - > EnergieTotaleElement ( ) ; // récup de l'énergie actuelle
double val_ref_energie = ( Dabs ( energieActuelle . EnergieElastique ( ) )
+ Dabs ( energieActuelle . DissipationPlastique ( ) )
+ Dabs ( energieActuelle . DissipationVisqueuse ( ) ) ) ;
double coef_a_appliquer = coefStabHourglass ;
if ( Dabs ( energieHourglassDeBase * coefStabHourglass ) > MaX ( ConstMath : : petit , val_ref_energie ) )
coef_a_appliquer * = val_ref_energie / energieHourglassDeBase ;
// on tiend compte du facteur de modération
res * = coef_a_appliquer ;
// on met à jour le résidu de l'élément principal
( * residu ) - = res ;
// pour la partie raideur: on met à jour la raideur de l'élément principal
( * raideur ) + = coef_a_appliquer * ( * raid_hourglass_transitoire ) ;
// calcul de l'énergie
E_Hourglass = energieHourglassDeBase * coef_a_appliquer ; //0.5 * (D_X * res);
*/
// premier calcul de l'accroissement de l'énergie
E_Hourglass = 0.5 * ( D_X * res ) ;
// on tiend compte du facteur de modération
double coefMultiplicatif = coefStabHourglass ;
if ( E_Hourglass < 0. ) // on veut un accroissement d'énergie positive (1/2 k deltaX ^2)
{ coefMultiplicatif = - coefStabHourglass ; } ;
// grandeurs finales avec le bon signe
res * = coefMultiplicatif ;
// calcul définitif: correspond à une grandeur qui peut être comparée avec l'énergie élastique totale
E_Hourglass = 0.5 * ( X * ( res ) ) ;
// on met à jour le résidu de l'élément principal
( * residu ) - = res ;
// pour la partie raideur: on met à jour la raideur de l'élément principal
( * raideur ) + = coefMultiplicatif * ( * raid_hourglass_transitoire ) ;
//debug
//if (E_Hourglass < 0.)
// { D_X.Affiche();
// res.Affiche();
// raid_hourglass_transitoire->Affiche();
// cout << "\n E_Hourglass= "<< E_Hourglass << endl ;
//
// };
//res.Affiche(); (raid_hourglass_transitoire)->Affiche();
//fin debug
} ;
//---- debug
//cout << "\n raideur locale après stabilisation d'hourglass ";
//matRaideur.Affiche();
//---- fin debug
} ;
break ;
case STABHOURGLASS_NON_DEFINIE :
// on ne fait rien
break ;
default :
cout < < " \n *** Erreur : cas de gestop, d'hourglass non defini ! \n " ;
cout < < " \n ElemMeca::Cal_implicit_hourglass() \n " ;
Sortie ( 1 ) ;
} ;
} ;
// stabilisation pour un calcul explicit
void ElemMeca : : Cal_explicit_hourglass ( bool atdt )
{ switch ( type_stabHourglass )
{ case STABHOURGLASS_PAR_COMPORTEMENT :
{ // --- calcul de la raideur de l'élément, dans le cas implicite ---
Vecteur * resHourglass = NULL ;
if ( atdt )
2023-05-03 17:23:49 +02:00
{ resHourglass = tab_elHourglass ( 1 ) - > CalculResidu_tdt ( ParaGlob : : param - > ParaAlgoControleActifs ( ) ) ; }
else
{ resHourglass = tab_elHourglass ( 1 ) - > CalculResidu_t ( ParaGlob : : param - > ParaAlgoControleActifs ( ) ) ; } ;
// on tiend compte du facteur de modération
( * resHourglass ) * = coefStabHourglass ;
// on met à jour le résidu de l'élément principal
( * residu ) + = ( * resHourglass ) ;
2021-09-27 12:42:13 +02:00
//---- debug
//cout << "\n raideur locale après stabilisation d'hourglass ";
//matRaideur.Affiche();
//---- fin debug
} ;
break ;
case STABHOURGLASS_PAR_COMPORTEMENT_REDUIT :
{ // --- on calcul la raideur de l'élément juste à la définition ---
if ( ( raid_hourglass_transitoire = = NULL ) | | ( prepa_niveau_precision > 18 ) ) // pour le debug
{ // 1) calcul de la partie complète
Element : : ResRaid resraid1 = tab_elHourglass ( 1 ) - > Calcul_implicit ( ParaGlob : : param - > ParaAlgoControleActifs ( ) ) ;
// on sauvegarde la raideur initiale
if ( raid_hourglass_transitoire = = NULL )
raid_hourglass_transitoire = new Mat_pleine ( ( * ( resraid1 . raid ) ) ) ;
// on commence par récupérer dans le vecteur résidu les incréments de ddl correspondant à l'élément
int NBN = tab_noeud . Taille ( ) ;
int dim = ParaGlob : : Dimension ( ) ;
if ( ParaGlob : : AxiSymetrie ( ) )
dim - - ; // on passe à une dimension 2 ! et il y a uniquement les x et les y qui sont concernés
Vecteur & D_X = vec_hourglass ( 1 ) ; // pour simplier
Vecteur & X = vec_hourglass ( 2 ) ; // pour simplier
D_X . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
X . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
for ( int ine = 1 ; ine < = NBN ; ine + + )
{ Noeud & noe = * ( tab_noeud ( ine ) ) ; // pour simplifier
Coordonnee delta_X = noe . Coord2 ( ) - noe . Coord1 ( ) ;
Coordonnee Xnoeud = noe . Coord2 ( ) ;
int ipos = ( ine - 1 ) * dim ;
switch ( dim )
{ case 3 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; D_X ( ipos + 2 ) = delta_X ( 2 ) ; D_X ( ipos + 3 ) = delta_X ( 3 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; X ( ipos + 2 ) = Xnoeud ( 2 ) ; X ( ipos + 3 ) = Xnoeud ( 3 ) ; break ; }
case 2 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; D_X ( ipos + 2 ) = delta_X ( 2 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; X ( ipos + 2 ) = Xnoeud ( 2 ) ; break ; }
case 1 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; break ; }
} ;
} ;
// on multiplie la matrice par le vecteur delta_X
( * ( resraid1 . res ) ) = resraid1 . raid - > Prod_mat_vec ( D_X ) ;
// premier calcul de l'accroissement de l'énergie
E_Hourglass = 0.5 * ( D_X * ( * ( resraid1 . res ) ) ) ;
// on tiend compte du facteur de modération
double coefMultiplicatif = coefStabHourglass ;
if ( E_Hourglass < 0. ) // on veut un accroissement d'énergie positive (1/2 k deltaX ^2)
{ coefMultiplicatif = - coefStabHourglass ; } ;
// grandeurs finales avec le bon signe
( * ( resraid1 . raid ) ) * = coefMultiplicatif ;
( * ( resraid1 . res ) ) * = coefMultiplicatif ;
// calcul définitif: correspond à une grandeur qui peut être comparée avec l'énergie élastique totale
E_Hourglass = 0.5 * ( X * ( * ( resraid1 . res ) ) ) ;
// on met à jour le résidu de l'élément principal
( * residu ) - = ( * ( resraid1 . res ) ) ;
// calcul de l'énergie
E_Hourglass = 0.5 * ( D_X * ( * ( resraid1 . res ) ) ) ;
}
else // cas courant
{ // on utilise la précédente raideur sauvegardée et on ne calcul qu'une partie résidue linéaire
// on commence par récupérer dans le vecteur résidu les incréments de ddl correspondant à l'élément
int NBN = tab_noeud . Taille ( ) ;
int dim = ParaGlob : : Dimension ( ) ;
if ( ParaGlob : : AxiSymetrie ( ) )
dim - - ; // on passe à une dimension 2 ! et il y a uniquement les x et les y qui sont concernés
Vecteur & D_X = vec_hourglass ( 1 ) ; // pour simplier
Vecteur & X = vec_hourglass ( 2 ) ; // pour simplier
D_X . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
X . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
for ( int ine = 1 ; ine < = NBN ; ine + + )
{ Noeud & noe = * ( tab_noeud ( ine ) ) ; // pour simplifier
Coordonnee delta_X = noe . Coord2 ( ) - noe . Coord1 ( ) ;
Coordonnee Xnoeud = noe . Coord2 ( ) ;
int ipos = ( ine - 1 ) * dim ;
switch ( dim )
{ case 3 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; D_X ( ipos + 2 ) = delta_X ( 2 ) ; D_X ( ipos + 3 ) = delta_X ( 3 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; X ( ipos + 2 ) = Xnoeud ( 2 ) ; X ( ipos + 3 ) = Xnoeud ( 3 ) ; break ; }
case 2 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; D_X ( ipos + 2 ) = delta_X ( 2 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; X ( ipos + 2 ) = Xnoeud ( 2 ) ; break ; }
case 1 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ;
X ( ipos + 1 ) = Xnoeud ( 1 ) ; break ; }
} ;
} ;
Vecteur & res = vec_hourglass ( 2 ) ; // pour simplier
res . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
// on multiplie le vecteur delta_X par la matrice
res = raid_hourglass_transitoire - > Prod_mat_vec ( D_X ) ;
// premier calcul de l'accroissement de l'énergie
E_Hourglass = 0.5 * ( D_X * res ) ;
// on tiend compte du facteur de modération
double coefMultiplicatif = coefStabHourglass ;
if ( E_Hourglass < 0. ) // on veut un accroissement d'énergie positive (1/2 k deltaX ^2)
{ coefMultiplicatif = - coefStabHourglass ; } ;
// grandeurs finales avec le bon signe
res * = coefMultiplicatif ;
// calcul définitif: correspond à une grandeur qui peut être comparée avec l'énergie élastique totale
E_Hourglass = 0.5 * ( X * ( res ) ) ;
// on met à jour le résidu de l'élément principal
( * residu ) - = res ;
} ;
} ;
break ;
case STABHOURGLASS_NON_DEFINIE :
// on ne fait rien
break ;
default :
cout < < " \n *** Erreur : cas de gestop, d'hourglass non defini ! \n " ;
cout < < " \n ElemMeca::Cal_implicit_hourglass() \n " ;
Sortie ( 1 ) ;
} ;
} ;
/*
// récupération de l'energie d'hourglass éventuelle
double ElemMeca : : Energie_Hourglass ( )
{ double enerHourglass = 0. ;
switch ( type_stabHourglass )
{ case STABHOURGLASS_PAR_COMPORTEMENT :
{ // on récupère les énergies stockées à l'élément
const EnergieMeca & energieTotale = tab_elHourglass ( 1 ) - > EnergieTotaleElement ( ) ;
enerHourglass = coefStabHourglass * ( energieTotale . EnergieElastique ( )
+ energieTotale . DissipationPlastique ( ) + energieTotale . DissipationVisqueuse ( ) ) ;
break ;
}
case STABHOURGLASS_PAR_COMPORTEMENT_REDUIT :
{ // --- l'énergie est purement une énergie élastique du type 0.5 <D_X> [K] (D_X) ---
if ( raid_hourglass_transitoire = = NULL )
{ enerHourglass = 0. ; } // car rien n'a été calculé
else
{ int NBN = tab_noeud . Taille ( ) ; // récup du nombre de noeuds
int dim = ParaGlob : : Dimension ( ) ;
Vecteur & D_X = vec_hourglass ( 1 ) ; // pour simplier
D_X . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
for ( int ine = 1 ; ine < = NBN ; ine + + )
{ Noeud & noe = * ( tab_noeud ( ine ) ) ; // pour simplifier
Coordonnee delta_X = noe . Coord2 ( ) - noe . Coord1 ( ) ;
int ipos = ( ine - 1 ) * dim ;
switch ( dim )
{ case 3 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; D_X ( ipos + 2 ) = delta_X ( 2 ) ; D_X ( ipos + 3 ) = delta_X ( 3 ) ; break ; }
case 2 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; D_X ( ipos + 2 ) = delta_X ( 2 ) ; break ; }
case 1 : { D_X ( ipos + 1 ) = delta_X ( 1 ) ; break ; }
} ;
} ;
Vecteur & res = vec_hourglass ( 2 ) ; // pour simplier
res . Change_taille ( NBN * dim ) ; // changement de la taille au cas où
// on multiplie la matrice par le vecteur delta_X
res = raid_hourglass_transitoire - > Prod_mat_vec ( D_X ) ;
enerHourglass = 0.5 * ( D_X * res ) ;
//debug
//cout << " , " << enerHourglass;
//fin debug
} ;
break ;
case STABHOURGLASS_NON_DEFINIE :
// on ne fait rien
break ;
}
} ;
return enerHourglass ;
} ; */
// -------------------- stabilisation transversale de membrane ou biel -------------------
// mise à jour de "a_calculer" en fonction du contexte
2023-05-03 17:23:49 +02:00
// NB: pour l'instant ne concerne que le calcul en implicite
2021-09-27 12:42:13 +02:00
void ElemMeca : : Mise_a_jour_A_calculer_force_stab ( )
{ if ( pt_StabMembBiel ! = NULL )
{ pt_StabMembBiel - > Aa_calculer ( ) = true ; // on init par défaut
Enum_StabMembraneBiel enu_stab = pt_StabMembBiel - > Type_stabMembrane ( ) ;
switch ( enu_stab )
{ case STABMEMBRANE_BIEL_PREMIER_ITER : case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER :
case STABMEMBRANE_BIEL_PREMIER_ITER_NORMALE_AU_NOEUD :
case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_NORMALE_AU_NOEUD :
{ // on récupère le pointeur correspondant aux iteration
const void * pointe = ( ParaGlob : : param - > GrandeurGlobal ( COMPTEUR_ITERATION_ALGO_GLOBAL ) ) ;
if ( pointe = = NULL )
{ cout < < " \n *** pb dans la stabilisation de membrane et/ou biel "
< < " le type de stabilisation est " < < Nom_StabMembraneBiel ( enu_stab )
< < " et la variable globale COMPTEUR_ITERATION_ALGO_GLOBAL n'est pas disponible "
< < " on ne peut pas continuer "
< < " \n ElemMeca::Mise_a_jour_A_calculer_force_stab(... " < < endl ;
Sortie ( 1 ) ;
} ;
// récup du conteneur du compteur
TypeQuelconque * gr_quelc = ( TypeQuelconque * ) ( pointe ) ;
Grandeur_scalaire_entier & gr
= * ( ( Grandeur_scalaire_entier * ) gr_quelc - > Grandeur_pointee ( ) ) ; // pour simplifier
int compteur_iter = * ( gr . ConteneurEntier ( ) ) ;
if ( compteur_iter > 1 ) // cela veut dire que la force de stabilisation a déjà été calculé
{ pt_StabMembBiel - > Aa_calculer ( ) = false ; }
else // sinon il faut la calculer
{ pt_StabMembBiel - > Aa_calculer ( ) = true ; } ;
break ;
}
case STABMEMBRANE_BIEL_PREMIER_ITER_INCR : case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR :
case STABMEMBRANE_BIEL_PREMIER_ITER_INCR_NORMALE_AU_NOEUD :
case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR_NORMALE_AU_NOEUD :
{ int compteur_iter = 0 ; // init
{ // on récupère le pointeur correspondant aux iteration
const void * pointe = ( ParaGlob : : param - > GrandeurGlobal ( COMPTEUR_ITERATION_ALGO_GLOBAL ) ) ;
if ( pointe = = NULL )
{ cout < < " \n *** pb dans la stabilisation de membrane et/ou biel "
< < " le type de stabilisation est " < < Nom_StabMembraneBiel ( enu_stab )
< < " et la variable globale COMPTEUR_ITERATION_ALGO_GLOBAL n'est pas disponible "
< < " on ne peut pas continuer "
< < " \n ElemMeca::Cal_implicit_StabMembBiel(... " < < endl ;
Sortie ( 1 ) ;
} ;
// récup du conteneur du compteur
TypeQuelconque * gr_quelc = ( TypeQuelconque * ) ( pointe ) ;
Grandeur_scalaire_entier & gr
= * ( ( Grandeur_scalaire_entier * ) gr_quelc - > Grandeur_pointee ( ) ) ; // pour simplifier
compteur_iter = * ( gr . ConteneurEntier ( ) ) ;
}
// on récupère le pointeur correspondant aux incréments
int compteur_incr = 0 ; // init
{ const void * pointe = ( ParaGlob : : param - > GrandeurGlobal ( COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL ) ) ;
if ( pointe = = NULL )
{ cout < < " \n *** pb dans la stabilisation de membrane et/ou biel "
< < " le type de stabilisation est STABMEMBRANE_BIEL_PREMIER_ITER_INCR "
< < " et la variable globale COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL n'est pas disponible "
< < " on ne peut pas continuer "
< < " \n ElemMeca::Cal_implicit_StabMembBiel(... " < < endl ;
Sortie ( 1 ) ;
} ;
// récup du conteneur du compteur
TypeQuelconque * gr_quelc = ( TypeQuelconque * ) ( pointe ) ;
Grandeur_scalaire_entier & gr
= * ( ( Grandeur_scalaire_entier * ) gr_quelc - > Grandeur_pointee ( ) ) ; // pour simplifier
compteur_iter = * ( gr . ConteneurEntier ( ) ) ;
}
if ( ( compteur_incr > 1 ) & & ( compteur_iter > 1 ) ) // cela veut dire que la force de stabilisation a déjà été calculé
{ pt_StabMembBiel - > Aa_calculer ( ) = false ; }
else // sinon il faut la calculer
{ pt_StabMembBiel - > Aa_calculer ( ) = true ; } ;
break ;
}
default :
cout < < " \n *** erreur le type de stabilisation de membrane et/ou biel n'est pas definie "
< < " on ne peut pas continuer "
< < " \n ElemMeca::Mise_a_jour_A_calculer_force_stab(... " < < endl ;
Sortie ( 1 ) ;
break ;
} ;
} ;
} ;
// stabilisation pour un calcul implicit,
// iteg -> donne le numéro du dernier pti sur lequel on a travaillé, auquel met est associé
// iteg == 0 : un seul calcul global, et on est à la suite d'une boucle
// iteg == -1 : fin d'un calcul avec boucle sur tous les pti, et on est à la suite de la boucle
// correspond au calcul d'alpha: == stab_ref
// iteg entre 1 et nbint: on est dans une boucle de pti
// nbint : le nombre maxi de pti
// poids_volume : si iteg > 0 : le poids d'intégration
// si iteg <= 0 : l'intégrale de l'élément (ici la surface totale)
void ElemMeca : : Cal_implicit_StabMembBiel ( int iteg , const Met_abstraite : : Impli & ex , const int & nbint
, const double & poids_volume
, Tableau < int > * noeud_a_prendre_en_compte )
{ // l'idée est de bloquer la direction normale à l'élément pour une plaque et
// les deux directions transverses pour une biel
2023-05-03 17:23:49 +02:00
int niveau_affichage = ParaGlob : : NiveauImpression ( ) ;
if ( pt_StabMembBiel - > Affichage_stabilisation ( ) > 0 )
niveau_affichage = pt_StabMembBiel - > Affichage_stabilisation ( ) ;
2021-09-27 12:42:13 +02:00
//---- détermination de l'intensité de stabilisation
// celle-ci n'est calculé que dans le cas ou iteg == 0
// c'est à dire que l'on n'est pas dans une boucle sur les pti, dans la méthode appelant
double alpha = 0. ; // init
bool stab_utilisable = false ; // init
// bool alpha_vient_etre_calculee = false;
Enum_StabMembraneBiel enu_stab = pt_StabMembBiel - > Type_stabMembrane ( ) ;
switch ( enu_stab )
{ case STABMEMBRANE_BIEL_PREMIER_ITER : case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER :
case STABMEMBRANE_BIEL_PREMIER_ITER_NORMALE_AU_NOEUD :
case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_NORMALE_AU_NOEUD :
{ // on récupère le pointeur correspondant aux iteration
const void * pointe = ( ParaGlob : : param - > GrandeurGlobal ( COMPTEUR_ITERATION_ALGO_GLOBAL ) ) ;
if ( pointe = = NULL )
{ cout < < " \n *** pb dans la stabilisation de membrane et/ou biel "
< < " le type de stabilisation est " < < Nom_StabMembraneBiel ( enu_stab )
< < " et la variable globale COMPTEUR_ITERATION_ALGO_GLOBAL n'est pas disponible "
< < " on ne peut pas continuer "
< < " \n ElemMeca::Cal_implicit_StabMembBiel(... " < < endl ;
Sortie ( 1 ) ;
} ;
// récup du conteneur du compteur
TypeQuelconque * gr_quelc = ( TypeQuelconque * ) ( pointe ) ;
Grandeur_scalaire_entier & gr
= * ( ( Grandeur_scalaire_entier * ) gr_quelc - > Grandeur_pointee ( ) ) ; // pour simplifier
int compteur_iter = * ( gr . ConteneurEntier ( ) ) ;
if ( ( iteg > 0 ) & & ( compteur_iter > 1 ) ) // cela veut dire que la la référence de stabilisation a déjà été calculée
{ alpha = pt_StabMembBiel - > Stab_ref ( ) ; // récup de l'intensité de stabilisation
// ce n'est pas un pointeur, alpha est une variable locale
// dans le cas d'une fonction multiplicatrice ,
// et il faut appeler la fonction qui peut varier pendant les itérations
if ( pt_StabMembBiel - > Pt_fct_gamma ( ) ! = NULL )
// cas où il y a une fonction nD
{ Fonction_nD * fct = pt_StabMembBiel - > Pt_fct_gamma ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = - 1 ; // signifie que l'on utilise la métrique qui vient juste d'être calculée
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi ( absolue , TEMPS_tdt , li_enu_scal , iteg , - cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , iteg , - cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// seule la première valeur est ok
alpha * = tab_ret ( 1 ) ;
pt_StabMembBiel - > Valgamma ( ) = alpha ; // on sauvegarde la valeur d'alpha
}
else // sinon cela veut dire que l'intensité de stabilisation est fixe
{ alpha * = pt_StabMembBiel - > Valgamma ( ) ; }
stab_utilisable = true ;
}
else if ( ( pt_StabMembBiel - > Aa_calculer ( ) )
& & ( iteg = = 0 )
) // cas ou il faut calculer la référence de stabilisation, mais on ne le fait qu'ici car
// on est sortie de la boucle sur des pti, de manière à utiliser la raideur totale
{ // on récupère la valeur maxi de la raideur
int i , j ; // def
double & maxival = pt_StabMembBiel - > Stab_ref ( ) ; // init
// là il s'agit d'un pointeur donc on travaille directement sur la grandeur stockée
if ( ( enu_stab = = STABMEMBRANE_BIEL_PREMIER_ITER )
| | ( enu_stab = = STABMEMBRANE_BIEL_PREMIER_ITER_NORMALE_AU_NOEUD ) )
{ maxival = ( * raideur ) . MaxiValAbs ( i , j ) / poids_volume ;
// alpha_vient_etre_calculee = true;
}
else if ( ( enu_stab = = STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER )
| | ( enu_stab = = STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_NORMALE_AU_NOEUD ) )
{ maxival = ( * raideur ) . Maxi_ligne_ValAbs ( i ) / poids_volume ;
// alpha_vient_etre_calculee = true;
}
else
{ cout < < " \n *** pb dans la stabilisation de membrane et/ou biel "
< < " le type de stabilisation est " < < Nom_StabMembraneBiel ( enu_stab )
< < " pas pris en compte pour l'instant .... "
< < " on ne peut pas continuer "
< < " \n ElemMeca::Cal_implicit_StabMembBiel(... " < < endl ;
Sortie ( 1 ) ;
} ;
// arrivée ici Stab_ref() a une valeur
// maintenant on va calculer la valeur locale d'alpha
// on récupère la stabilisation pour lui donner une valeur
if ( pt_StabMembBiel - > Pt_fct_gamma ( ) ! = NULL )
// cas où il y a une fonction nD
{ Fonction_nD * fct = pt_StabMembBiel - > Pt_fct_gamma ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = - 1 ; // signifie que l'on utilise la métrique qui vient juste d'être calculée
// ici c'est la métrique correspondante au dernier pti
// comme iteg == 0 on est donc sortie de la boucle sur les pti, le dernier pti c'est nbint
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi ( absolue , TEMPS_tdt , li_enu_scal , nbint , - cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , nbint , - cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// seule la première valeur est ok
alpha = tab_ret ( 1 ) * maxival ;
pt_StabMembBiel - > Valgamma ( ) = alpha ; // on sauvegarde la valeur d'alpha
}
else // sinon cela veut dire que l'intensité de stabilisation est fixe ;
{ alpha = pt_StabMembBiel - > Valgamma ( ) * maxival ; // on prend une proportion du maxi
} ;
stab_utilisable = true ;
}
else if ( iteg = = - 1 ) // cas de la fin d'une boucle, on passe car il n'y a rien à faire
{ if ( pt_StabMembBiel - > Pt_fct_gamma ( ) ! = NULL )
{ alpha = pt_StabMembBiel - > Valgamma ( ) ; } // on récupère la dernière valeur d'alpha calculée
else // sinon il s'agit d'une valeur fixe, on recalcule alpha
{ alpha = pt_StabMembBiel - > Valgamma ( ) * pt_StabMembBiel - > Stab_ref ( ) ; } ;
stab_utilisable = true ;
}
else // ce n'est pas prévu
{ cout < < " \n *** erreur, iteg = " < < iteg //<< " local= "<< local
< < " cas non prevu dans la stabilisation "
< < " \n ElemMeca::Cal_implicit_StabMembBiel(.. " ;
Sortie ( 1 ) ;
} ;
break ;
}
case STABMEMBRANE_BIEL_PREMIER_ITER_INCR : case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR :
case STABMEMBRANE_BIEL_PREMIER_ITER_INCR_NORMALE_AU_NOEUD :
case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR_NORMALE_AU_NOEUD :
{ int compteur_iter = 0 ; // init
// on récupère le pointeur correspondant aux iteration
{ const void * pointe = ( ParaGlob : : param - > GrandeurGlobal ( COMPTEUR_ITERATION_ALGO_GLOBAL ) ) ;
if ( pointe = = NULL )
{ cout < < " \n *** pb dans la stabilisation de membrane et/ou biel "
< < " le type de stabilisation est " < < Nom_StabMembraneBiel ( enu_stab )
< < " et la variable globale COMPTEUR_ITERATION_ALGO_GLOBAL n'est pas disponible "
< < " on ne peut pas continuer "
< < " \n ElemMeca::Cal_implicit_StabMembBiel(... " < < endl ;
Sortie ( 1 ) ;
} ;
// récup du conteneur du compteur
TypeQuelconque * gr_quelc = ( TypeQuelconque * ) ( pointe ) ;
Grandeur_scalaire_entier & gr
= * ( ( Grandeur_scalaire_entier * ) gr_quelc - > Grandeur_pointee ( ) ) ; // pour simplifier
compteur_iter = * ( gr . ConteneurEntier ( ) ) ;
}
// on récupère le pointeur correspondant aux incréments
int compteur_incr = 0 ; // init
{ const void * pointe = ( ParaGlob : : param - > GrandeurGlobal ( COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL ) ) ;
if ( pointe = = NULL )
{ cout < < " \n *** pb dans la stabilisation de membrane et/ou biel "
< < " le type de stabilisation est " < < Nom_StabMembraneBiel ( enu_stab )
< < " et la variable globale COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL n'est pas disponible "
< < " on ne peut pas continuer "
< < " \n ElemMeca::Cal_implicit_StabMembBiel(... " < < endl ;
Sortie ( 1 ) ;
} ;
// récup du conteneur du compteur
TypeQuelconque * gr_quelc = ( TypeQuelconque * ) ( pointe ) ;
Grandeur_scalaire_entier & gr
= * ( ( Grandeur_scalaire_entier * ) gr_quelc - > Grandeur_pointee ( ) ) ; // pour simplifier
compteur_iter = * ( gr . ConteneurEntier ( ) ) ;
}
if ( ( iteg ! = 0 ) & & ( compteur_incr > 1 ) & & ( compteur_iter > 1 ) ) // cela veut dire que la la référence de stabilisation a déjà été calculée
{ alpha = pt_StabMembBiel - > Stab_ref ( ) ; // init de l'intensité de stabilisation
// dans le cas d'une fonction, la partie liée à la raideur est stockée dans Valgamma()
// et il faut appeler la fonction qui peut varier pendant les itérations
if ( pt_StabMembBiel - > Pt_fct_gamma ( ) ! = NULL )
// cas où il y a une fonction nD
{ Fonction_nD * fct = pt_StabMembBiel - > Pt_fct_gamma ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = - 1 ; // signifie que l'on utilise la métrique qui vient juste d'être calculée
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi ( absolue , TEMPS_tdt , li_enu_scal , iteg , - cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , iteg , - cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// seule la première valeur est ok
alpha * = tab_ret ( 1 ) ;
pt_StabMembBiel - > Valgamma ( ) = alpha ; // on sauvegarde la valeur d'alpha
}
else // sinon cela veut dire que l'intensité de stabilisation est fixe
{ alpha * = pt_StabMembBiel - > Valgamma ( ) ; }
stab_utilisable = true ;
}
else if ( ( pt_StabMembBiel - > Aa_calculer ( ) )
& & ( iteg = = 0 )
) // cas ou il faut calculer la référence de stabilisation, mais on ne le fait qu'ici car
// on est sortie de la boucle sur des pti, de manière à utiliser la raideur totale
{ // on récupère la valeur maxi de la raideur
int i , j ; // def
double & maxival = pt_StabMembBiel - > Stab_ref ( ) ; // init
if ( ( enu_stab = = STABMEMBRANE_BIEL_PREMIER_ITER_INCR )
| | ( enu_stab = = STABMEMBRANE_BIEL_PREMIER_ITER_INCR_NORMALE_AU_NOEUD ) )
{ maxival = ( * raideur ) . MaxiValAbs ( i , j ) ;
// alpha_vient_etre_calculee = true;
}
else if ( ( enu_stab = = STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR )
| | ( enu_stab = = STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR_NORMALE_AU_NOEUD ) )
{ maxival = ( * raideur ) . Maxi_ligne_ValAbs ( i ) ;
// alpha_vient_etre_calculee = true;
}
else
{ cout < < " \n *** pb dans la stabilisation de membrane et/ou biel "
< < " le type de stabilisation est " < < Nom_StabMembraneBiel ( enu_stab )
< < " pas pris en compte pour l'instant .... "
< < " on ne peut pas continuer "
< < " \n ElemMeca::Cal_implicit_StabMembBiel(... " < < endl ;
Sortie ( 1 ) ;
} ;
// on récupère la stabilisation
if ( pt_StabMembBiel - > Pt_fct_gamma ( ) ! = NULL )
// cas où il y a une fonction nD
{ Fonction_nD * fct = pt_StabMembBiel - > Pt_fct_gamma ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = - 1 ; // signifie que l'on utilise la métrique qui vient juste d'être calculée
// ici c'est la métrique correspondante au dernier pti
// comme iteg == 0 on est donc sortie de la boucle sur les pti, le dernier pti c'est nbint
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi ( absolue , TEMPS_tdt , li_enu_scal , nbint , - cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , nbint , - cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// seule la première valeur est ok
alpha = tab_ret ( 1 ) * maxival ;
pt_StabMembBiel - > Valgamma ( ) = alpha ; // on sauvegarde la valeur d'alpha
}
else // sinon cela veut dire que l'intensité de stabilisation est fixe ;
{ alpha = pt_StabMembBiel - > Valgamma ( ) * maxival ; // on prend une proportion du maxi
} ;
stab_utilisable = true ;
}
else if ( iteg = = - 1 ) // cas de la fin d'une boucle, on passe car il n'y a rien à faire
{ if ( pt_StabMembBiel - > Pt_fct_gamma ( ) ! = NULL )
{ alpha = pt_StabMembBiel - > Valgamma ( ) ; } // on récupère la dernière valeur d'alpha calculée
else // sinon il s'agit d'une valeur fixe, on recalcule alpha
{ alpha = pt_StabMembBiel - > Valgamma ( ) * pt_StabMembBiel - > Stab_ref ( ) ; } ;
stab_utilisable = true ;
}
else // ce n'est pas prévu
{ cout < < " \n *** erreur, iteg = " < < iteg //<< " local= "<< local
< < " cas non prevu dans la stabilisation "
< < " \n ElemMeca::Cal_implicit_StabMembBiel(.. " ;
Sortie ( 1 ) ;
} ;
break ;
}
/*
// case STAB_M_B_ITER_1_VIA_F_EXT: case STAB_M_B_ITER_1_VIA_F_EXT_NORMALE_AU_NOEUD:
// {// on récupère le pointeur correspondant aux iteration
// const void* pointe = (ParaGlob::param->GrandeurGlobal(COMPTEUR_ITERATION_ALGO_GLOBAL));
// if (pointe == NULL)
// { cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
// << " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
// << "et la variable globale COMPTEUR_ITERATION_ALGO_GLOBAL n'est pas disponible"
// << " on ne peut pas continuer "
// << "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
// Sortie(1);
// };
// // récup du conteneur du compteur
// TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
// Grandeur_scalaire_entier& gr
// = *((Grandeur_scalaire_entier*) gr_quelc->Grandeur_pointee()); // pour simplifier
// int compteur_iter = *(gr.ConteneurEntier());
// if (compteur_iter > 1) // cela veut dire que la force de stabilisation a déjà été calculé
// {// dans le cas d'une fonction, la partie liée à la raiseur est stockée dans Valgamma()
// // et il faut appeler la fonction qui peut varier pendant les itérations
// if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
// // cas où il y a une fonction nD
// { Tableau <double>& tab_val = pt_StabMembBiel->Pt_fct_gamma()->Valeur_pour_variables_globales();
// // seule la première valeur est ok
// alpha = tab_val(1) * pt_StabMembBiel->Valgamma();
// pt_StabMembBiel->FF_StabMembBiel() = alpha; // on sauvegarde
// }
// else // sinon cela veut dire que l'intensité de stabilisation est fixe
// {alpha = pt_StabMembBiel->FF_StabMembBiel();}
// alpha_utilisable=true;
// }
// else if (!local) // sinon il faut la calculer, mais on ne le fait que
// // si on est sortie de la boucle sur des pti,
// {// ***** à faire ***** on récupère la valeur maxi de la raideur
// { cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
// << " partie pas operationnelle pour l'instant !! "
// << " on ne peut pas continuer "
// << "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
// Sortie(1);
// };
// int i,j; // def
// double maxival = 0.; // init
// if ( (enu_stab == STABMEMBRANE_BIEL_PREMIER_ITER)
// ||(enu_stab == STABMEMBRANE_BIEL_PREMIER_ITER_NORMALE_AU_NOEUD))
// {maxival = (*raideur).MaxiValAbs(i,j);
// alpha_vient_etre_calculee = true;
// }
// else if ( (enu_stab == STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER)
// ||(enu_stab == STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_NORMALE_AU_NOEUD))
// {maxival = (*raideur).Maxi_ligne_ValAbs(i);
// alpha_vient_etre_calculee = true;
// }
// else
// { cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
// << " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
// << " pas pris en compte pour l'instant .... "
// << " on ne peut pas continuer "
// << "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
// Sortie(1);
// };
// // on récupère la stabilisation
// alpha=0.; // init
// if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
// // cas où il y a une fonction nD
// { Tableau <double>& tab_val = pt_StabMembBiel->Pt_fct_gamma()->Valeur_pour_variables_globales();
// // seule la première valeur est ok
// alpha = tab_val(1) * maxival;
// pt_StabMembBiel->Valgamma() = maxival; // sauvegarde
// }
// else // sinon cela veut dire que l'intensité de stabilisation est fixe ;
// {alpha = pt_StabMembBiel->Valgamma() * maxival; // on prend une proportion du maxi
// };
// pt_StabMembBiel->FF_StabMembBiel() = alpha; // on sauvegarde
// alpha_utilisable = true;
// };
// break;
// }
*/
default :
cout < < " \n *** erreur le type de stabilisation de membrane et/ou biel n'est pas definie "
< < " on ne peut pas continuer " < < endl ;
Sortie ( 1 ) ;
break ;
} ;
// si la stabilisation est utilisable on s'en sert
if ( stab_utilisable )
// if ( (alpha_utilisable)
// // et soit la stabilisation vient juste d'être calculée, donc on est en dehors de la boucle sur les pti
// // donc on calcule globalement l'impact sur la raideur
// // ou soit on est en local : et on calcul localement l'impact sur la raideur, mais à la sortie de la boucle
// // on ne recalculera pas en globa
// && (alpha_vient_etre_calculee || local)
// )
{ // choix en fonction de la géométrie
Enum_type_geom enu_type_geom = Type_geom_generique ( this - > Id_geometrie ( ) ) ;
switch ( enu_type_geom )
{ case SURFACE :
2023-05-03 17:23:49 +02:00
//#ifdef MISE_AU_POINT
// if (num_elt==25) // -----debug à virer
// {cout << "\n ElemMeca::Cal_implicit_StabMembBiel(..";
// cout << "\n iteg= "<< iteg << flush;
// };
//#endif
2021-09-27 12:42:13 +02:00
if ( ! Contient_Normale_au_noeud ( enu_stab ) )
// cas où on utilise la normale au pti
{ // la normale vaut le produit vectoriel des 2 premiers vecteurs
Coordonnee normale = Util : : ProdVec_coor ( ( * ex . giB_tdt ) . Coordo ( 1 ) , ( * ex . giB_tdt ) . Coordo ( 2 ) ) ;
// calcul de la variation de la normale
// 1) variation du produit vectoriel
Tableau < Coordonnee > D_pasnormale =
Util : : VarProdVect_coor ( ( * ex . giB_tdt ) . Coordo ( 1 ) , ( * ex . giB_tdt ) . Coordo ( 2 ) , ( * ex . d_giB_tdt ) ) ;
// 2) de la normale
Tableau < Coordonnee > D_normale = Util : : VarUnVect_coor ( normale , D_pasnormale , normale . Norme ( ) ) ;
// calcul de la normale normée et pondérée par la stabilisation
normale . Normer ( ) ;
int nbddl = ( * residu ) . Taille ( ) ;
int dimf = 3 ; // ne fonctionne qu'en dim 3
if ( ParaGlob : : AxiSymetrie ( ) )
// cas axisymétrique, dans ce cas on ne prend en compte que les
// dimf-1 coordonnées donc on décrémente
dimf - - ;
int nbne = tab_noeud . Taille ( ) ;
2023-05-03 17:23:49 +02:00
double ponder_F_t = 1. ; // init
2021-09-27 12:42:13 +02:00
// si iteg == 0 ==> un seul calcul, global: correspond au calcul d'alpha: == stab_ref
// si teg == 1 ==> premier calcul correspondant au premier pti
if ( ( iteg = = 0 ) | | ( iteg = = 1 ) )
{ matD - > Initialise ( 0. ) ;
resD - > Inita ( 0. ) ;
maxi_F_t = 0. ;
// on initialise les forces et énergie
int taille = pt_StabMembBiel - > Taille ( ) ;
for ( int i = 1 ; i < = taille ; i + + )
{ pt_StabMembBiel - > FF_StabMembBiel ( i ) = 0. ;
pt_StabMembBiel - > EE_StabMembBiel ( i ) = pt_StabMembBiel - > EE_StabMembBiel_t ( i ) ;
}
} ;
// on récupère l'élément géométrique
ElemGeomC0 & elemGeom = this - > ElementGeometrique ( ) ;
// récup des fonctions d'interpolation
const Tableau < Vecteur > & taphi = elemGeom . TaPhi ( ) ;
2023-05-03 17:23:49 +02:00
int nb_noeud_pris_en_compte = 0 ; //pour faire la moyenne dans le cas iteg == 0
2021-09-27 12:42:13 +02:00
// on applique une force de stabilisation
// uniquement suivant la direction de la normale
int num_pti = iteg ; // init
if ( iteg < = 0 ) num_pti = nbint ;
// la force est proportionelle au déplacement suivant la normale
if ( iteg > - 1 ) // c-a-d == 0 ou plus
{ int ix = 1 ;
for ( int ne = 1 ; ne < = nbne ; ne + + )
{ bool continuer = true ;
if ( noeud_a_prendre_en_compte ! = NULL )
if ( ! ( * noeud_a_prendre_en_compte ) ( ne ) )
continuer = false ;
if ( continuer )
{ Noeud & noe = * ( tab_noeud ( ne ) ) ; // pour simplifier
2023-05-03 17:23:49 +02:00
nb_noeud_pris_en_compte + + ;
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
// on va bloquer le déplacement dans la direction de la normale calculée au pti
2021-09-27 12:42:13 +02:00
//dans un calcul implicit on utilise la position à l'itération 0
// celle-ci peut-être différente de celle à t=0, lorsque l'on a une succession d'algo sur le même incrément
// on commence par récupérer les coordonnées pour l'itération 0
TypeQuelconque_enum_etendu enu ( XI_ITER_0 ) ;
// récupération du conteneur pour lecture uniquement
const TypeQuelconque & typquel = noe . Grandeur_quelconque ( enu ) ;
const Grandeur_coordonnee & cofo = * ( ( Grandeur_coordonnee * ) typquel . Const_Grandeur_pointee ( ) ) ;
Coordonnee delta_X = noe . Coord2 ( ) - cofo . ConteneurCoordonnee_const ( ) ;
2023-05-03 17:23:49 +02:00
double normal_fois_delta_X = normale * delta_X ; // le déplacement au noeud suivant la normale au pti
// dans le cas où on veut pondérer F_t, on calcule la fonction
if ( pt_StabMembBiel - > Pt_fct_ponder_Ft ( ) ! = NULL )
{ Fonction_nD * fct = pt_StabMembBiel - > Pt_fct_ponder_Ft ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = - 1 ; // signifie que l'on utilise la métrique qui vient juste d'être calculée
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi
( absolue , TEMPS_tdt , li_enu_scal , num_pti , - cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , num_pti , - cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// seule la première valeur est ok
ponder_F_t = tab_ret ( 1 ) ;
}
2021-09-27 12:42:13 +02:00
// un - car la force doit s'opposer au déplacement
double intensite = - alpha * normal_fois_delta_X
2023-05-03 17:23:49 +02:00
+ ponder_F_t * pt_StabMembBiel - > FF_StabMembBiel_t ( num_pti ) ;
2021-09-27 12:42:13 +02:00
// on enregistre le maxi de la force à t
maxi_F_t = DabsMaX ( maxi_F_t , pt_StabMembBiel - > FF_StabMembBiel_t ( num_pti ) ) ;
// d'où une variation réelle d'intensité
2023-05-03 17:23:49 +02:00
double delta_intensite = intensite - ponder_F_t * pt_StabMembBiel - > FF_StabMembBiel_t ( num_pti ) ;
2021-09-27 12:42:13 +02:00
pt_StabMembBiel - > EE_StabMembBiel ( num_pti ) + = 0.5 * delta_intensite * normal_fois_delta_X ;
# ifdef MISE_AU_POINT
2023-05-03 17:23:49 +02:00
if ( niveau_affichage > 5 )
// if (num_elt==25) // -----debug à virer
2021-09-27 12:42:13 +02:00
{ cout < < " \n ElemMeca::Cal_implicit_StabMembBiel(.. " ;
2023-05-03 17:23:49 +02:00
cout < < " \n noeud ( " < < ne < < " ), alpha= " < < alpha < < " intensite= " < < intensite
2021-09-27 12:42:13 +02:00
< < " delta_intensite= " < < delta_intensite
< < " delta_X= " < < delta_X
< < " \n normale= " < < normale < < " normal_fois_delta_X= " < < normal_fois_delta_X < < flush ;
} ;
# endif
2023-05-03 17:23:49 +02:00
// //----- debug -----
// if (num_elt==25) // -----debug à virer
// {int nb_pti = pt_StabMembBiel->Taille();
// for (int i=1;i<=nb_pti;i++)
// {cout << "\n avant update F_stab_t("<<i<<")= "<<pt_StabMembBiel->FF_StabMembBiel_t(i)<< " ";
// cout << " F_stab("<<i<<")= "<<pt_StabMembBiel->FF_StabMembBiel(i)<< " ";
// cout << " ponder_F_t= " << ponder_F_t;
// }
// };
// //--- fin debug ----
2021-09-27 12:42:13 +02:00
// calcul de la raideur: on traite les cas :
// iteg >0 et iteg == 0
// le cas iteg == -1 n'est pas à prendre en compte car il correspond
// au cas où on vient de balayer tous les pti
if ( iteg > 0 ) // cas où on est dans une boucle d'intégration
{ // un premier coef pour optimiser
double phi_jacob_poids = taphi ( iteg ) ( ne ) * ( * ex . jacobien_tdt ) * poids_volume ;
2023-05-03 17:23:49 +02:00
2021-09-27 12:42:13 +02:00
pt_StabMembBiel - > FF_StabMembBiel ( iteg ) + = intensite * taphi ( iteg ) ( ne ) ; // * phi_jacob_poids;
// taphi(iteg)(ne) * intensite * (*ex.jacobien_tdt) * poids_volume ;
2023-05-03 17:23:49 +02:00
// NB: 1) comme on multiplie par taphi(iteg)(ne), et que l'on boucle sur les noeuds (ne = 1 à nbne)
// au sortir de la boucle on a l'intensité totale qui est due au déplacement de tous les noeuds
// 2) si alpha * normal_fois_delta_X = 0 c-a-d pas d'accroissement d'intensité, on va obtenir
// somme_noeud(ponder_F_t * pt_StabMembBiel->FF_StabMembBiel_t(num_pti) * taphi(iteg)(ne))
// et comme somme_noeud(taphi(iteg)(ne))=1, on devrait obtenir une intensité globale identique
// à celle à t
2021-09-27 12:42:13 +02:00
pt_StabMembBiel - > EE_StabMembBiel ( iteg ) + = 0.5 * intensite * normal_fois_delta_X * phi_jacob_poids ;
// 0.5 * delta_intensite * normal_fois_delta_X * taphi(iteg)(ne) * (*ex.jacobien_tdt) * poids_volume ;
// un second coef pour optimiser
double phi_poids = taphi ( iteg ) ( ne ) * poids_volume ;
// pour mémoire: intensite = -alpha * normal_fois_delta_X
// + pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
for ( int i = 1 ; i < = dimf ; i + + , ix + + )
{ // = (- normale(i) * alpha * normal_fois_delta_X + F_t)
// * (taphi(iteg)(ne) * poids * jacobien) ;
( * resD ) ( ix ) + = normale ( i ) * intensite * phi_jacob_poids ;
// taphi(iteg)(ne)* normale(i) * intensite
// * (*ex.jacobien_tdt) * poids_volume ; // ok
// contribution à la raideur: de signe inverse à celui du résidu
int j = 1 ;
int me = 1 ; // numéro de noeud dans la direction des ddl
for ( int jx = 1 ; jx < = nbddl ; jx + + , j + + )
{ if ( j > dimf ) { j = 1 ; me + + ; } ;
( * matD ) ( ix , jx ) + = phi_poids // == taphi(iteg)(ne) * poids_volume
* ( ( * ex . jacobien_tdt ) *
( D_normale ( jx ) ( i ) * alpha * normal_fois_delta_X
+ normale ( i ) * alpha * ( D_normale ( jx ) * delta_X )
)
+ ( * ex . d_jacobien_tdt ) ( jx ) * normale ( i ) * ( - intensite )
) ;
if ( ne = = me )
{ ( * matD ) ( ix , jx ) + = phi_poids * normale ( i ) * alpha * normale ( j ) * ( * ex . jacobien_tdt ) ;
// taphi(iteg)(ne)* alpha * normale(i) * normale(j) * poids_volume;
} ;
} ;
} ;
}
else if ( iteg = = 0 ) // on est hors boucle poids_volume == toute la surface
{ // pour le cas particulier d'iteg == 0, on n'a qu'un calcul global
// et on affecte la même intensité à tous les pti
// du coup on sauvegarde uniquement pour le dernier
// ensuite après on répartira sur l'ensemble de pti
// --> en fait on cumulle l'influence de tous les normales
// il faudra ensuite prendre la moyenne
pt_StabMembBiel - > FF_StabMembBiel ( num_pti ) + = intensite ; //* poids_volume;
// taphi(iteg)(ne) * intensite * poids_volume ;
pt_StabMembBiel - > EE_StabMembBiel ( num_pti ) + = 0.5 * intensite * normal_fois_delta_X * poids_volume ;
// 0.5 * intensite * normal_fois_delta_X * taphi(iteg)(ne) * poids_volume ;
for ( int i = 1 ; i < = dimf ; i + + , ix + + )
{ // - normale(i) * alpha * normal_fois_delta_X;
( * resD ) ( ix ) + = normale ( i ) * intensite * poids_volume ; // ok
// contribution à la raideur: de signe inverse à celui du résidu
int j = 1 ;
int me = 1 ; // numéro de noeud dans la direction des ddl
for ( int jx = 1 ; jx < = nbddl ; jx + + , j + + )
{ if ( j > dimf ) { j = 1 ; me + + ; } ;
( * matD ) ( ix , jx ) + = poids_volume
* ( ( D_normale ( jx ) ( i ) * alpha * normal_fois_delta_X
+ normale ( i ) * alpha * ( D_normale ( jx ) * delta_X )
)
) ;
if ( ne = = me )
{ ( * matD ) ( ix , jx ) + = poids_volume * normale ( i ) * alpha * normale ( j ) ;
} ;
} ;
} ;
} ;
} ; // fin du test sur les noeuds à considérer
} ; // fin de la boucle sur les noeuds
} ;
// dans le cas particulier de iteg == 0, c-a-d correspond au calcul d'alpha: == stab_ref
// on n'a pas bouclé sur les pti, c'est uniquement un calcul global qui a été effectué
if ( iteg = = 0 )
2023-05-03 17:23:49 +02:00
{ // on commence par moyenner: se rappeller que l'on a bouclé aussi sur les noeuds !
2021-09-27 12:42:13 +02:00
int taille = pt_StabMembBiel - > Taille ( ) ;
2023-05-03 17:23:49 +02:00
pt_StabMembBiel - > FF_StabMembBiel ( num_pti ) / = ( nb_noeud_pris_en_compte ) ;
pt_StabMembBiel - > EE_StabMembBiel ( num_pti ) / = ( nb_noeud_pris_en_compte ) ;
2021-09-27 12:42:13 +02:00
// ensuite on répartie la même moyenne sur tous les pti
for ( int i = 1 ; i < = taille ; i + + )
{ pt_StabMembBiel - > FF_StabMembBiel ( i ) = pt_StabMembBiel - > FF_StabMembBiel ( num_pti ) ;
pt_StabMembBiel - > EE_StabMembBiel ( i ) = pt_StabMembBiel - > EE_StabMembBiel ( num_pti ) ;
} ;
2023-05-03 17:23:49 +02:00
// également sur le résidu et la raideur, car due à la boucle sur les noeuds
if ( nb_noeud_pris_en_compte > 0 )
{ double inverse = 1. / nb_noeud_pris_en_compte ;
( * matD ) * = inverse ;
( * resD ) * = inverse ;
} ;
2021-09-27 12:42:13 +02:00
} ;
2023-05-03 17:23:49 +02:00
////----- debug -----
// if (num_elt==25) // -----debug à virer
// {int nb_pti = pt_StabMembBiel->Taille();
// for (int i=1;i<=nb_pti;i++)
// cout << "\n final F_stab("<<i<<")= "<<pt_StabMembBiel->FF_StabMembBiel(i)<< " ";
// };
////--- fin debug ----
2021-09-27 12:42:13 +02:00
// on est en dehors de la boucle, on effectue 2 choses
// 1) application d'une limitation éventuelle des efforts
// 2) prise en compte dans la raideur et résidu de l'élément, de la stabilisation
if ( ( iteg = = 0 ) | | ( iteg = = - 1 ) ) // signifie que l'on est en dehors de la boucle
{
2023-05-03 17:23:49 +02:00
if ( pt_StabMembBiel - > Gestion_maxi_mini ( ) )
{
# ifdef MISE_AU_POINT
if ( niveau_affichage > 9 )
// if (num_elt==1) // -----debug à virer
{ cout < < " \n residu de stabilisation avant limitation: " < < ( * resD )
< < " \n residu initial (sans stabilisation): " < < ( * residu )
< < " \n maxi_F_t= " < < maxi_F_t ;
} ;
if ( niveau_affichage > 5 )
// if (num_elt==25) // -----debug à virer
{ int nb_pti = pt_StabMembBiel - > Taille ( ) ;
for ( int i = 1 ; i < = nb_pti ; i + + )
cout < < " \n F_stab( " < < i < < " )= " < < pt_StabMembBiel - > FF_StabMembBiel ( i ) < < " " ;
} ;
# endif
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
// 1) ---- on regarde s'il faut limiter la valeur de la stabilisation
// et on prend en compte dans la raideur globale et le résidu global
// la stabilisation
bool limitation = false ;
double coef = 1. ; // coef de limitation éventuelle si limitation devient true
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
// récup des limitations éventuelles
double beta = pt_StabMembBiel - > Beta ( ) ;
double F_mini = pt_StabMembBiel - > F_mini ( ) ;
double d_maxi = pt_StabMembBiel - > D_maxi ( ) ;
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
// on fait une boucle sur les noeuds
int ix = 1 ;
for ( int ne = 1 ; ne < = nbne ; ne + + )
{ bool continuer = true ;
if ( noeud_a_prendre_en_compte ! = NULL )
if ( ! ( * noeud_a_prendre_en_compte ) ( ne ) )
continuer = false ;
if ( continuer )
{ Noeud & noe = * ( tab_noeud ( ne ) ) ; // pour simplifier
// récup de la force de stabilisation généralisée au noeud
// c'est différent de celle au pti car elle inclue la surface de l'élément
Coordonnee force_stab ( dimf ) ; // init
for ( int i = 1 ; i < = dimf ; i + + , ix + + )
force_stab ( i ) = ( * resD ) ( ix ) ;
double intensite_generalise = force_stab . Norme ( ) ;
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
// on récupère la force externe au noeud à l'instant t
TypeQuelconque_enum_etendu enu_ext ( FORCE_GENE_EXT_t ) ;
// récupération du conteneur pour lecture uniquement
const TypeQuelconque & typquel_ext = noe . Grandeur_quelconque ( enu_ext ) ;
const Grandeur_coordonnee & cofo_ext = * ( ( Grandeur_coordonnee * ) typquel_ext . Const_Grandeur_pointee ( ) ) ;
// on calcule l'intensité de la force externe selon la normale
double intensite_Fext = Dabs ( cofo_ext . ConteneurCoordonnee_const ( ) * normale ) ;
// si l'intensité de la stabilisation est supérieure à beta * intensite_Fext
// on limite
double maxi_intensite = beta * intensite_Fext ;
double max_intens_via_Fext = maxi_intensite ; // on sauve pour le cas où on est dep max
// on applique l'algo théorique cf. doc théorique
// erreur, maxi_intensité est une limitation, mais ce que l'on cherche à voir c'est
// est-ce que les forces externes sont non-nulles, donc c'est intensite_Fext qu'il faut tester
// if (maxi_intensite >= F_mini)
if ( intensite_Fext > = F_mini )
// on est dans le cas où les forces externes sont à prendre en compte
{ if ( ( intensite_generalise > maxi_intensite )
& & ( intensite_generalise > ConstMath : : petit ) // au cas où maxi_intensite et intensite très petites
)
{ limitation = true ;
if ( niveau_affichage > 3 )
{ cout < < " \n limitation stabilisation due aux forces externes, F_stab gene= " < < intensite_generalise
< < " ==> " < < maxi_intensite ;
if ( niveau_affichage > 5 )
cout < < " \n F_mini= " < < F_mini < < " intensite_Fext= " < < intensite_Fext ;
} ;
coef = DabsMiN ( coef , maxi_intensite / intensite_generalise ) ;
} ;
// sinon pas de pb donc on ne fait rien
}
else // cas où l'intensité dans la direction normale, des forces externes n'existent pas
// la limitation va s'exercée via un déplacement maxi
{ maxi_intensite = Dabs ( - alpha * d_maxi + maxi_F_t ) ;
{ if ( ( intensite_generalise > maxi_intensite )
2021-09-27 12:42:13 +02:00
& & ( intensite_generalise > ConstMath : : petit ) // au cas où maxi_intensite et intensite très petites
2023-05-03 17:23:49 +02:00
& & ( Dabs ( alpha ) > ConstMath : : pasmalpetit ) // il faut qu'alpha ne soit pas nul
2021-09-27 12:42:13 +02:00
)
2023-05-03 17:23:49 +02:00
{ limitation = true ;
if ( ( niveau_affichage > 3 ) )
{ cout < < " \n limitation stabilisation due deplace maxi acceptable, F_stab gene= " < < intensite_generalise
< < " ==> " < < maxi_intensite
< < " (max_intens_via_Fext= " < < max_intens_via_Fext < < " ) " ;
if ( ( niveau_affichage > 4 ) )
cout < < " \n maxi_F_t= " < < maxi_F_t < < " alpha*d_maxi= "
< < ( alpha * d_maxi )
< < " alpha= " < < alpha
< < " d_maxi= " < < d_maxi ;
} ;
coef = DabsMiN ( coef , maxi_intensite / intensite_generalise ) ;
} ;
2021-09-27 12:42:13 +02:00
} ;
2023-05-03 17:23:49 +02:00
} ;
} ; // fin du test sur les noeuds à considérer
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
} ; // fin de la boucle sur les noeuds
if ( limitation )
{ ( * resD ) * = coef ; // mise à jour de la force de stabilisation
// idem au niveau des pti
int taille = pt_StabMembBiel - > Taille ( ) ;
for ( int i = 1 ; i < = taille ; i + + )
{ pt_StabMembBiel - > FF_StabMembBiel ( i ) * = coef ;
// cas des énergies
double delta_energie_initiale =
pt_StabMembBiel - > EE_StabMembBiel ( i ) - pt_StabMembBiel - > EE_StabMembBiel_t ( i ) ;
// on coefficiente le delta
double delta_energie = delta_energie_initiale * coef ;
//ce qui donne au final:
pt_StabMembBiel - > EE_StabMembBiel ( i ) =
delta_energie + pt_StabMembBiel - > EE_StabMembBiel_t ( i ) ;
} ;
} ;
2021-09-27 12:42:13 +02:00
} ;
// 2) ---- on met à jour la raideur et résidu global
// en tenant compte de la stabilisation
( * residu ) + = ( * resD ) ;
( * raideur ) + = ( * matD ) ;
//et du coup on ne boucle pas sur les pti et on peut limiter la force sauf qu'une fois le calcul des forces
//aux noeuds fait, il faut ramener par interpolation, la force au pti pour la sauvegarder... bordel
//non... on sauvegarde dans un tableau indicé par les num de noeud, et seulement pour le post-traitement
//on calcul la valeur aux pti
//
# ifdef MISE_AU_POINT
2023-05-03 17:23:49 +02:00
if ( ( niveau_affichage > 9 ) )
// if (num_elt==1) // -----debug à virer
2021-09-27 12:42:13 +02:00
{ cout < < " \n residu de stabilisation: " < < ( * resD ) ;
cout < < " \n raideur de stabilisation: " ; matD - > Affiche ( ) ;
cout < < " \n fin ElemMeca::Cal_implicit_StabMembBiel(.. " ;
} ;
# endif
2023-05-03 17:23:49 +02:00
} ; // fin du cas iteg == 0 ou -1
2021-09-27 12:42:13 +02:00
} // fin du cas où on utilise la normale aux pti
else
{ if ( iteg < 1 )
// on ne calcul que si iteg < 1, c-a-d qu'on est sortie
// de la boucle sur les pti
// sinon cas où on utilise la normale aux noeuds
{ int nbddl = ( * residu ) . Taille ( ) ;
int dimf = 3 ; // ne fonctionne qu'en dim 3
if ( ParaGlob : : AxiSymetrie ( ) )
// cas axisymétrique, dans ce cas on ne prend en compte que les
// dimf-1 coordonnées donc on décrémente
dimf - - ;
int nbne = tab_noeud . Taille ( ) ;
matD - > Initialise ( 0. ) ;
resD - > Inita ( 0. ) ;
maxi_F_t = 0. ;
# ifdef MISE_AU_POINT
// on vérifie le bon dimensionnement
int taille = pt_StabMembBiel - > Taille ( ) ;
if ( taille ! = nbne )
{ cout < < " \n erreur** ElemMeca::Cal_implicit_StabMembBiel(.. "
< < " le stockage est mal dimensionne: taille = " < < taille
< < " alors qu'il devrait etre : nbe = " < < nbne < < flush ;
Sortie ( 1 ) ;
} ;
# endif
// on applique une force de stabilisation aux noeuds
// uniquement suivant la direction de la normale
// la force est proportionelle au déplacement suivant la normale au noeud
int ix = 1 ;
for ( int ne = 1 ; ne < = nbne ; ne + + )
{ Noeud & noe = * ( tab_noeud ( ne ) ) ; // pour simplifier
// on initialise les forces et énergie
pt_StabMembBiel - > FF_StabMembBiel ( ne ) = 0. ;
pt_StabMembBiel - > EE_StabMembBiel ( ne ) = pt_StabMembBiel - > EE_StabMembBiel_t ( ne ) ;
// pour chaque noeud on va bloquer le déplacement dans la direction de la normale au noeud à t
// on récupère la normale au noeud
const TypeQuelconque & tiq = noe . Grandeur_quelconque ( NN_SURF_t ) ;
const Grandeur_coordonnee & gr = * ( ( const Grandeur_coordonnee * ) ( tiq . Const_Grandeur_pointee ( ) ) ) ;
const Coordonnee & normale = gr . ConteneurCoordonnee_const ( ) ;
//dans un calcul implicit on utilise la position à l'itération 0
// celle-ci peut-être différente de celle à t=0, lorsque l'on a une succession d'algo sur le même incrément
// on commence par récupérer les coordonnées pour l'itération 0
TypeQuelconque_enum_etendu enu ( XI_ITER_0 ) ;
// récupération du conteneur pour lecture uniquement
const TypeQuelconque & typquel = noe . Grandeur_quelconque ( enu ) ;
const Grandeur_coordonnee & cofo = * ( ( Grandeur_coordonnee * ) typquel . Const_Grandeur_pointee ( ) ) ;
Coordonnee delta_X = noe . Coord2 ( ) - cofo . ConteneurCoordonnee_const ( ) ;
double normal_fois_delta_X = normale * delta_X ;
2023-05-03 17:23:49 +02:00
double ponder_F_t = 1. ; // init
// dans le cas où on veut pondérer F_t, on calcule la fonction
if ( pt_StabMembBiel - > Pt_fct_ponder_Ft ( ) ! = NULL )
{ Fonction_nD * fct = pt_StabMembBiel - > Pt_fct_ponder_Ft ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = - 1 ; // signifie que l'on utilise la métrique qui vient juste d'être calculée
// au dernier pti !
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi
( absolue , TEMPS_tdt , li_enu_scal , nbint , - cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , nbint , - cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// seule la première valeur est ok
ponder_F_t = tab_ret ( 1 ) ;
}
2021-09-27 12:42:13 +02:00
// un - car la force doit s'opposer au déplacement
double intensite = - alpha * normal_fois_delta_X
2023-05-03 17:23:49 +02:00
+ ponder_F_t * pt_StabMembBiel - > FF_StabMembBiel_t ( ne ) ;
2021-09-27 12:42:13 +02:00
// on enregistre le maxi de la force à t
maxi_F_t = DabsMaX ( maxi_F_t , pt_StabMembBiel - > FF_StabMembBiel_t ( ne ) ) ;
// d'où une variation réelle d'intensité
2023-05-03 17:23:49 +02:00
double delta_intensite = intensite - ponder_F_t * pt_StabMembBiel - > FF_StabMembBiel_t ( ne ) ;
2021-09-27 12:42:13 +02:00
pt_StabMembBiel - > EE_StabMembBiel ( ne ) + = 0.5 * delta_intensite * normal_fois_delta_X ;
# ifdef MISE_AU_POINT
2023-05-03 17:23:49 +02:00
if ( ( niveau_affichage > 5 ) )
2021-09-27 12:42:13 +02:00
{ cout < < " \n ElemMeca::Cal_implicit_StabMembBiel(.. " ;
cout < < " \n alpha= " < < alpha < < " intensite= " < < intensite
< < " delta_intensite= " < < delta_intensite
< < " delta_X= " < < delta_X
< < " \n normale= " < < normale < < " normal_fois_delta_X= " < < normal_fois_delta_X < < flush ;
} ;
# endif
// pour mémoire: intensite = -alpha * normal_fois_delta_X
// + pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
for ( int i = 1 ; i < = dimf ; i + + , ix + + )
{ // - normale(i) * stabili * normal_fois_delta_X;
( * resD ) ( ix ) + = normale ( i ) * intensite ;
// contribution à la raideur: de signe inverse à celui du résidu
int j = 1 ;
int me = 1 ; // numéro de noeud dans la direction des ddl
for ( int jx = 1 ; jx < = nbddl ; jx + + , j + + )
{ if ( j > dimf ) { j = 1 ; me + + ; } ;
if ( ne = = me )
{ ( * matD ) ( ix , jx ) + = alpha * normale ( i ) * normale ( j ) ;
} ;
} ;
} ;
} ; // fin de la boucle sur les noeuds
// on est en dehors de la boucle, on effectue 2 choses
// 1) application d'une limitation éventuelle des efforts
// 2) prise en compte dans la raideur et résidu de l'élément, de la stabilisation
2023-05-03 17:23:49 +02:00
if ( pt_StabMembBiel - > Gestion_maxi_mini ( ) )
{
# ifdef MISE_AU_POINT
if ( ( niveau_affichage > 9 ) )
{ cout < < " \n residu de stabilisation avant limitation: " < < ( * resD )
< < " \n residu initial (sans stabilisation): " < < ( * residu )
< < " \n maxi_F_t= " < < maxi_F_t ;
} ;
if ( ( niveau_affichage > 5 ) )
{ int taille = pt_StabMembBiel - > Taille ( ) ;
for ( int i = 1 ; i < = taille ; i + + )
cout < < " \n F_stab( " < < i < < " )= " < < pt_StabMembBiel - > FF_StabMembBiel ( i ) < < " " ;
} ;
# endif
// 1) ---- on regarde s'il faut limiter la valeur de la stabilisation
// et on prend en compte dans la raideur globale et le résidu global
// la stabilisation
bool limitation = false ;
double coef = 1. ; // coef de limitation éventuelle si limitation devient true
// récup des limitations éventuelles
double beta = pt_StabMembBiel - > Beta ( ) ;
double F_mini = pt_StabMembBiel - > F_mini ( ) ;
double d_maxi = pt_StabMembBiel - > D_maxi ( ) ;
// on refait une boucle sur les noeuds
int ix = 1 ;
for ( int ne = 1 ; ne < = nbne ; ne + + )
{ bool continuer = true ;
if ( noeud_a_prendre_en_compte ! = NULL )
if ( ! ( * noeud_a_prendre_en_compte ) ( ne ) )
continuer = false ;
if ( continuer )
{ Noeud & noe = * ( tab_noeud ( ne ) ) ; // pour simplifier
// récup de la force de stabilisation généralisée au noeud
// c'est différent de celle au pti car elle inclue la surface de l'élément
Coordonnee force_stab ( dimf ) ; // init
for ( int i = 1 ; i < = dimf ; i + + , ix + + )
force_stab ( i ) = ( * resD ) ( ix ) ;
double intensite_generalise = force_stab . Norme ( ) ;
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
// on récupère la force externe au noeud à l'instant t
TypeQuelconque_enum_etendu enu_ext ( FORCE_GENE_EXT_t ) ;
// récupération du conteneur pour lecture uniquement
const TypeQuelconque & typquel_ext = noe . Grandeur_quelconque ( enu_ext ) ;
const Grandeur_coordonnee & cofo_ext = * ( ( Grandeur_coordonnee * ) typquel_ext . Const_Grandeur_pointee ( ) ) ;
// on récupère la normale au noeud
const TypeQuelconque & tiq = noe . Grandeur_quelconque ( NN_SURF_t ) ;
const Grandeur_coordonnee & gr = * ( ( const Grandeur_coordonnee * ) ( tiq . Const_Grandeur_pointee ( ) ) ) ;
const Coordonnee & normale = gr . ConteneurCoordonnee_const ( ) ;
// on calcule l'intensité de la force externe selon la normale
double intensite_Fext = Dabs ( cofo_ext . ConteneurCoordonnee_const ( ) * normale ) ;
// si l'intensité de la stabilisation est supérieure à beta * intensite_Fext
// on limite
double maxi_intensite = beta * intensite_Fext ;
double max_intens_via_Fext = maxi_intensite ; // on sauve pour le cas où on est dep max
// on applique l'algo théorique cf. doc théorique
if ( intensite_Fext > = F_mini )
// on est dans le cas où les forces externes existent
{ if ( ( intensite_generalise > maxi_intensite )
& & ( intensite_generalise > ConstMath : : petit ) // au cas où maxi_intintensite très petites
)
{ limitation = true ;
if ( ( niveau_affichage > 3 ) )
{ cout < < " \n limitation stabilisation due aux forces externes, F_stab gene= " < < intensite_generalise
< < " ==> " < < maxi_intensite ;
if ( ( niveau_affichage > 4 )
| | ( pt_StabMembBiel - > Affichage_stabilisation ( ) > 4 ) )
cout < < " \n F_mini= " < < F_mini < < " intensite_Fext= " < < intensite_Fext ;
} ;
coef = DabsMiN ( coef , maxi_intensite / intensite_generalise ) ;
} ;
// sinon pas de pb donc on ne fait rien
}
else // cas où les forces externes n'existent pas
// la limitation va s'exercée via un déplacement maxi
{ maxi_intensite = Dabs ( - alpha * d_maxi + maxi_F_t ) ;
{ if ( ( intensite_generalise > maxi_intensite )
2021-09-27 12:42:13 +02:00
& & ( intensite_generalise > ConstMath : : petit ) // au cas où maxi_intintensite très petites
2023-05-03 17:23:49 +02:00
& & ( Dabs ( alpha ) > ConstMath : : pasmalpetit ) // il faut qu'alpha ne soit pas nul
2021-09-27 12:42:13 +02:00
)
2023-05-03 17:23:49 +02:00
{ limitation = true ;
if ( ( niveau_affichage > 3 ) )
{ cout < < " \n limitation stabilisation due deplace maxi acceptable, F_s " < < intensite_generalise
< < " ==> " < < maxi_intensite
< < " (max_intens_via_Fext= " < < max_intens_via_Fext < < " ) " ;
if ( ( niveau_affichage > 4 ) )
cout < < " \n maxi_F_t= " < < maxi_F_t < < " alpha*d_maxi= "
< < ( alpha * d_maxi )
< < " alpha= " < < alpha
< < " d_maxi= " < < d_maxi ;
} ;
coef = DabsMiN ( coef , maxi_intensite / intensite_generalise ) ;
} ;
2021-09-27 12:42:13 +02:00
} ;
2023-05-03 17:23:49 +02:00
} ;
} ; // fin du test sur les noeuds à considérer
} ; // fin de la boucle sur les noeuds
if ( limitation )
{ ( * resD ) * = coef ; // mise à jour de la force de stabilisation
// idem au niveau des noeuds
int taille = pt_StabMembBiel - > Taille ( ) ;
for ( int i = 1 ; i < = taille ; i + + )
{ pt_StabMembBiel - > FF_StabMembBiel ( i ) * = coef ;
// cas des énergies
double delta_energie_initiale =
pt_StabMembBiel - > EE_StabMembBiel ( i ) - pt_StabMembBiel - > EE_StabMembBiel_t ( i ) ;
// on coefficiente le delta
double delta_energie = delta_energie_initiale * coef ;
//ce qui donne au final:
pt_StabMembBiel - > EE_StabMembBiel ( i ) =
delta_energie + pt_StabMembBiel - > EE_StabMembBiel_t ( i ) ;
2021-09-27 12:42:13 +02:00
} ;
2023-05-03 17:23:49 +02:00
// 2) ---- on met à jour la raideur et résidu global
// en tenant compte de la stabilisation
( * residu ) + = ( * resD ) ;
( * raideur ) + = ( * matD ) ;
} ; // fin du cas if (limitation)
2021-09-27 12:42:13 +02:00
} ; // fin du cas ou on applique éventuellement une limitation
# ifdef MISE_AU_POINT
2023-05-03 17:23:49 +02:00
if ( ( niveau_affichage > 9 ) )
2021-09-27 12:42:13 +02:00
{ cout < < " \n residu de stabilisation: " < < ( * resD ) ;
cout < < " \n raideur de stabilisation: " ; matD - > Affiche ( ) ;
cout < < " \n fin debug ElemMeca::Cal_implicit_StabMembBiel(.. " ;
} ;
# endif
} // fin du cas où on utilise la normale aux noeuds
} ;
break ;
case LIGNE :
break ;
default : // on ne fait rien pour l'instant pour les autres types
break ;
} ;
} ;
} ;
// stabilisation pour un calcul explicit
void ElemMeca : : Cal_explicit_StabMembBiel ( int iteg , const Met_abstraite : : Expli_t_tdt ex
, const int & nbint , const double & poids_volume
, Tableau < int > * noeud_a_prendre_en_compte )
{
// la stabilisation en explicite n'est possible que s'il existe soit:
// - une force de stabilisation déjà enregistrée
// - un coefficient alpha non nulle que l'on peut calculer avec les données déjà enregistrées
// l'idée est de bloquer la direction normale à l'élément pour une plaque et
// les deux directions transverses pour une biel
2023-05-03 17:23:49 +02:00
int niveau_affichage = ParaGlob : : NiveauImpression ( ) ;
if ( pt_StabMembBiel - > Affichage_stabilisation ( ) > 0 )
niveau_affichage = pt_StabMembBiel - > Affichage_stabilisation ( ) ;
2021-09-27 12:42:13 +02:00
//---- détermination de l'intensité de stabilisation
// ici on ne s'occupe pas de l'itération ni de l'incrément, dans tous les cas
// on utilise les grandeurs déjà stockées
// ceci quelque soit le type de stabilisation
double alpha = 0. ; // init
bool stab_utilisable = false ; // init
Enum_StabMembraneBiel enu_stab = pt_StabMembBiel - > Type_stabMembrane ( ) ;
2023-05-03 17:23:49 +02:00
2021-09-27 12:42:13 +02:00
if ( pt_StabMembBiel ! = NULL )
2023-05-03 17:23:49 +02:00
if ( pt_StabMembBiel - > Activite_en_explicite ( ) )
{ //alpha = pt_StabMembBiel->Stab_ref(); // récup de l'intensité de stabilisation
2021-09-27 12:42:13 +02:00
// ce n'est pas un pointeur, alpha est une variable locale
2023-05-03 17:23:49 +02:00
// stab_ref, par défaut == 1, mais dans le cas d'un passage
// implicit -> explicite : vaut la dernière valeur implicite calculée
// cela permet de profiter du passage en implicite, en récupérant l'intensité
// de raideur réelle
alpha = pt_StabMembBiel - > Stab_ref ( ) ; // init
2021-09-27 12:42:13 +02:00
// dans le cas d'une fonction multiplicatrice ,
2023-05-03 17:23:49 +02:00
// et il faut appeler la fonction qui peut varier
2021-09-27 12:42:13 +02:00
if ( pt_StabMembBiel - > Pt_fct_gamma ( ) ! = NULL )
// cas où il y a une fonction nD
{ Fonction_nD * fct = pt_StabMembBiel - > Pt_fct_gamma ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = - 1 ; // signifie que l'on utilise la métrique qui vient juste d'être calculée
2023-05-03 17:23:49 +02:00
int num_iteg_a_considerer ; // init
if ( iteg > 0 ) // on est dans la boucle
num_iteg_a_considerer = iteg ;
else // sinon on est à la sortie de la boucle
num_iteg_a_considerer = nbint ;
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi
( absolue , TEMPS_tdt , li_enu_scal , num_iteg_a_considerer , - cas ) ) ;
2021-09-27 12:42:13 +02:00
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
2023-05-03 17:23:49 +02:00
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , num_iteg_a_considerer , - cas ) ;
2021-09-27 12:42:13 +02:00
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// seule la première valeur est ok
alpha * = tab_ret ( 1 ) ;
pt_StabMembBiel - > Valgamma ( ) = alpha ; // on sauvegarde la valeur d'alpha
}
else // sinon cela veut dire que l'intensité de stabilisation est fixe
2023-05-03 17:23:49 +02:00
// en explicite c'est directement la valeur lue
{ alpha * = pt_StabMembBiel - > Valgamma ( ) ; }
2021-09-27 12:42:13 +02:00
stab_utilisable = true ;
} ;
// si la stabilisation est utilisable on s'en sert
if ( stab_utilisable )
{ // choix en fonction de la géométrie
Enum_type_geom enu_type_geom = Type_geom_generique ( this - > Id_geometrie ( ) ) ;
switch ( enu_type_geom )
{ case SURFACE :
if ( ! Contient_Normale_au_noeud ( enu_stab ) )
// --- cas où on utilise la normale au pti
{ // la normale vaut le produit vectoriel des 2 premiers vecteurs
Coordonnee normale = Util : : ProdVec_coor ( ( * ex . giB_tdt ) . Coordo ( 1 ) , ( * ex . giB_tdt ) . Coordo ( 2 ) ) ;
// calcul de la normale normée et pondérée par la stabilisation
normale . Normer ( ) ;
int nbddl = ( * residu ) . Taille ( ) ;
int dimf = 3 ; // ne fonctionne qu'en dim 3
if ( ParaGlob : : AxiSymetrie ( ) )
// cas axisymétrique, dans ce cas on ne prend en compte que les
// dimf-1 coordonnées donc on décrémente
dimf - - ;
int nbne = tab_noeud . Taille ( ) ;
2023-05-03 17:23:49 +02:00
double ponder_F_t = 1. ; // init
// si iteg == 1 ==> premier calcul correspondant au premier pti
2021-09-27 12:42:13 +02:00
if ( iteg = = 1 )
{ resD - > Inita ( 0. ) ;
maxi_F_t = 0. ;
// on initialise les forces et énergie
int taille = pt_StabMembBiel - > Taille ( ) ;
for ( int i = 1 ; i < = taille ; i + + )
{ pt_StabMembBiel - > FF_StabMembBiel ( i ) = 0. ;
pt_StabMembBiel - > EE_StabMembBiel ( i ) = pt_StabMembBiel - > EE_StabMembBiel_t ( i ) ;
}
} ;
// on récupère l'élément géométrique
ElemGeomC0 & elemGeom = this - > ElementGeometrique ( ) ;
// récup des fonctions d'interpolation
const Tableau < Vecteur > & taphi = elemGeom . TaPhi ( ) ;
// on applique une force de stabilisation
// uniquement suivant la direction de la normale
// ici iteg peut soit == -1 : on est sortie de la boucle sur les pti
// soit > 0 : on et dans la boucle des pti
int num_pti = iteg ; // init
if ( iteg = = - 1 ) num_pti = nbint ;
// la force est proportionelle au déplacement suivant la normale
if ( iteg > 0 )
{ int ix = 1 ;
for ( int ne = 1 ; ne < = nbne ; ne + + )
{ bool continuer = true ;
if ( noeud_a_prendre_en_compte ! = NULL )
if ( ! ( * noeud_a_prendre_en_compte ) ( ne ) )
continuer = false ;
if ( continuer )
{ Noeud & noe = * ( tab_noeud ( ne ) ) ; // pour simplifier
// pour chaque noeud on va rectifier la normale par rapport à celle calculée au pti
// on va bloquer le déplacement dans la direction de la normale au noeud à t
// celle-ci peut-être différente de celle à t=0, lorsque l'on a une succession d'algo sur le même incrément
// on commence par récupérer les coordonnées pour l'itération 0
TypeQuelconque_enum_etendu enu ( XI_ITER_0 ) ;
// récupération du conteneur pour lecture uniquement
const TypeQuelconque & typquel = noe . Grandeur_quelconque ( enu ) ;
const Grandeur_coordonnee & cofo = * ( ( Grandeur_coordonnee * ) typquel . Const_Grandeur_pointee ( ) ) ;
Coordonnee delta_X = noe . Coord2 ( ) - cofo . ConteneurCoordonnee_const ( ) ;
double normal_fois_delta_X = normale * delta_X ; // le déplacement suivant la normale
2023-05-03 17:23:49 +02:00
// dans le cas où on veut pondérer F_t, on calcule la fonction
if ( pt_StabMembBiel - > Pt_fct_ponder_Ft ( ) ! = NULL )
{ Fonction_nD * fct = pt_StabMembBiel - > Pt_fct_ponder_Ft ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = - 1 ; // signifie que l'on utilise la métrique qui vient juste d'être calculée
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi
( absolue , TEMPS_tdt , li_enu_scal , num_pti , - cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , num_pti , - cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// seule la première valeur est ok
ponder_F_t = tab_ret ( 1 ) ;
}
2021-09-27 12:42:13 +02:00
// un - car la force doit s'opposer au déplacement
double intensite = - alpha * normal_fois_delta_X
2023-05-03 17:23:49 +02:00
+ ponder_F_t * pt_StabMembBiel - > FF_StabMembBiel_t ( num_pti ) ;
2021-09-27 12:42:13 +02:00
// on enregistre le maxi de la force à t
maxi_F_t = DabsMaX ( maxi_F_t , pt_StabMembBiel - > FF_StabMembBiel_t ( num_pti ) ) ;
// d'où une variation réelle d'intensité
2023-05-03 17:23:49 +02:00
double delta_intensite = intensite - ponder_F_t * pt_StabMembBiel - > FF_StabMembBiel_t ( num_pti ) ;
2021-09-27 12:42:13 +02:00
pt_StabMembBiel - > EE_StabMembBiel ( num_pti ) + = 0.5 * delta_intensite * normal_fois_delta_X ;
# ifdef MISE_AU_POINT
2023-05-03 17:23:49 +02:00
if ( ( niveau_affichage > 5 ) )
2021-09-27 12:42:13 +02:00
{ cout < < " \n ElemMeca::Cal_explicit_StabMembBiel(.. " ;
cout < < " \n alpha= " < < alpha < < " intensite= " < < intensite
< < " delta_intensite= " < < delta_intensite
< < " delta_X= " < < delta_X
< < " \n normale= " < < normale < < " normal_fois_delta_X= " < < normal_fois_delta_X < < flush ;
} ;
# endif
// calcul de la raideur: on traite les cas :
// iteg >0 et iteg == 0
// le cas iteg == -1 n'est pas à prendre en compte car il correspond
// au cas où on vient de balayer tous les pti
if ( iteg > 0 ) // cas où on est dans une boucle d'intégration
{ // un premier coef pour optimiser
double phi_jacob_poids = taphi ( iteg ) ( ne ) * ( * ex . jacobien_tdt ) * poids_volume ;
pt_StabMembBiel - > FF_StabMembBiel ( iteg ) + = intensite * taphi ( iteg ) ( ne ) ; // * phi_jacob_poids;
// taphi(iteg)(ne) * intensite * (*ex.jacobien_tdt) * poids_volume ;
pt_StabMembBiel - > EE_StabMembBiel ( iteg ) + = 0.5 * intensite * normal_fois_delta_X * phi_jacob_poids ;
// 0.5 * delta_intensite * normal_fois_delta_X * taphi(iteg)(ne) * (*ex.jacobien_tdt) * poids_volume ;
// un second coef pour optimiser
2023-05-03 17:23:49 +02:00
// double phi_poids = taphi(iteg)(ne) * poids_volume ;
2021-09-27 12:42:13 +02:00
// pour mémoire: intensite = -alpha * normal_fois_delta_X
// + pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
for ( int i = 1 ; i < = dimf ; i + + , ix + + )
{ // = (- normale(i) * alpha * normal_fois_delta_X + F_t)
// * (taphi(iteg)(ne) * poids * jacobien) ;
( * resD ) ( ix ) + = normale ( i ) * intensite * phi_jacob_poids ;
// taphi(iteg)(ne)* normale(i) * intensite
// * (*ex.jacobien_tdt) * poids_volume ; // ok
} ;
}
} ; // fin du test sur les noeuds à considérer
} ; // fin de la boucle sur les noeuds
} ;
// on est en dehors de la boucle, on effectue 2 choses
// 1) application d'une limitation éventuelle des efforts
// 2) prise en compte dans la raideur et résidu de l'élément, de la stabilisation
// ici seul le cas iteg == -1 est pris en compte contrairement à l'implicite
if ( iteg = = - 1 ) // signifie que l'on est en dehors de la boucle
{
2023-05-03 17:23:49 +02:00
if ( pt_StabMembBiel - > Gestion_maxi_mini ( ) )
{
# ifdef MISE_AU_POINT
if ( ( niveau_affichage > 5 ) )
{ cout < < " \n residu de stabilisation avant limitation: " < < ( * resD )
< < " \n residu initial (sans stabilisation): " < < ( * residu )
< < " \n maxi_F_t= " < < maxi_F_t ;
} ;
if ( ( niveau_affichage > 5 ) )
{ int nb_pti = pt_StabMembBiel - > Taille ( ) ;
for ( int i = 1 ; i < = nb_pti ; i + + )
cout < < " \n F_stab( " < < i < < " )= " < < pt_StabMembBiel - > FF_StabMembBiel ( i ) < < " " ;
} ;
# endif
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
// 1) ---- on regarde s'il faut limiter la valeur de la stabilisation
// et on prend en compte dans la raideur globale et le résidu global
// la stabilisation
bool limitation = false ;
double coef = 1. ; // coef de limitation éventuelle si limitation devient true
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
// récup des limitations éventuelles
double beta = pt_StabMembBiel - > Beta ( ) ;
double F_mini = pt_StabMembBiel - > F_mini ( ) ;
double d_maxi = pt_StabMembBiel - > D_maxi ( ) ;
// on refait une boucle sur les pti
int ix = 1 ;
for ( int ne = 1 ; ne < = nbne ; ne + + )
{ bool continuer = true ;
if ( noeud_a_prendre_en_compte ! = NULL )
if ( ! ( * noeud_a_prendre_en_compte ) ( ne ) )
continuer = false ;
if ( continuer )
{ Noeud & noe = * ( tab_noeud ( ne ) ) ; // pour simplifier
// récup de la force de stabilisation généralisée au noeud
// c'est différent de celle au pti car elle inclue la surface de l'élément
Coordonnee force_stab ( dimf ) ; // init
for ( int i = 1 ; i < = dimf ; i + + , ix + + )
force_stab ( i ) = ( * resD ) ( ix ) ;
double intensite_generalise = force_stab . Norme ( ) ;
// on récupère la force externe au noeud à l'instant t
TypeQuelconque_enum_etendu enu_ext ( FORCE_GENE_EXT_t ) ;
// récupération du conteneur pour lecture uniquement
const TypeQuelconque & typquel_ext = noe . Grandeur_quelconque ( enu_ext ) ;
const Grandeur_coordonnee & cofo_ext = * ( ( Grandeur_coordonnee * ) typquel_ext . Const_Grandeur_pointee ( ) ) ;
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
// on calcule l'intensité de la force externe selon la normale
double intensite_Fext = Dabs ( cofo_ext . ConteneurCoordonnee_const ( ) * normale ) ;
// si l'intensité de la stabilisation est supérieure à beta * intensite_Fext
// on limite
double maxi_intensite = beta * intensite_Fext ;
double max_intens_via_Fext = maxi_intensite ; // on sauve pour le cas où on est dep max
// on applique l'algo théorique cf. doc théorique
if ( intensite_Fext > = F_mini )
// on est dans le cas où les forces externes existent
{ if ( ( intensite_generalise > maxi_intensite )
& & ( intensite_generalise > ConstMath : : petit ) // au cas où maxi_intensite et intensite très petites
)
{ limitation = true ;
if ( ( niveau_affichage > 3 ) )
{ cout < < " \n limitation stabilisation due aux forces externes, F_stab gene= " < < intensite_generalise
< < " ==> " < < maxi_intensite ;
if ( ( niveau_affichage > 5 ) )
cout < < " \n F_mini= " < < F_mini < < " intensite_Fext= " < < intensite_Fext ;
} ;
coef = DabsMiN ( coef , maxi_intensite / intensite_generalise ) ;
2021-09-27 12:42:13 +02:00
} ;
2023-05-03 17:23:49 +02:00
// sinon pas de pb donc on ne fait rien
}
else // cas où les forces externes n'existent pas
// la limitation va s'exercée via un déplacement maxi
{ maxi_intensite = Dabs ( - alpha * d_maxi + maxi_F_t ) ;
{ if ( ( intensite_generalise > maxi_intensite )
& & ( intensite_generalise > ConstMath : : petit ) // au cas où maxi_intensite et intensite très petites
& & ( Dabs ( alpha ) > ConstMath : : pasmalpetit ) // il faut qu'alpha ne soit pas nul
)
{ limitation = true ;
if ( ( niveau_affichage > 3 ) )
{ cout < < " \n limitation stabilisation due deplace maxi acceptable, F_stab gene= " < < intensite_generalise
< < " ==> " < < maxi_intensite
< < " (max_intens_via_Fext= " < < max_intens_via_Fext < < " ) " ;
if ( ( niveau_affichage > 4 ) )
cout < < " \n maxi_F_t= " < < maxi_F_t < < " alpha*d_maxi= "
< < ( alpha * d_maxi )
< < " alpha= " < < alpha
< < " d_maxi= " < < d_maxi ;
} ;
coef = DabsMiN ( coef , maxi_intensite / intensite_generalise ) ;
2021-09-27 12:42:13 +02:00
} ;
2023-05-03 17:23:49 +02:00
} ;
2021-09-27 12:42:13 +02:00
} ;
2023-05-03 17:23:49 +02:00
} ; // fin du test sur les noeuds à considérer
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
} ; // fin de la boucle sur les noeuds
if ( limitation )
{ ( * resD ) * = coef ; // mise à jour de la force de stabilisation
// idem au niveau des pti
int taille = pt_StabMembBiel - > Taille ( ) ;
for ( int i = 1 ; i < = taille ; i + + )
{ pt_StabMembBiel - > FF_StabMembBiel ( i ) * = coef ;
// cas des énergies
double delta_energie_initiale =
pt_StabMembBiel - > EE_StabMembBiel ( i ) - pt_StabMembBiel - > EE_StabMembBiel_t ( i ) ;
// on coefficiente le delta
double delta_energie = delta_energie_initiale * coef ;
//ce qui donne au final:
pt_StabMembBiel - > EE_StabMembBiel ( i ) =
delta_energie + pt_StabMembBiel - > EE_StabMembBiel_t ( i ) ;
} ;
} ;
} ; // fin mise en place d'une limitation éventuelle
2021-09-27 12:42:13 +02:00
// 2) ---- on met à jour le résidu global
// en tenant compte de la stabilisation
( * residu ) + = ( * resD ) ;
# ifdef MISE_AU_POINT
2023-05-03 17:23:49 +02:00
if ( ( niveau_affichage > 9 ) )
2021-09-27 12:42:13 +02:00
{ cout < < " \n residu de stabilisation: " < < ( * resD ) ;
cout < < " \n fin ElemMeca::Cal_explicit_StabMembBiel(.. " ;
} ;
# endif
} ; // fin du cas iteg == -1
} // fin du cas où on utilise la normale aux pti
else //--- cas où on utilise la normale au noeud
{ if ( iteg = = - 1 )
// on ne calcul que si iteg == -1, c-a-d qu'on est sortie
// de la boucle sur les pti
// sinon cas où on utilise la normale aux noeuds
{ int nbddl = ( * residu ) . Taille ( ) ;
int dimf = 3 ; // ne fonctionne qu'en dim 3
if ( ParaGlob : : AxiSymetrie ( ) )
// cas axisymétrique, dans ce cas on ne prend en compte que les
// dimf-1 coordonnées donc on décrémente
dimf - - ;
int nbne = tab_noeud . Taille ( ) ;
resD - > Inita ( 0. ) ;
maxi_F_t = 0. ;
# ifdef MISE_AU_POINT
// on vérifie le bon dimensionnement
int taille = pt_StabMembBiel - > Taille ( ) ;
if ( taille ! = nbne )
{ cout < < " \n erreur** ElemMeca::Cal_explicit_StabMembBiel(.. "
< < " le stockage est mal dimensionne: taille = " < < taille
< < " alors qu'il devrait etre : nbe = " < < nbne < < flush ;
Sortie ( 1 ) ;
} ;
# endif
// on applique une force de stabilisation aux noeuds
// uniquement suivant la direction de la normale
// la force est proportionelle au déplacement suivant la normale au noeud
int ix = 1 ;
for ( int ne = 1 ; ne < = nbne ; ne + + )
{ Noeud & noe = * ( tab_noeud ( ne ) ) ; // pour simplifier
// on initialise les forces et énergie
pt_StabMembBiel - > FF_StabMembBiel ( ne ) = 0. ;
pt_StabMembBiel - > EE_StabMembBiel ( ne ) = pt_StabMembBiel - > EE_StabMembBiel_t ( ne ) ;
// pour chaque noeud on va bloquer le déplacement dans la direction de la normale au noeud à t
// on récupère la normale au noeud
const TypeQuelconque & tiq = noe . Grandeur_quelconque ( NN_SURF_t ) ;
const Grandeur_coordonnee & gr = * ( ( const Grandeur_coordonnee * ) ( tiq . Const_Grandeur_pointee ( ) ) ) ;
const Coordonnee & normale = gr . ConteneurCoordonnee_const ( ) ;
//dans un calcul implicit on utilise la position à l'itération 0
// celle-ci peut-être différente de celle à t=0, lorsque l'on a une succession d'algo sur le même incrément
// on commence par récupérer les coordonnées pour l'itération 0
TypeQuelconque_enum_etendu enu ( XI_ITER_0 ) ;
// récupération du conteneur pour lecture uniquement
const TypeQuelconque & typquel = noe . Grandeur_quelconque ( enu ) ;
const Grandeur_coordonnee & cofo = * ( ( Grandeur_coordonnee * ) typquel . Const_Grandeur_pointee ( ) ) ;
Coordonnee delta_X = noe . Coord2 ( ) - cofo . ConteneurCoordonnee_const ( ) ;
double normal_fois_delta_X = normale * delta_X ;
2023-05-03 17:23:49 +02:00
double ponder_F_t = 1. ; // init
// dans le cas où on veut pondérer F_t, on calcule la fonction
if ( pt_StabMembBiel - > Pt_fct_ponder_Ft ( ) ! = NULL )
{ Fonction_nD * fct = pt_StabMembBiel - > Pt_fct_ponder_Ft ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = - 1 ; // signifie que l'on utilise la métrique qui vient juste d'être calculée
// au dernier pti !
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi
( absolue , TEMPS_tdt , li_enu_scal , nbint , - cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , nbint , - cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// seule la première valeur est ok
ponder_F_t = tab_ret ( 1 ) ;
}
2021-09-27 12:42:13 +02:00
// un - car la force doit s'opposer au déplacement
double intensite = - alpha * normal_fois_delta_X
2023-05-03 17:23:49 +02:00
+ ponder_F_t * pt_StabMembBiel - > FF_StabMembBiel_t ( ne ) ;
2021-09-27 12:42:13 +02:00
// on enregistre le maxi de la force à t
maxi_F_t = DabsMaX ( maxi_F_t , pt_StabMembBiel - > FF_StabMembBiel_t ( ne ) ) ;
// d'où une variation réelle d'intensité
2023-05-03 17:23:49 +02:00
double delta_intensite = intensite - ponder_F_t * pt_StabMembBiel - > FF_StabMembBiel_t ( ne ) ;
2021-09-27 12:42:13 +02:00
pt_StabMembBiel - > EE_StabMembBiel ( ne ) + = 0.5 * delta_intensite * normal_fois_delta_X ;
# ifdef MISE_AU_POINT
2023-05-03 17:23:49 +02:00
if ( ( niveau_affichage > 5 ) )
2021-09-27 12:42:13 +02:00
{ cout < < " \n ElemMeca::Cal_explicit_StabMembBiel(.. " ;
cout < < " \n alpha= " < < alpha < < " intensite= " < < intensite
< < " delta_intensite= " < < delta_intensite
< < " delta_X= " < < delta_X
< < " \n normale= " < < normale < < " normal_fois_delta_X= " < < normal_fois_delta_X < < flush ;
} ;
# endif
// pour mémoire: intensite = -alpha * normal_fois_delta_X
// + pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
for ( int i = 1 ; i < = dimf ; i + + , ix + + )
{ // - normale(i) * stabili * normal_fois_delta_X;
( * resD ) ( ix ) + = normale ( i ) * intensite ;
} ;
} ; // fin de la boucle sur les noeuds
// on est en dehors de la boucle, on effectue 2 choses
// 1) application d'une limitation éventuelle des efforts
// 2) prise en compte dans la raideur et résidu de l'élément, de la stabilisation
2023-05-03 17:23:49 +02:00
if ( pt_StabMembBiel - > Gestion_maxi_mini ( ) )
{
# ifdef MISE_AU_POINT
if ( ( niveau_affichage > 4 ) )
{ cout < < " \n residu de stabilisation avant limitation: " < < ( * resD )
< < " \n residu initial (sans stabilisation): " < < ( * residu )
< < " \n maxi_F_t= " < < maxi_F_t ;
} ;
if ( ( niveau_affichage > 4 ) )
{ int taille = pt_StabMembBiel - > Taille ( ) ;
for ( int i = 1 ; i < = taille ; i + + )
cout < < " \n F_stab( " < < i < < " )= " < < pt_StabMembBiel - > FF_StabMembBiel ( i ) < < " " ;
} ;
# endif
// 1) ---- on regarde s'il faut limiter la valeur de la stabilisation
// et on prend en compte dans la raideur globale et le résidu global
// la stabilisation
bool limitation = false ;
double coef = 1. ; // coef de limitation éventuelle si limitation devient true
// récup des limitations éventuelles
double beta = pt_StabMembBiel - > Beta ( ) ;
double F_mini = pt_StabMembBiel - > F_mini ( ) ;
double d_maxi = pt_StabMembBiel - > D_maxi ( ) ;
// on refait une boucle sur les noeuds
int ix = 1 ;
for ( int ne = 1 ; ne < = nbne ; ne + + )
{ bool continuer = true ;
if ( noeud_a_prendre_en_compte ! = NULL )
if ( ! ( * noeud_a_prendre_en_compte ) ( ne ) )
continuer = false ;
if ( continuer )
{ Noeud & noe = * ( tab_noeud ( ne ) ) ; // pour simplifier
// récup de la force de stabilisation généralisée au noeud
// c'est différent de celle au pti car elle inclue la surface de l'élément
Coordonnee force_stab ( dimf ) ; // init
for ( int i = 1 ; i < = dimf ; i + + , ix + + )
force_stab ( i ) = ( * resD ) ( ix ) ;
double intensite_generalise = force_stab . Norme ( ) ;
// on récupère la force externe au noeud à l'instant t
TypeQuelconque_enum_etendu enu_ext ( FORCE_GENE_EXT_t ) ;
// récupération du conteneur pour lecture uniquement
const TypeQuelconque & typquel_ext = noe . Grandeur_quelconque ( enu_ext ) ;
const Grandeur_coordonnee & cofo_ext = * ( ( Grandeur_coordonnee * ) typquel_ext . Const_Grandeur_pointee ( ) ) ;
// on récupère la normale au noeud
const TypeQuelconque & tiq = noe . Grandeur_quelconque ( NN_SURF_t ) ;
const Grandeur_coordonnee & gr = * ( ( const Grandeur_coordonnee * ) ( tiq . Const_Grandeur_pointee ( ) ) ) ;
const Coordonnee & normale = gr . ConteneurCoordonnee_const ( ) ;
// on calcule l'intensité de la force externe selon la normale
double intensite_Fext = Dabs ( cofo_ext . ConteneurCoordonnee_const ( ) * normale ) ;
// si l'intensité de la stabilisation est supérieure à beta * intensite_Fext
// on limite
double maxi_intensite = beta * intensite_Fext ;
double max_intens_via_Fext = maxi_intensite ; // on sauve pour le cas où on est en dep max
// on applique l'algo théorique cf. doc théorique
if ( intensite_Fext > = F_mini )
// on est dans le cas où les forces externes existent
{ if ( ( intensite_generalise > maxi_intensite )
& & ( intensite_generalise > ConstMath : : petit ) // au cas où maxi_intintensite très petites
)
{ limitation = true ;
if ( ( niveau_affichage > 3 ) )
{ cout < < " \n limitation stabilisation due aux forces externes, F_stab gene= " < < intensite_generalise
< < " ==> " < < maxi_intensite ;
if ( ( niveau_affichage > 4 ) )
cout < < " \n F_mini= " < < F_mini < < " intensite_Fext= " < < intensite_Fext ;
} ;
coef = DabsMiN ( coef , maxi_intensite / intensite_generalise ) ;
2021-09-27 12:42:13 +02:00
} ;
2023-05-03 17:23:49 +02:00
// sinon pas de pb donc on ne fait rien
}
else // cas où les forces externes n'existent pas
// la limitation va s'exercée via un déplacement maxi
{ maxi_intensite = Dabs ( - alpha * d_maxi + maxi_F_t ) ;
{ if ( ( intensite_generalise > maxi_intensite )
& & ( intensite_generalise > ConstMath : : petit ) // au cas où maxi_intintensite très petites
& & ( Dabs ( alpha ) > ConstMath : : pasmalpetit ) // il faut qu'alpha ne soit pas nul
)
{ limitation = true ;
if ( ( niveau_affichage > 3 ) )
{ cout < < " \n limitation stabilisation due deplace maxi acceptable, F_s " < < intensite_generalise
< < " ==> " < < maxi_intensite < < " (max_intens_via_Fext= " < < max_intens_via_Fext < < " ) " ;
if ( ( niveau_affichage > 4 ) )
cout < < " \n maxi_F_t= " < < maxi_F_t < < " alpha*d_maxi= "
< < ( alpha * d_maxi )
< < " alpha= " < < alpha
< < " d_maxi= " < < d_maxi ;
} ;
coef = DabsMiN ( coef , maxi_intensite / intensite_generalise ) ;
2021-09-27 12:42:13 +02:00
} ;
2023-05-03 17:23:49 +02:00
} ;
2021-09-27 12:42:13 +02:00
} ;
2023-05-03 17:23:49 +02:00
} ; // fin du test sur les noeuds à considérer
} ; // fin de la boucle sur les noeuds
if ( limitation )
{ ( * resD ) * = coef ; // mise à jour de la force de stabilisation
// idem au niveau des noeuds
int taille = pt_StabMembBiel - > Taille ( ) ;
for ( int i = 1 ; i < = taille ; i + + )
{ pt_StabMembBiel - > FF_StabMembBiel ( i ) * = coef ;
// cas des énergies
double delta_energie_initiale =
pt_StabMembBiel - > EE_StabMembBiel ( i ) - pt_StabMembBiel - > EE_StabMembBiel_t ( i ) ;
// on coefficiente le delta
double delta_energie = delta_energie_initiale * coef ;
//ce qui donne au final:
pt_StabMembBiel - > EE_StabMembBiel ( i ) =
delta_energie + pt_StabMembBiel - > EE_StabMembBiel_t ( i ) ;
} ;
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
} ; // fin du cas if (limitation)
} ; // fin du cas ou on applique éventuellement une limitation
// 2) ---- on met à jour la raideur et résidu global
// en tenant compte de la stabilisation
( * residu ) + = ( * resD ) ;
2021-09-27 12:42:13 +02:00
2023-05-03 17:23:49 +02:00
# ifdef MISE_AU_POINT
if ( ( niveau_affichage > 9 ) )
{ cout < < " \n residu de stabilisation: " < < ( * resD ) ;
cout < < " \n fin debug ElemMeca::Cal_explicit_StabMembBiel(.. " ;
} ;
# endif
2021-09-27 12:42:13 +02:00
} // fin du cas où on utilise la normale aux noeuds
} ;
break ;
case LIGNE :
break ;
default : // on ne fait rien pour l'instant pour les autres types
break ;
} ;
} ;
} ;
// fonction a renseigner par les classes dérivées, concernant les répercutions
// éventuelles due à la suppression de tous les frontières
// nums_i : donnent les listes de frontières supprimées
void ElemMeca : : Prise_en_compte_des_consequences_suppression_tous_frontieres ( )
{ // il faut supprimer toutes les déformations liés aux frontières
int tail_S = defSurf . Taille ( ) ;
for ( int i = 1 ; i < = tail_S ; i + + )
{ delete defSurf ( i ) ;
defSurf ( i ) = NULL ;
} ;
int tail_A = defArete . Taille ( ) ;
for ( int i = 1 ; i < = tail_A ; i + + )
{ delete defArete ( i ) ;
defArete ( i ) = NULL ;
} ;
} ;
// idem pour une frontière (avant qu'elle soit supprimée)
void ElemMeca : : Prise_en_compte_des_consequences_suppression_une_frontiere ( ElFrontiere * elemFront )
{ int taille = tabb . Taille ( ) ;
// on choisit en fonction du type de frontière
if ( ( elemFront - > Type_geom_front ( ) = = LIGNE ) & & ( defArete . Taille ( ) ! = 0 ) )
{ int taille_ligne = taille - posi_tab_front_lin ; // a priori = cas sans les points
if ( posi_tab_front_point ! = 0 ) // cas où il y a des points
taille_ligne = posi_tab_front_point - posi_tab_front_lin ;
for ( int i = 1 ; i < = taille_ligne ; i + + )
{ if ( ( tabb ( i + posi_tab_front_lin ) = = elemFront ) & & ( defArete ( i ) ! = NULL ) )
{ delete defArete ( i ) ; defArete ( i ) = NULL ; break ; }
}
}
else if ( ( elemFront - > Type_geom_front ( ) = = SURFACE ) & & ( defSurf . Taille ( ) ! = 0 ) )
{ if ( posi_tab_front_lin ! = 0 ) // si == 0 cela signifie qu'il n'y a pas de surface à supprimer !!
{ for ( int i = 1 ; i < = posi_tab_front_lin ; i + + ) // posi_tab_front_lin == le nombre de surfaces, qui sont obligatoirement en début de tableau
if ( ( tabb ( i ) = = elemFront ) & & ( defSurf ( i ) ! = NULL ) )
{ delete defSurf ( i ) ; defSurf ( i ) = NULL ; break ; }
} ;
} ;
} ;
// Calcul des frontieres de l'element
// creation des elements frontieres et retour du tableau de ces elements
// la création n'a lieu qu'au premier appel
// ou lorsque l'on force le paramètre force a true
// dans ce dernier cas seul les frontière effacées sont recréée
// cas :
// = 0 -> on veut toutes les frontières
// = 1 -> on veut uniquement les surfaces
// = 2 -> on veut uniquement les lignes
// = 3 -> on veut uniquement les points
// = 4 -> on veut les surfaces + les lignes
// = 5 -> on veut les surfaces + les points
// = 6 -> on veut les lignes + les points
Tableau < ElFrontiere * > const & ElemMeca : : Frontiere_elemeca ( int cas , bool force )
{ // le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
int taille = tabb . Taille ( ) ; // la taille initiales des frontières
if ( force ) // dans ce cas on commence par tout effacer
{ // on efface les surfaces (s'il y en a)
for ( int i = 1 ; i < = posi_tab_front_lin ; i + + )
{ if ( tabb ( i ) ! = NULL )
{ delete tabb ( i ) ; //on commence par supprimer
tabb ( i ) = NULL ;
// on supprime également éventuellement la déformation associée
if ( defSurf . Taille ( ) ! = 0 )
if ( defSurf ( i ) ! = NULL )
{ delete defSurf ( i ) ; defSurf ( i ) = NULL ; } ;
ind_front_surf = 0 ; // on indique qu'il ne reste plus de frontière surface
} ;
} ;
// on efface les lignes (s'il y en a)
for ( int i = 1 ; i < = posi_tab_front_point - posi_tab_front_lin ; i + + )
{ if ( tabb ( i + posi_tab_front_lin ) ! = NULL )
{ delete tabb ( i + posi_tab_front_lin ) ; //on commence par supprimer
tabb ( i + posi_tab_front_lin ) = NULL ;
// on supprime également éventuellement la déformation associée
if ( defArete . Taille ( ) ! = 0 )
if ( defArete ( i ) ! = NULL )
{ delete defArete ( i ) ; defArete ( i ) = NULL ; } ;
ind_front_lin = 0 ; // on indique qu'il ne reste plus de frontière ligne
} ;
} ;
// on efface les points (s'il y en a)
for ( int i = 1 ; i < = taille - posi_tab_front_point ; i + + )
{ if ( tabb ( i + posi_tab_front_point ) ! = NULL )
{ delete tabb ( i + posi_tab_front_point ) ; //on commence par supprimer
tabb ( i + posi_tab_front_point ) = NULL ;
ind_front_point = 0 ; // on indique qu'il ne reste plus de frontière ligne
} ;
} ;
} ;
// -- maintenant on s'occupe de la construction conditionnelle
bool built_surf = false ; bool built_ligne = false ; bool built_point = false ;
switch ( cas )
{ case 0 : built_surf = built_ligne = built_point = true ; break ;
case 1 : built_surf = true ; break ;
case 2 : built_ligne = true ; break ;
case 3 : built_point = true ; break ;
case 4 : built_surf = built_ligne = true ; break ;
case 5 : built_surf = built_point = true ; break ;
case 6 : built_ligne = built_point = true ; break ;
} ;
if ( ( ( ind_front_surf = = 0 ) & & ( ind_front_lin = = 0 ) & & ( ind_front_point = = 0 ) ) | | force )
{
// récup de l'élément géométrique
ElemGeomC0 & el = ElementGeometrique ( ) ;
int tail_ar = el . NbSe ( ) ; // nombre potentiel d'arêtes
int tail_fa = el . NbFe ( ) ; // nombre potentiel de faces
int tail_po = el . Nbne ( ) ; // nombre potentiel de points
// récup du tableau de ddl actuel
const DdlElement & tdd = TableauDdl ( ) ;
// --- on va construire en fonction des indicateurs des tableaux intermédiaires
int new_posi_tab_front_point = 0 ; //init par défaut
int new_posi_tab_front_lin = 0 ; //init par défaut
int new_ind_front_point = 0 ;
int new_ind_front_lin = 0 ;
int new_ind_front_surf = 0 ;
// -- pour les surfaces
Tableau < ElFrontiere * > tabb_surf ;
if ( ( built_surf ) & & ( ( ind_front_surf = = 0 ) | | force ) )
{ tabb_surf . Change_taille ( tail_fa , NULL ) ; // init par défaut
for ( int num = 1 ; num < = tail_fa ; num + + )
{ int nbnoe = el . Nonf ( ) ( num ) . Taille ( ) ; // nb noeud de la surface
Tableau < Noeud * > tab ( nbnoe ) ; // les noeuds de la frontiere
DdlElement ddelem ( nbnoe ) ; // les ddlelements des noeuds frontieres
for ( int i = 1 ; i < = nbnoe ; i + + )
{ tab ( i ) = tab_noeud ( el . Nonf ( ) ( num ) ( i ) ) ;
ddelem . Change_un_ddlNoeudElement ( i , tdd ( el . Nonf ( ) ( num ) ( i ) ) ) ;
} ;
tabb_surf ( num ) = new_frontiere_surf ( num , tab , ddelem ) ;
} ;
// nouveau indicateur d'existence
new_ind_front_surf = 1 ;
// on positionne les nouvelles positions
new_posi_tab_front_point + = tail_fa ;
new_posi_tab_front_lin + = tail_fa ;
} ;
// -- pour les lignes
Tableau < ElFrontiere * > tabb_ligne ;
if ( ( built_ligne ) & & ( ( ind_front_lin = = 0 ) | | force ) )
{ tabb_ligne . Change_taille ( tail_ar , NULL ) ; // init par défaut
for ( int num = 1 ; num < = tail_ar ; num + + )
{ int nbnoe = el . NonS ( ) ( num ) . Taille ( ) ; // nb noeud de l'arête
Tableau < Noeud * > tab ( nbnoe ) ; // les noeuds de l'arête frontiere
DdlElement ddelem ( nbnoe ) ; // les ddlelements des noeuds frontieres
for ( int i = 1 ; i < = nbnoe ; i + + )
{ tab ( i ) = tab_noeud ( el . NonS ( ) ( num ) ( i ) ) ;
ddelem . Change_un_ddlNoeudElement ( i , tdd ( el . NonS ( ) ( num ) ( i ) ) ) ;
} ;
tabb_ligne ( num ) = new_frontiere_lin ( num , tab , ddelem ) ;
} ;
// nouveau indicateur d'existence
new_ind_front_lin = 1 ;
// on positionne les nouvelles positions
new_posi_tab_front_point + = tail_ar ;
} ;
// -- pour les points
Tableau < ElFrontiere * > tabb_point ;
if ( ( built_point ) & & ( ( ind_front_point = = 0 ) | | force ) )
{ tabb_point . Change_taille ( tail_po , NULL ) ; // init par défaut
// maintenant création des frontière point éventuellement
Tableau < Noeud * > tab ( 1 ) ; // les noeuds de la frontiere (tab de travail)
DdlElement ddelem ( 1 ) ; // les ddlelements des points frontieres (tab de travail)
for ( int i = 1 ; i < = tail_po ; i + + )
if ( tabb_point ( i ) = = NULL )
{ tab ( 1 ) = tab_noeud ( i ) ;
ddelem . Change_un_ddlNoeudElement ( 1 , tdd ( i ) ) ;
tabb_point ( i ) = new FrontPointF ( tab , ddelem ) ;
} ;
// nouveau indicateur d'existence
new_ind_front_point = 1 ;
} ;
// --- mise à jour du tableau globale et des indicateurs ad hoc
int taille_finale = tabb_surf . Taille ( ) + tabb_ligne . Taille ( ) + tabb_point . Taille ( ) ;
tabb . Change_taille ( taille_finale ) ;
// cas des points
if ( new_ind_front_point ) // là is s'agit de nouveaux éléments
{ for ( int i = tail_po ; i > 0 ; i - - ) // transfert pour les noeuds
{ tabb ( i + new_posi_tab_front_point ) = tabb_point ( i ) ; }
}
else if ( ind_front_point ) // là il s'agit d'anciens éléments
{ for ( int i = tail_po ; i > 0 ; i - - ) // transfert pour les noeuds en descendant
{ tabb ( i + new_posi_tab_front_point ) = tabb ( i + posi_tab_front_point ) ; }
} ;
// cas des lignes
if ( new_ind_front_lin ) // là il s'agit de nouveaux éléments
{ for ( int i = 1 ; i < = tail_ar ; i + + ) // transfert
{ tabb ( i + new_posi_tab_front_lin ) = tabb_ligne ( i ) ; }
}
else if ( ind_front_lin ) // là il s'agit d'anciens éléments
{ for ( int i = tail_ar ; i > 0 ; i - - ) // transfert en descendant
{ tabb ( i + new_posi_tab_front_lin ) = tabb ( i + posi_tab_front_lin ) ; }
} ;
// cas des surfaces
if ( new_ind_front_surf ) // là is s'agit de nouveaux éléments
{ for ( int i = 1 ; i < = tail_fa ; i + + ) // transfert
{ tabb ( i ) = tabb_surf ( i ) ; }
} ;
// dans le cas où il y avait des anciens éléments, il n'y a rien n'a faire
// car le redimentionnement de tabb ne change pas les premiers éléments
// mis à jour des indicateurs
ind_front_surf = new_ind_front_surf ;
posi_tab_front_lin = new_posi_tab_front_lin ;
ind_front_lin = new_ind_front_lin ;
posi_tab_front_point = new_posi_tab_front_point ;
ind_front_point = new_ind_front_point ;
} ;
// retour du tableau
return ( Tableau < ElFrontiere * > & ) tabb ;
} ;
// ramène la frontière point
// éventuellement création des frontieres points de l'element et stockage dans l'element
// si c'est la première fois sinon il y a seulement retour de l'elements
// a moins que le paramètre force est mis a true
// dans ce dernier cas la frontière effacéee est recréée
// num indique le numéro du point à créer (numérotation EF)
ElFrontiere * const ElemMeca : : Frontiere_points ( int num , bool force )
2023-05-03 17:23:49 +02:00
{ // -- tout d'abord on évacue le cas où il n'y a pas de frontière point à calculer
2021-09-27 12:42:13 +02:00
// récup de l'élément géométrique
ElemGeomC0 & el = ElementGeometrique ( ) ;
int tail_po = el . Nbne ( ) ; // nombre potentiel de points
if ( num > tail_po )
return NULL ;
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
// on regarde si les frontières points existent sinon on les crée
if ( ind_front_point = = 1 )
{ return ( ElFrontiere * ) tabb ( posi_tab_front_point + num ) ; }
else if ( ind_front_point = = 2 )
// cas où certaines frontières existent
{ if ( tabb ( posi_tab_front_point + num ) ! = NULL )
return ( ElFrontiere * ) tabb ( posi_tab_front_point + num ) ;
} ;
// arrivée ici cela veut dire que la frontière point n'existe pas
// on l'a reconstruit éventuellement
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
if ( ( ind_front_point = = 0 ) | | force )
{ // récup du tableau de ddl
const DdlElement & tdd = TableauDdl ( ) ;
int taille = tabb . Taille ( ) ; // la taille initiales des frontières
int tail_fa = el . NbFe ( ) ; // nombre potentiel de faces
int tail_ar = el . NbSe ( ) ; // nombre potentiel d'arêtes
// dimensionnement du tableau de frontières ligne si nécessaire
if ( ind_front_point = = 0 )
{ if ( ( ind_front_lin > 0 ) & & ( ind_front_surf = = 0 ) )
// cas où les frontières lignes existent seules, on ajoute les points
{ int taille_finale = tail_ar + tail_po ;
tabb . Change_taille ( taille_finale ) ;
for ( int i = 1 ; i < = tail_ar ; i + + ) // transfert pour les lignes
tabb ( i ) = tabb ( i + tail_ar ) ;
posi_tab_front_point = tail_ar ;
posi_tab_front_lin = 0 ; // car on n'a pas de surface
}
else if ( ( ind_front_lin > 0 ) & & ( ind_front_surf > 0 ) )
// cas où les frontières lignes existent et surfaces et pas de points, donc on les rajoutes
{ int taille_finale = tail_fa + tail_po + tail_ar ;
tabb . Change_taille ( taille_finale ) ; // les grandeurs pour les surfaces et les lignes sont
// conservées, donc pas de transferts à prévoir
posi_tab_front_point = tail_ar + tail_fa ; // après les faces et les lignes
posi_tab_front_lin = tail_fa ; // après les surfaces
}
else
{ // cas où il n'y a pas de frontières lignes
if ( ind_front_surf = = 0 ) // cas où il n'y a pas de surface
{ tabb . Change_taille ( tail_po , NULL ) ; // on n'a pas de ligne, pas de point et pas de surface
posi_tab_front_point = posi_tab_front_lin = 0 ; }
else { tabb . Change_taille ( tail_po + tail_fa ) ; // cas où les surfaces existent
// le redimensionnement n'affecte pas les surfaces qui sont en début de tableau
posi_tab_front_lin = posi_tab_front_point = tail_fa ; // après les surfaces
} ;
} ;
// et on met les pointeurs de points en NULL
for ( int i = 1 ; i < = tail_po ; i + + )
{ tabb ( i + posi_tab_front_point ) = NULL ; }
} ;
// maintenant création de la frontière point
Tableau < Noeud * > tab ( 1 ) ; // les noeuds de la frontiere
tab ( 1 ) = tab_noeud ( num ) ;
DdlElement ddelem ( 1 ) ; // les ddlelements des points frontieres
ddelem . Change_un_ddlNoeudElement ( 1 , tdd ( num ) ) ;
tabb ( posi_tab_front_point + num ) = new FrontPointF ( tab , ddelem ) ;
// on met à jour l'indicateur ind_front_point
ind_front_point = 1 ; // a priori
for ( int npoint = 1 ; npoint < = tail_po ; npoint + + )
if ( tabb ( posi_tab_front_point + npoint ) = = NULL )
{ ind_front_point = 2 ; break ; } ;
} ;
// maintenant normalement la frontière est créé on la ramène
return ( ElFrontiere * ) tabb ( posi_tab_front_point + num ) ;
} ;
// ramène la frontière linéique
// éventuellement création des frontieres linéique de l'element et stockage dans l'element
// si c'est la première fois et en 3D sinon il y a seulement retour de l'elements
// a moins que le paramètre force est mis a true
// dans ce dernier cas la frontière effacéee est recréée
// num indique le numéro de l'arête à créer (numérotation EF)
// nbneA: nombre de noeuds des segments frontières
// el : l'élément
ElFrontiere * const ElemMeca : : Frontiere_lineique ( int num , bool force )
{ // -- tout d'abord on évacue le cas où il n'y a pas de frontière linéique à calculer
// récup de l'élément géométrique
ElemGeomC0 & el = ElementGeometrique ( ) ;
int tail_ar = el . NbSe ( ) ; // nombre potentiel d'arêtes
if ( num > tail_ar )
return NULL ;
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
// on regarde si les frontières linéiques existent sinon on les crée
if ( ind_front_lin = = 1 )
{ return ( ElFrontiere * ) tabb ( posi_tab_front_lin + num ) ; }
else if ( ind_front_lin = = 2 )
// cas où certaines frontières existent
{ if ( tabb ( posi_tab_front_lin + num ) ! = NULL )
return ( ElFrontiere * ) tabb ( posi_tab_front_lin + num ) ;
} ;
// arrivée ici cela veut dire que la frontière ligne n'existe pas
// on l'a reconstruit
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
if ( ( ind_front_lin = = 0 ) | | force )
{ // récup du tableau de ddl
const DdlElement & tdd = TableauDdl ( ) ;
int taille = tabb . Taille ( ) ; // la taille initiales des frontières
int tail_fa = el . NbFe ( ) ; // nombre potentiel de faces
int tail_po = el . Nbne ( ) ; // nombre potentiel de points
// dimensionnement du tableau de frontières ligne si nécessaire
if ( ind_front_lin = = 0 )
2023-05-03 17:23:49 +02:00
{ if ( ( ind_front_point > 0 ) & & ( ind_front_surf = = 0 ) )
// cas où les frontières points existent seules, on ajoute les lignes
{ int taille_finale = tail_ar + tail_po ;
tabb . Change_taille ( taille_finale ) ;
for ( int i = 1 ; i < = tail_po ; i + + )
tabb ( i + tail_ar ) = tabb ( i ) ;
posi_tab_front_point = tail_ar ;
posi_tab_front_lin = 0 ; // car on n'a pas de surface
}
else if ( ( ind_front_point > 0 ) & & ( ind_front_surf > 0 ) )
// cas où les frontières points existent et surfaces et pas de ligne, donc on les rajoutes
{ int taille_finale = tail_fa + tail_po + tail_ar ;
tabb . Change_taille ( taille_finale ) ; // les grandeurs pour les surfaces sont conservées
for ( int i = 1 ; i < = tail_po ; i + + ) // transfert pour les noeuds
{ tabb ( i + tail_ar + tail_fa ) = tabb ( i + tail_fa ) ; } ;
posi_tab_front_point = tail_ar + tail_fa ; // après les faces et les lignes
posi_tab_front_lin = tail_fa ; // après les surfaces
}
else
{ // cas où il n'y a pas de frontières points
if ( ind_front_surf = = 0 ) // cas où il n'y a pas de surface
{ tabb . Change_taille ( tail_ar , NULL ) ; // on n'a pas de ligne, pas de point et pas de surface
posi_tab_front_lin = 0 ; posi_tab_front_point = tail_ar ;
}
else { tabb . Change_taille ( tail_ar + tail_fa ) ; // cas où les surfaces existent
// le redimensionnement n'affecte pas les surfaces qui sont en début de tableau
posi_tab_front_lin = tail_fa ; // après les surfaces
posi_tab_front_point = tail_fa + tail_ar ;
} ;
} ;
// et on met les pointeurs de lignes en NULL
for ( int i = 1 ; i < = tail_ar ; i + + ) // transfert pour les noeuds
{ tabb ( i + posi_tab_front_lin ) = NULL ; }
} ;
2021-09-27 12:42:13 +02:00
// maintenant création de la ligne
int nbnoe = el . NonS ( ) ( num ) . Taille ( ) ; // nb noeud de l'arête
Tableau < Noeud * > tab ( nbnoe ) ; // les noeuds de l'arête frontiere
DdlElement ddelem ( nbnoe ) ; // les ddlelements des noeuds frontieres
for ( int i = 1 ; i < = nbnoe ; i + + )
{ tab ( i ) = tab_noeud ( el . NonS ( ) ( num ) ( i ) ) ;
ddelem . Change_un_ddlNoeudElement ( i , tdd ( el . NonS ( ) ( num ) ( i ) ) ) ;
} ;
tabb ( posi_tab_front_lin + num ) = new_frontiere_lin ( num , tab , ddelem ) ;
ind_front_lin = 1 ; // a priori
for ( int nlign = 1 ; nlign < = tail_ar ; nlign + + )
if ( tabb ( posi_tab_front_lin + nlign ) = = NULL )
{ ind_front_lin = 2 ; break ; } ;
} ;
// maintenant normalement la frontière est créé on la ramène
return ( ElFrontiere * ) tabb ( posi_tab_front_lin + num ) ;
} ;
// ramène la frontière surfacique
// éventuellement création des frontieres surfacique de l'element et stockage dans l'element
// si c'est la première fois sinon il y a seulement retour de l'elements
// a moins que le paramètre force est mis a true
// dans ce dernier cas la frontière effacéee est recréée
// num indique le numéro de la surface à créer (numérotation EF)
ElFrontiere * const ElemMeca : : Frontiere_surfacique ( int num , bool force )
{ // -- tout d'abord on évacue le cas où il n'y a pas de frontière surfacique à calculer
// récup de l'élément géométrique
ElemGeomC0 & el = ElementGeometrique ( ) ;
int tail_fa = 0 ; // init
// test --- retour rapide si l'existance n'est pas possible
if ( ParaGlob : : AxiSymetrie ( ) )
// dans le cas axiSymétrique, le test est un peu différent car les surfaces sont générées
// par les lignes,
{ tail_fa = el . NbSe ( ) ; // nombre potentiel de segments
if ( num > tail_fa )
return NULL ;
}
else // dans le cas où on n'est pas en axi, test simple
{ tail_fa = el . NbFe ( ) ; // nombre potentiel de faces
if ( num > tail_fa )
return NULL ;
} ;
// --- la suite concerne le cas où on peut faire des surfaces (non axi ou axi !)
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
// on regarde si les frontières surfacique existent sinon on les crée
if ( ind_front_surf = = 1 )
{ return ( ElFrontiere * ) tabb ( num ) ; }
else if ( ind_front_surf = = 2 )
// cas où certaines frontières existent
{ if ( tabb ( num ) ! = NULL )
return ( ElFrontiere * ) tabb ( num ) ;
} ;
// arrivée ici cela veut dire que la frontière surface n'existe pas
// on l'a reconstruit
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
if ( ( ind_front_surf = = 0 ) | | force )
{ // récup du tableau de ddl
const DdlElement & tdd = TableauDdl ( ) ;
int taille = tabb . Taille ( ) ; // la taille initiales des frontières
int tail_ar = el . NbSe ( ) ; // nombre potentiel d'arêtes
int tail_po = el . Nbne ( ) ; // nombre potentiel de points
// dimensionnement du tableau de frontières surfaces si nécessaire
if ( ind_front_surf = = 0 )
{ if ( ( ind_front_point > 0 ) & & ( ind_front_lin = = 0 ) )
// cas où les frontières points existent seules, on ajoute les surfaces
{ int taille_finale = tail_fa + tail_po ;
tabb . Change_taille ( taille_finale ) ;
for ( int i = 1 ; i < = tail_po ; i + + )
tabb ( i + tail_fa ) = tabb ( i ) ;
posi_tab_front_lin = posi_tab_front_point = tail_fa ;
}
else if ( ( ind_front_point > 0 ) & & ( ind_front_lin > 0 ) )
// cas où les frontières points existent et lignes et pas de surfaces, donc on les rajoutes
{ int taille_finale = tail_fa + tail_po + tail_ar ;
tabb . Change_taille ( taille_finale ) ;
// --transfert pour les noeuds et les lignes
for ( int i = 1 ; i < = tail_po ; i + + ) // transfert pour les noeuds
{ tabb ( i + tail_ar + tail_fa ) = tabb ( i + tail_ar ) ; } ;
for ( int i = 1 ; i < = tail_ar ; i + + ) // transfert pour les lignes
{ tabb ( i + tail_fa ) = tabb ( i ) ; } ;
// --def des indicateurs
posi_tab_front_point = tail_ar + tail_fa ; // après les faces et les lignes
posi_tab_front_lin = tail_fa ; // après les surfaces
}
else
{ // cas où il n'y a pas de frontières points
if ( ind_front_lin = = 0 ) // cas où il n'y a pas de lignes
{ tabb . Change_taille ( tail_fa , NULL ) ; // on n'a pas de ligne, pas de point et pas de surface
posi_tab_front_lin = posi_tab_front_point = tail_fa ;
} // on peut tout mettre à NULL
else { tabb . Change_taille ( tail_ar + tail_fa ) ; // cas où les lignes existent
for ( int i = 1 ; i < = tail_ar ; i + + ) // transfert pour les lignes
tabb ( i + tail_fa ) = tabb ( i ) ;
posi_tab_front_lin = tail_fa ; // après les surfaces
posi_tab_front_point = tail_fa + tail_ar ;
} ;
} ;
// --et on met les pointeurs de surfaces en NULL
for ( int i = 1 ; i < = tail_fa ; i + + )
tabb ( i ) = NULL ;
} ;
// maintenant création de la surface
int nbnoe = el . Nonf ( ) ( num ) . Taille ( ) ; // b noeud de la surface
Tableau < Noeud * > tab ( nbnoe ) ; // les noeuds de la frontiere
DdlElement ddelem ( nbnoe ) ; // les ddlelements des noeuds frontieres
for ( int i = 1 ; i < = nbnoe ; i + + )
{ tab ( i ) = tab_noeud ( el . Nonf ( ) ( num ) ( i ) ) ;
ddelem . Change_un_ddlNoeudElement ( i , tdd ( el . Nonf ( ) ( num ) ( i ) ) ) ;
} ;
tabb ( num ) = new_frontiere_surf ( num , tab , ddelem ) ;
ind_front_surf = 1 ; // a priori
for ( int nsurf = 1 ; nsurf < = tail_fa ; nsurf + + )
if ( tabb ( nsurf ) = = NULL )
{ ind_front_surf = 2 ; break ; } ;
} ;
// maintenant normalement la frontière est créé on la ramène
return ( ElFrontiere * ) tabb ( num ) ;
} ;
// calcul des intégrales de volume et de volume + temps
// cas = 1 : pour un appel après un calcul implicit
// cas = 2 : pour un appel après un calcul explicit
void ElemMeca : : Calcul_Integ_vol_et_temps ( int cas , const double & poid_jacobien , int iteg )
{ // 1) intégration de volume uniquement
//List_io <TypeQuelconque> integ_vol_typeQuel,integ_vol_typeQuel_t;
if ( integ_vol_typeQuel ! = NULL )
// on balaie les conteneurs
{ int taille = integ_vol_typeQuel - > Taille ( ) ;
Tableau < TypeQuelconque > & tab_integ = * integ_vol_typeQuel ; // pour simplifier
for ( int il = 1 ; il < = taille ; il + + )
// l'intégrale n'est effectuée que si l'index est positif, sinon, il s'agit
// d'une valeur figée
if ( ( * index_Integ_vol_typeQuel ) ( il ) )
{ // différent choix
switch ( tab_integ ( il ) . Const_Grandeur_pointee ( ) - > Type_enumGrandeurParticuliere ( ) )
{ case PARTICULIER_DDL_ETENDU :
{ Grandeur_Ddl_etendu & gr
= * ( ( Grandeur_Ddl_etendu * ) tab_integ ( il ) . Grandeur_pointee ( ) ) ; // pour simplifier
Ddl_etendu & ddl = * ( gr . ConteneurDdl_etendu ( ) ) ; // récup du ddl
// on va utiliser la méthode Valeur_multi ce qui n'est pas optimal en vitesse
// mais 1) c'est cohérent avec l'implantation existante
// 2) identique avec les ddl déjà accessibles
List_io < Ddl_enum_etendu > inter ; // une liste intermédiaire
inter . push_back ( ddl . DdlEnumEtendu ( ) ) ;
bool absolue = true ; // on se place systématiquement en absolu
// calcul de la valeur et retour dans tab_d
Tableau < double > tab_d ( ElemMeca : : Valeur_multi ( absolue , TEMPS_tdt , inter , iteg , - cas ) ) ;
// on cumule l'intégrale
ddl . Valeur ( ) + = tab_d ( 1 ) * poid_jacobien ;
////---- debug
//cout << "\n debug: ElemMeca::Calcul_Integ_vol_et_temps(";
//cout << " tab_d(1)= "<< tab_d(1) << " ";ddl.Affiche();
//tab_noeud(1)->Affiche();
////---- fin debug
break ;
}
case PARTICULIER_VECTEUR_NOMMER :
{ Grandeur_Vecteur_Nommer & gr
= * ( ( Grandeur_Vecteur_Nommer * ) tab_integ ( il ) . Grandeur_pointee ( ) ) ; // pour simplifier
Fonction_nD * fct = gr . Fct ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi ( absolue , TEMPS_tdt , li_enu_scal , iteg , - cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , iteg , - cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// on cumule l'intégrale
Vecteur & V = gr . ConteneurVecteur ( ) ;
int taille = V . Taille ( ) ;
for ( int i = 1 ; i < = taille ; i + + )
V ( i ) + = tab_ret ( i ) * poid_jacobien ;
////----- debug
//cout << "\n debug *** ElemMeca::Grandeur_particuliere( ";
//cout << "\n li_enu_scal= "<< li_enu_scal << ", val_ddl_enum= " << val_ddl_enum;
//cout << "\n li_quelc= "<< li_quelc;
//V.Affiche();cout << " tab_ret= " << tab_ret << endl;;
////----- fin debug
break ;
}
break ;
default :
cout < < " \n *** pas d'integration prevu pour la grandeur " < < tab_integ ( il ) ;
if ( ParaGlob : : NiveauImpression ( ) > 0 )
cout < < " \n ElemMeca::Calcul_Integ_vol_et_temps(... " < < endl ;
} ;
} ;
} ;
// 2) intégration de volume et en temps
//List_io <TypeQuelconque> integ_vol_t_typeQuel,integ_vol_t_typeQuel_t;
if ( integ_vol_t_typeQuel ! = NULL )
// on balaie les conteneurs
{ int taille = integ_vol_t_typeQuel - > Taille ( ) ;
Tableau < TypeQuelconque > & tab_integ = * integ_vol_t_typeQuel ; // pour simplifier
// recup de l'incrément de temps
double deltat = ParaGlob : : Variables_de_temps ( ) . IncreTempsCourant ( ) ;
for ( int il = 1 ; il < = taille ; il + + )
// l'intégrale n'est effectuée que si l'index est positif, sinon, il s'agit
// d'une valeur figée
if ( ( * index_Integ_vol_t_typeQuel ) ( il ) )
{ // différent choix
switch ( tab_integ ( il ) . Const_Grandeur_pointee ( ) - > Type_enumGrandeurParticuliere ( ) )
{ case PARTICULIER_DDL_ETENDU :
{ Grandeur_Ddl_etendu & gr
= * ( ( Grandeur_Ddl_etendu * ) tab_integ ( il ) . Grandeur_pointee ( ) ) ; // pour simplifier
Ddl_etendu & ddl = * ( gr . ConteneurDdl_etendu ( ) ) ; // récup du ddl
// on va utiliser la méthode Valeur_multi ce qui n'est pas optimal en vitesse
// mais 1) c'est cohérent avec l'implantation existante
// 2) identique avec les ddl déjà accessibles
List_io < Ddl_enum_etendu > inter ; // une liste intermédiaire
inter . push_back ( ddl . DdlEnumEtendu ( ) ) ;
bool absolue = true ; // on se place systématiquement en absolu
// calcul de la valeur et retour dans tab_d
Tableau < double > tab_d ( ElemMeca : : Valeur_multi ( absolue , TEMPS_tdt , inter , iteg , - cas ) ) ;
// on cumule l'intégrale
ddl . Valeur ( ) + = deltat * tab_d ( 1 ) * poid_jacobien ;
break ;
}
case PARTICULIER_VECTEUR_NOMMER :
{ Grandeur_Vecteur_Nommer & gr
= * ( ( Grandeur_Vecteur_Nommer * ) tab_integ ( il ) . Grandeur_pointee ( ) ) ; // pour simplifier
Fonction_nD * fct = gr . Fct ( ) ; // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi ( absolue , TEMPS_tdt , li_enu_scal , iteg , - cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_tdt , li_quelc , iteg , - cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// on cumule l'intégrale
Vecteur & V = gr . ConteneurVecteur ( ) ;
int taille = V . Taille ( ) ;
for ( int i = 1 ; i < = taille ; i + + )
V ( i ) + = deltat * tab_ret ( i ) * poid_jacobien ;
break ;
}
break ;
default :
cout < < " \n *** pas d'integration prevu pour la grandeur " < < tab_integ ( il ) ;
if ( ParaGlob : : NiveauImpression ( ) > 0 )
cout < < " \n ElemMeca::Calcul_Integ_vol_et_temps(... " < < endl ;
} ;
} ;
} ;
} ;
// au premier pti init éventuel des intégrales de volume et temps
void ElemMeca : : Init_Integ_vol_et_temps ( )
{ // 1) intégration de volume uniquement
if ( integ_vol_typeQuel ! = NULL )
{ int taille = integ_vol_typeQuel - > Taille ( ) ;
Tableau < TypeQuelconque > & tab_integ = * integ_vol_typeQuel ; // pour simplifier
for ( int il = 1 ; il < = taille ; il + + )
// l'intégrale n'est effectuée que si l'index est positif, sinon, il s'agit
// d'une valeur figée
if ( ( * index_Integ_vol_typeQuel ) ( il ) )
{ tab_integ ( il ) . Grandeur_pointee ( ) - > InitParDefaut ( ) ; // on initialise le conteneur
} ;
} ;
// 2) intégration de volume et en temps
if ( integ_vol_t_typeQuel ! = NULL )
{ int taille = integ_vol_t_typeQuel - > Taille ( ) ;
Tableau < TypeQuelconque > & tab_integ = * integ_vol_t_typeQuel ; // pour simplifier
Tableau < TypeQuelconque > & tab_integ_t = * integ_vol_t_typeQuel_t ; // ""
for ( int il = 1 ; il < = taille ; il + + )
// l'intégrale n'est effectuée que si l'index est positif, sinon, il s'agit
// d'une valeur figée
if ( ( * index_Integ_vol_t_typeQuel ) ( il ) )
{ ( * tab_integ ( il ) . Grandeur_pointee ( ) ) = ( * tab_integ_t ( il ) . Grandeur_pointee ( ) ) ;
} ;
} ;
} ;
// définition d'un repère d'anisotropie à un point d'intégration
BaseH ElemMeca : : DefRepereAnisotropie
( int iteg , LesFonctions_nD * lesFonctionsnD , const BlocGen & bloc )
{ // la définition du repère dépend du type de définition
// (cf. la def du bloc dans LesMaillages::Completer .... dans LesMaillages.cc)
if ( bloc . Nom ( 4 ) = = " par_fonction_nD_ " )
{ // on doit simplement appeler la fonction nD
// Nom(5) est le nom de la fonction nD
Fonction_nD * fct = lesFonctionsnD - > Trouve ( bloc . Nom ( 5 ) ) ;
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io < Ddl_enum_etendu > & li_enu_scal = fct - > Li_enu_etendu_scalaire ( ) ;
List_io < TypeQuelconque > & li_quelc = fct - > Li_equi_Quel_evolue ( ) ;
bool absolue = true ; // on se place systématiquement en absolu
const Met_abstraite : : Impli * ex_impli = NULL ;
const Met_abstraite : : Expli_t_tdt * ex_expli_tdt = NULL ;
const Met_abstraite : : Expli * ex_expli = NULL ;
// pour l'instant on ne peut s'occuper que des grandeurs à t au maximum: mais pas
// à tdt
bool premier_calcul = true ;
// def->ChangeNumInteg(iteg);
// const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);
// ex_expli = &ex;
// // on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
// // pour les grandeurs strictement scalaire
// Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
// (absolue,TEMPS_t,li_enu_scal,*def,ex_impli,ex_expli_tdt,ex_expli)
// );
// // on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
// // pour les Coordonnees et Tenseur
// Valeurs_Tensorielles_interpoler_ou_calculer
// (absolue,TEMPS_t,li_quelc,*def,ex_impli,ex_expli_tdt,ex_expli);
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = 1 ; // on se met dans un cas normal avec utilisation de certaine partie de la métrique implicite
// a priori ce serait idem avec l'explicite
Tableau < double > val_ddl_enum ( ElemMeca : : Valeur_multi ( absolue , TEMPS_t , li_enu_scal , iteg , cas ) ) ;
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles ( absolue , TEMPS_t , li_quelc , iteg , cas ) ;
// calcul de la valeur et retour dans tab_ret
Tableau < double > & tab_ret = fct - > Valeur_FnD_Evoluee ( & val_ddl_enum , & li_enu_scal , & li_quelc , NULL , NULL ) ;
// on doit récupérer 6 composantes en 3D (pour le reste c'est actuellement en attente)
int dim = ParaGlob : : Dimension ( ) ;
if ( dim ! = 3 )
{ cout < < " \n arret: le type de construction de repere: par_fonction_nD_ "
< < " n'est implante que pour la dimension 3 pour l'instant !!! "
< < " affaire a suivre ....(et se plaindre! ) " ;
Sortie ( 1 ) ;
} ;
int nb_composante = tab_ret . Taille ( ) ;
if ( nb_composante ! = 6 )
{ cout < < " \n erreur dans la construction d'un repere d'orthotropie par_fonction_nD_ "
< < " la fonction nD ramene " < < nb_composante < < " composante(s) "
< < " or il en faut 6 uniquement ! " ;
Sortie ( 1 ) ;
} ;
// arrivée ici, cela veut dire que tout est ok
// --- construction de la base dans I_a -------
BaseH baseH ( 3 , 3 ) ;
// premier vecteur:
CoordonneeH & c1 = baseH . CoordoH ( 1 ) ;
for ( int i = 1 ; i < 4 ; i + + )
c1 ( i ) = tab_ret ( i ) ;
// on norme si c'est demandé
if ( bloc . Nom ( 3 ) = = " orthonorme_ " )
{ c1 . Normer ( ) ;
}
else
{ cout < < " \n arret: le type de repere: " < < bloc . Nom ( 3 )
< < " n'est pas encore implante !!! affaire a suivre ....(et se plaindre! ) " ;
Sortie ( 1 ) ;
} ;
// le deuxième vecteur:
CoordonneeH & c2 = baseH . CoordoH ( 2 ) ;
for ( int i = 1 ; i < 4 ; i + + )
c2 ( i ) = tab_ret ( i + 3 ) ;
// on norme
c2 . Normer ( ) ;
// le troisième vecteur
CoordonneeH & c3 = baseH . CoordoH ( 3 ) ;
c3 = Util : : ProdVec_coorH ( c1 , c2 ) ;
/* ***** la suite est supprimée: elle correspond au calcul des coordonnées locales alpha_a^{.i}
du repère d ' anisotropie . Le pb est que le repère local g_i à t = 0 peut changer pendant le calcul
par exemple dans le cas de l ' utilisation des contraintes planes on va avoir un repère de travail
qui change . Sans doute que ce sera un pb identique avec les Umat
donc on passe en paramètre le repère en coordonnée absolue et au moment de l ' utilisation
les coordonnées locales seront calculées
si c ' est ok par la suite , il faudra supprimer la partie commentée qui suit
// --- maintenant on s'occupe de calculer la base en local -----
// en fait il s'agit de définir les coordonnées locales par projection de la base
// globale en local
def - > ChangeNumInteg ( iteg ) ;
const Met_abstraite : : InfoExp_t & ex_infoExp_t = def - > RemontExp_t ( ) ;
int nb_vec = ex_infoExp_t . giB_0 - > NbVecteur ( ) ; // nb vecteur local
BaseH base_ret_H ( 3 , 3 ) ;
// les coordonnées locales : théta^te = V * g^te
// soit O_a la base d'anisotropie: O_a = beta_a^{.i} g_i
// ona beta_a^{.i} = O_a . g^i
// et on stocke dans base_ret_H de la manière suivante
// base_ret_H(a)(i) = beta_a^{.i}
for ( int a = 1 ; a < = 3 ; a + + )
{ CoordonneeH & inter = base_ret_H . CoordoH ( a ) ;
for ( int i = 1 ; i < nb_vec ; i + + )
inter ( i ) = baseB . CoordoB ( a ) * ( * ex_infoExp_t . giH_0 ) ( i ) ;
} ;
// dans le cas où nb_vec n'est pas égale à 3 il faut complèter la base
switch ( nb_vec )
{ case 3 : break ; // cas où c'est déjà ok
case 2 : // cas où l'élément est une surface, le troisième vecteur est la normale
{ CoordonneeH N_H = Util : : ProdVec_coorH ( ( * ex_infoExp_t . giH_0 ) ( 1 ) , ( * ex_infoExp_t . giH_0 ) ( 2 ) ) ;
for ( int a = 1 ; a < = 3 ; a + + )
base_ret_H . CoordoH ( a ) ( 3 ) = baseB . CoordoB ( a ) * N_H ;
break ;
}
case 1 : // cas où l'élément est un segment, il faut définir deux autres vecteurs
{ // on génère
la suite est problèmatique car lorsque l ' on changera de repère local , les coordonnées locales
devraient changer donc on arrête
CoordonneeH N_H = Util : : ProdVec_coorH ( ( * ex_infoExp_t . giH_0 ) ( 1 ) , ( * ex_infoExp_t . giH_0 ) ( 2 ) ) ;
for ( int a = 1 ; a < = 3 ; a + + )
base_ret_H . CoordoH ( a ) ( 3 ) = baseB . CoordoB ( a ) * N_H ;
break ;
}
default :
break ;
}
for ( int i = 1 ; i < 4 ; i + + )
{ CoordonneeH & inter = base_ret_H . CoordoH ( i ) ;
for ( int te = 1 ; te < = nb_vec ; te + + )
inter ( te ) = baseB . CoordoB ( i ) * ( * ex_infoExp_t . giH_0 ) ( te ) ;
} ;
return base_ret_H ;
*/
//cout << "\n debug ElemMeca::DefRepereAnisotropie ";
//for (int a = 1;a < 4; a++)
// {cout << "\n baseH("<<a<<"): ";
// baseH(a).Affiche();
// };
return baseH ;
}
else if ( bloc . Nom ( 4 ) = = " Par_projection_et_normale_au_plan_ " )
{ cout < < " \n arret: le type de construction de repere: " < < bloc . Nom ( 3 )
< < " n'est pas encore implante !!! affaire a suivre ....(et se plaindre! ) " ;
Sortie ( 1 ) ;
}
else
{ cout < < " \n **** erreur: type de construction de repere: " < < bloc . Nom ( 3 )
< < " inconnu ?? " ;
Sortie ( 1 ) ;
} ;
return BaseH ( ) ; // pour faire taire le compilo
} ;
// calcul éventuel de la normale à un noeud
// ce calcul existe pour les éléments 2D, 1D axi, et aussi pour les éléments 1D
// qui possède un repère d'orientation
// en retour coor = la normale si coor.Dimension() est = à la dimension de l'espace
// si le calcul n'existe pas --> coor.Dimension() = 0
// ramène un entier :
// == 1 : calcul normal
// == 0 : problème de calcul -> coor.Dimension() = 0
// == 2 : indique que le calcul n'est pas licite pour le noeud passé en paramètre
// c'est le cas par exemple des noeuds exterieurs pour les éléments SFE
// mais il n'y a pas d'erreur, c'est seulement que l'élément n'est pas ad hoc pour
// calculer la normale à ce noeud là
// temps: indique à quel moment on veut le calcul
// pour des éléments particulier (ex: SFE) la méthode est surchargée
int ElemMeca : : CalculNormale_noeud ( Enum_dure temps , const Noeud & noe , Coordonnee & coor )
{
int retour = 1 ; // init du retour : on part d'un bon a priori
Enum_type_geom enutygeom = Type_geom_generique ( this - > Id_geometrie ( ) ) ;
// if ((enutygeom == LIGNE) || (enutygeom == SURFACE))
if ( enutygeom ! = SURFACE ) // dans une première étape on ne s'occupe que des surfaces
{ coor . Libere ( ) ; // pas de calcul possible
retour = 0 ;
}
else // sinon le calcul est possible
{ // on commence par repérer le noeud dans la numérotation locale
int nuoe = 0 ;
int borne_nb_noeud = tab_noeud . Taille ( ) + 1 ;
for ( int i = 1 ; i < borne_nb_noeud ; i + + )
{ Noeud & noeloc = * tab_noeud ( i ) ;
if ( ( noe . Num_noeud ( ) = = noeloc . Num_noeud ( ) )
& & ( noe . Num_Mail ( ) = = noeloc . Num_Mail ( ) )
)
{ nuoe = i ; break ;
} ;
} ;
// on ne continue que si on a trouvé le noeud
if ( nuoe ! = 0 )
{ ElemGeomC0 & elemgeom = ElementGeometrique ( ) ; // récup de la géométrie
// récup des coordonnées locales du noeuds
const Coordonnee & theta_noeud = elemgeom . PtelemRef ( ) ( nuoe ) ;
// récup des phi et dphi au noeud
2023-05-03 17:23:49 +02:00
const Vecteur & phi = elemgeom . Phi_point ( theta_noeud ) ;
const Mat_pleine & dphi = elemgeom . Dphi_point ( theta_noeud ) ;
2021-09-27 12:42:13 +02:00
switch ( temps )
{ case TEMPS_0 :
{ const BaseB & baseB = met - > BaseNat_0 ( tab_noeud , dphi , phi ) ;
coor = Util : : ProdVec_coor ( baseB . Coordo ( 1 ) , baseB . Coordo ( 2 ) ) ;
coor . Normer ( ) ;
break ;
}
case TEMPS_t :
{ const BaseB & baseB = met - > BaseNat_t ( tab_noeud , dphi , phi ) ;
coor = Util : : ProdVec_coor ( baseB . Coordo ( 1 ) , baseB . Coordo ( 2 ) ) ;
coor . Normer ( ) ;
break ;
}
case TEMPS_tdt :
{ const BaseB & baseB = met - > BaseNat_tdt ( tab_noeud , dphi , phi ) ;
coor = Util : : ProdVec_coor ( baseB . Coordo ( 1 ) , baseB . Coordo ( 2 ) ) ;
coor . Normer ( ) ;
break ;
}
default :
cout < < " \n Erreur : valeur incorrecte du temps demande = "
< < Nom_dure ( temps ) < < " ! \n " ;
cout < < " \n ElemMeca::CalculNormale_noeud(Enum_dure temps... \n " ;
retour = 0 ;
Sortie ( 1 ) ;
} ;
}
else
{ cout < < " \n *** erreur le noeud demande num= " < < noe . Num_noeud ( )
< < " du maillage " < < noe . Num_Mail ( )
< < " ne fait pas parti de l'element "
< < num_elt < < " du maillage " < < noe . Num_Mail ( )
< < " on ne peut pas calculer la normale au noeud !! "
< < " \n ElemMeca::CalculNormale_noeud(... " ;
retour = 0 ;
Sortie ( 1 ) ;
}
} ;
// retour
return retour ;
} ;
