2194 lines
114 KiB
C++
Executable file
2194 lines
114 KiB
C++
Executable file
// FICHIER : LoiCritere.cp
|
|
// CLASSE : LoiCritere
|
|
|
|
|
|
// This file is part of the Herezh++ application.
|
|
//
|
|
// The finite element software Herezh++ is dedicated to the field
|
|
// of mechanics for large transformations of solid structures.
|
|
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
|
|
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
|
|
//
|
|
// Herezh++ is distributed under GPL 3 license ou ultérieure.
|
|
//
|
|
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
|
|
// AUTHOR : Gérard Rio
|
|
// E-MAIL : gerardrio56@free.fr
|
|
//
|
|
// This program is free software: you can redistribute it and/or modify
|
|
// it under the terms of the GNU General Public License as published by
|
|
// the Free Software Foundation, either version 3 of the License,
|
|
// or (at your option) any later version.
|
|
//
|
|
// This program is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty
|
|
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
|
// See the GNU General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU General Public License
|
|
// along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
//
|
|
// For more information, please consult: <https://herezh.irdl.fr/>.
|
|
|
|
|
|
//#include "Debug.h"
|
|
#include "LesLoisDeComp.h"
|
|
|
|
# include <iostream>
|
|
using namespace std; //introduces namespace std
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include "Sortie.h"
|
|
#include "TypeConsTens.h"
|
|
#include "TypeQuelconqueParticulier.h"
|
|
#include "NevezTenseurQ.h"
|
|
|
|
#include "LoiCritere.h"
|
|
#include "LoiContraintesPlanesDouble.h"
|
|
#include "LoiContraintesPlanes.h"
|
|
#include "Util.h"
|
|
|
|
|
|
|
|
|
|
// fonction critère de plissement de membrane
|
|
// en entrée:
|
|
// implicit : si oui on est en implicite, sinon en explicite
|
|
// retour:
|
|
// ** ancienne mouture
|
|
// = -2 : erreur en calcul des valeurs propres de déformation, on ne tiens pas compte des plis
|
|
// = -1 : erreur en calcul des valeurs propres de contrainte, on ne tiens pas compte des plis
|
|
// : ou erreur en calcul des vecteurs propres de contrainte
|
|
// = 1 : tout est ok, il n'y a pas d'apparition de plis, rien n'est en compression
|
|
// = 2 : il y a des plis dans les deux sens, aucune raideur ni contrainte
|
|
// = 3 : on a des plis dans une direction
|
|
|
|
// ** nouvelle mouture
|
|
// = -1 : pas de calcul de valeur propre possible en contrainte
|
|
// = 1 : pas de plis (pas de calcul de nouvelle direction )
|
|
// = -2 : pas de calcul de valeur propre de déformation
|
|
// = -3 : plis dans les deux sens, mais pas de calcul de direction propre valide
|
|
// = 2 : plis dans les deux sens, calcul des directions de plis
|
|
// = -4 : pas de calcul de vecteurs propres possible pour les contraintes
|
|
// = 3 : plis dans une seule directions, calcul ok
|
|
|
|
int LoiCritere::Critere_plis_membrane(bool en_base_orthonormee,TenseurHH & sigHH_t,TenseurBB& DepsBB
|
|
,TenseurBB & epsBB_tdt,TenseurBB & delta_epsBB,double& jacobien_0,double& jacobien
|
|
,TenseurHH& sigHH,TenseurHHHH& d_sigma_deps_inter
|
|
,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite
|
|
,double& module_cisaillement
|
|
,bool implicit,const Met_abstraite::Umat_cont& ex)
|
|
{// on sauvegarde les modules de compressibilité et de cisaillement
|
|
// car l'appel à la loi de contrainte plane annulle les valeurs stockées
|
|
// et dans le cas ou il y a un pb, le mieux est de garder les anciennes valeurs
|
|
// double module_compressibilite_save = module_compressibilite;
|
|
// double module_cisaillement_save = module_cisaillement;
|
|
|
|
// dim et vérification
|
|
int dim_tens = abs(sigHH.Dimension());
|
|
// normalement le critère ne s'applique que pour des tenseurs 2D
|
|
#ifdef MISE_AU_POINT
|
|
if (dim_tens != 2)
|
|
{ cout << "\n erreur le critere de plissement de membrane ne peut s'appliquer que pour des lois 2D"
|
|
<< " \n LoiCritere::Critere_plis_membrane_expli(..";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
SaveResul_LoiCritere & save_resul = *((SaveResul_LoiCritere*) saveResul);
|
|
int retour = 0; // init par défaut
|
|
|
|
// pilotage éventuelle du recalcul des directions de plis
|
|
int calcul_dir_plis = 1; // initialisation
|
|
// si rien n'a été auparavant calculé, on est obligé de faire un premier calcul
|
|
// de direction sinon on appelle éventuellement l'indicateur
|
|
|
|
// on récupère l'indicateur: qui est déterminé par une fonction nD
|
|
// si cas_cal_plis == 0 cela veut dire soit qu'il n'y a pas eu déjà un calcul
|
|
// soit que ce précédant calcul ne s'est pas bien passé, dans ces deux il faut
|
|
// imposer de faire le calcul
|
|
if ((recalcul_dir_plis != NULL) && (save_resul.cas_cal_plis != 0))
|
|
{// opération de transmission de la métrique
|
|
const Met_abstraite::Impli* ex_impli = NULL;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Umat_cont* ex_expli = &ex ;
|
|
|
|
// ici on utilise les variables connues aux pti, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal
|
|
= recalcul_dir_plis->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = recalcul_dir_plis->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// 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_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val
|
|
= recalcul_dir_plis->Valeur_FnD_Evoluee
|
|
(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
#ifdef MISE_AU_POINT
|
|
if (tab_val.Taille() != 1)
|
|
{ cout << "\nErreur : la fonction nD relative au recalcul de la direction des plis "
|
|
<< " doit calculer un scalaire or le tableau de retour est de taille "
|
|
<< tab_val.Taille() << " ce n'est pas normal !\n";
|
|
cout << " LoiCritere::Critere_plis_membrane(... \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// on récupère le premier élément du tableau uniquement
|
|
calcul_dir_plis = (int) tab_val(1);
|
|
};
|
|
|
|
// en fonction de l'indicateur calcul conditionnel des plis
|
|
// ou utilisation des grandeurs sauvegardées
|
|
Tableau <Coordonnee2H> V_Pr_H(2);
|
|
Coordonnee2 valPropreSig;
|
|
Coordonnee2 valPropreEps;
|
|
|
|
switch (calcul_dir_plis) // test en fonction d'un pilotage éventuel fct de la fonction nD
|
|
// recalcul_dir_plis->Valeur_FnD_Evoluee
|
|
// rappel de la gestion des différents cas via la valeur de: calcul_dir_plis
|
|
// si == 0 on utilise la valeur stockee actuelle, c-a-d ici à l'itération précédente
|
|
// si aucun calcul de direction n'a ete effectue, il y a un premier calcul qui est fait
|
|
// si == 1 on recalcule les directions des plis,
|
|
// si == 2 on recalcule les directions des plis uniquement pour la zone non complètement
|
|
// relâchée à l'incrément précédent, pour cette dernière on maintient le relâchement
|
|
// si == -1 on utilise la valeur stockee a t (celle de l'incrément precedent qui a convergé),
|
|
// si aucun calcul de direction n'a ete effectue, il y a un premier calcul qui est fait
|
|
// si == -2 idem -1 , sauf que s'il n'y a pas de direction existante, on ne recalcule pas une
|
|
// nouvelle direction, on continue avec un calcul sans plis
|
|
// --> cela signifie que la zone de plis est totalement figée / au précédent incrément
|
|
// si == -3 idem 0 , sauf que s'il n'y a pas de direction existante, on ne recalcule pas une
|
|
// nouvelle direction, on continue avec un calcul sans plis
|
|
// --> cela signifie que la zone de plis est totalement figée / à la précédente itération
|
|
|
|
{ case 2: // cas ou on doit recalculer la direction des plis uniquement pour la zone non relâché
|
|
{ if (save_resul.cas_cal_plis_t != 2)
|
|
{ int retour = Calcul_directions_plis(epsBB_tdt,sigHH,valPropreEps,V_Pr_H,ex,valPropreSig);
|
|
#ifdef MISE_AU_POINT
|
|
// vérif normalement inutile car le retour doit toujours être non nul
|
|
if (!retour)
|
|
{cout << "\n **** erreur a priori: normalement on ne devrait jamais avoir de retour nul "
|
|
<< " or on veut l'utiliser: ce n'est pas possible !! arret !! "
|
|
<< "\n LoiCritere::Critere_plis_membrane(...";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
}
|
|
else // sinon on maintien l'indicateur de relâchement
|
|
{save_resul.cas_cal_plis = save_resul.cas_cal_plis_t;
|
|
};
|
|
|
|
// retour dans le save_resul.cas_cal_plis
|
|
// = -1 : pas de calcul de valeur propre possible en contrainte
|
|
// = 1 : pas de plis (pas de calcul de nouvelle direction )
|
|
// = -2 : pas de calcul de valeur propre de déformation
|
|
// = -3 : plis dans les deux sens, mais pas de calcul de direction propre valide
|
|
// = 2 : plis dans les deux sens, calcul des directions de plis
|
|
// = -4 : pas de calcul de vecteurs propres possible pour les contraintes
|
|
// = 3 : plis dans une seule directions, calcul ok
|
|
// si cela s'est bien passé, les résultats sont dans V_Pr_H
|
|
// sinon on utilisera pas de directions: calculées ou stockées
|
|
break; // fin du cas où on recalcule la direction des plis
|
|
}
|
|
case 1: // cas ou on doit recalculer la direction des plis
|
|
{ int retour = Calcul_directions_plis(epsBB_tdt,sigHH,valPropreEps,V_Pr_H,ex,valPropreSig);
|
|
#ifdef MISE_AU_POINT
|
|
// vérif normalement inutile car le retour doit toujours être non nul
|
|
if (!retour)
|
|
{cout << "\n **** erreur a priori: normalement on ne devrait jamais avoir de retour nul "
|
|
<< " or on veut l'utiliser: ce n'est pas possible !! arret !! "
|
|
<< "\n LoiCritere::Critere_plis_membrane(...";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
// retour dans le save_resul.cas_cal_plis
|
|
// = -1 : pas de calcul de valeur propre possible en contrainte
|
|
// = 1 : pas de plis (pas de calcul de nouvelle direction )
|
|
// = -2 : pas de calcul de valeur propre de déformation
|
|
// = -3 : plis dans les deux sens, mais pas de calcul de direction propre valide
|
|
// = 2 : plis dans les deux sens, calcul des directions de plis
|
|
// = -4 : pas de calcul de vecteurs propres possible pour les contraintes
|
|
// = 3 : plis dans une seule directions, calcul ok
|
|
// si cela s'est bien passé, les résultats sont dans V_Pr_H
|
|
// sinon on utilisera pas de directions: calculées ou stockées
|
|
break; // fin du cas où on recalcule la direction des plis
|
|
}
|
|
case 0: // cas où on se sert de ce qui a été sauvegardé à l'itération précédente
|
|
// si aucun calcul de direction n'a ete effectue, il y a un premier calcul qui est fait
|
|
{ // en fait seul les cas 3 et 2 sont exploitables et seront exploités
|
|
if ((save_resul.cas_cal_plis == 3)||(save_resul.cas_cal_plis == 2))
|
|
{ // V_Pr_H(be)(j) correspond à la coordonnée locale ^j du vecteur be
|
|
// donc V_Pr_H(be)(j) = V_Pr(be) * giH(j)
|
|
// les vecteurs V_Pr_H(be) sont supposés être dans le plan des g_alpha
|
|
for (int be=1;be<3;be++) // < 3 donc de 1 à 2
|
|
{for (int al=1;al<3;al++) // < 3 donc de 1 à 2
|
|
V_Pr_H(be)(al) = (*save_resul.V_P_sig)(be) * ex.giH_tdt->Coordo(al);
|
|
};
|
|
}
|
|
else // sinon cela veut dire que le calcul précédent n'est pas exploitable
|
|
// on refait un calcul
|
|
{int retour = Calcul_directions_plis(epsBB_tdt,sigHH,valPropreEps,V_Pr_H,ex,valPropreSig);
|
|
#ifdef MISE_AU_POINT
|
|
// vérif normalement inutile car le retour doit toujours être non nul
|
|
if (!retour)
|
|
{cout << "\n **** erreur a priori: normalement on ne devrait jamais avoir de retour nul "
|
|
<< " or on veut l'utiliser: ce n'est pas possible !! arret !! "
|
|
<< "\n LoiCritere::Critere_plis_membrane(...";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
};
|
|
break; // fin du cas où on utilise la direction des plis sauvegardée à l'itération prec
|
|
}
|
|
case -1: // cas où on se sert de ce qui a été sauvegardé à l'incrément précédent
|
|
// si aucun calcul de direction n'a ete effectue, il y a un premier calcul qui est fait
|
|
{ // en fait seul les cas 3 et 2 sont exploitables et seront exploités
|
|
if ((save_resul.cas_cal_plis_t == 3)||(save_resul.cas_cal_plis_t == 2))
|
|
{ // V_Pr_H(be)(j) correspond à la coordonnée locale ^j du vecteur be
|
|
// donc V_Pr_H(be)(j) = V_Pr(be) * giH(j)
|
|
// les vecteurs V_Pr_H(be) sont supposés être dans le plan des g_alpha
|
|
// on utilise aussi les vecteurs de la base à t, ce qui permet de convecter
|
|
// les vecteurs des directions de plis
|
|
#ifdef MISE_AU_POINT
|
|
// on vérifie que les données à t existent
|
|
if (save_resul.V_P_sig_t == NULL)
|
|
{cout << "\n **** erreur la direction des plis a t n'est pas encore sauvegardee "
|
|
<< " or on veut l'utiliser: ce n'est pas possible !! arret !! "
|
|
<< "\n LoiCritere::Critere_plis_membrane(...";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
for (int be=1;be<3;be++) // < 3 donc de 1 à 2
|
|
{for (int al=1;al<3;al++) // < 3 donc de 1 à 2
|
|
V_Pr_H(be)(al) = (*save_resul.V_P_sig_t)(be) * ex.giH_t->Coordo(al);
|
|
};
|
|
// on met à jour save_resul.cas_cal_plis en fonction
|
|
save_resul.cas_cal_plis = save_resul.cas_cal_plis_t;
|
|
}
|
|
else // sinon cela veut dire que le calcul précédent n'est pas exploitable
|
|
// on refait un calcul
|
|
{int retour = Calcul_directions_plis(epsBB_tdt,sigHH,valPropreEps,V_Pr_H,ex,valPropreSig);
|
|
#ifdef MISE_AU_POINT
|
|
// vérif normalement inutile car le retour doit toujours être non nul
|
|
if (!retour)
|
|
{cout << "\n **** erreur a priori: normalement on ne devrait jamais avoir de retour nul "
|
|
<< " or on veut l'utiliser: ce n'est pas possible !! arret !! "
|
|
<< "\n LoiCritere::Critere_plis_membrane(...";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
};
|
|
|
|
break; // fin du cas où on utilise la direction des plis sauvegardée à l'increment prec
|
|
}
|
|
case -2: // cas où la zone de plis de l'incrément précédent
|
|
// est totalement figée
|
|
{ // en fait seul les cas 3 et 2 sont exploitables et seront exploités
|
|
if ((save_resul.cas_cal_plis_t == 3)||(save_resul.cas_cal_plis_t == 2))
|
|
{ // V_Pr_H(be)(j) correspond à la coordonnée locale ^j du vecteur be
|
|
// donc V_Pr_H(be)(j) = V_Pr(be) * giH(j)
|
|
// les vecteurs V_Pr_H(be) sont supposés être dans le plan des g_alpha
|
|
// on utilise aussi les vecteurs de la base à t, ce qui permet de convecter
|
|
// les vecteurs des directions de plis
|
|
#ifdef MISE_AU_POINT
|
|
// on vérifie que les données à t existent
|
|
if (save_resul.V_P_sig_t == NULL)
|
|
{cout << "\n **** erreur la direction des plis a t n'est pas encore sauvegardee "
|
|
<< " or on veut l'utiliser: ce n'est pas possible !! arret !! "
|
|
<< "\n LoiCritere::Critere_plis_membrane(...";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
for (int be=1;be<3;be++) // < 3 donc de 1 à 2
|
|
{for (int al=1;al<3;al++) // < 3 donc de 1 à 2
|
|
V_Pr_H(be)(al) = (*save_resul.V_P_sig_t)(be) * ex.giH_t->Coordo(al);
|
|
};
|
|
// on met à jour save_resul.cas_cal_plis en fonction
|
|
save_resul.cas_cal_plis = save_resul.cas_cal_plis_t;
|
|
}
|
|
else // sinon on continue avec un cas sans plis
|
|
{ int retour = Calcul_directions_plis(epsBB_tdt,sigHH,valPropreEps,V_Pr_H,ex,valPropreSig,false);
|
|
// on met à jour save_resul.cas_cal_plis en fonction
|
|
save_resul.cas_cal_plis = retour;
|
|
};
|
|
break; // fin du cas où on utilise la direction des plis sauvegardée à l'increment prec
|
|
}
|
|
case -3: // cas où on se sert de ce qui a été sauvegardé à l'itération précédente
|
|
// idem 0 , sauf que s'il n'y a pas de direction existante, on ne recalcule pas une
|
|
// nouvelle direction, on continue avec un calcul sans plis
|
|
// --> cela signifie que la zone de plis est totalement figée / à la précédente itération
|
|
{ // en fait seul les cas 3 et 2 sont exploitables et seront exploités
|
|
if ((save_resul.cas_cal_plis == 3)||(save_resul.cas_cal_plis == 2))
|
|
{ // V_Pr_H(be)(j) correspond à la coordonnée locale ^j du vecteur be
|
|
// donc V_Pr_H(be)(j) = V_Pr(be) * giH(j)
|
|
// les vecteurs V_Pr_H(be) sont supposés être dans le plan des g_alpha
|
|
for (int be=1;be<3;be++) // < 3 donc de 1 à 2
|
|
{for (int al=1;al<3;al++) // < 3 donc de 1 à 2
|
|
V_Pr_H(be)(al) = (*save_resul.V_P_sig)(be) * ex.giH_tdt->Coordo(al);
|
|
};
|
|
}
|
|
else // sinon on continue avec un cas sans plis
|
|
{ int retour = Calcul_directions_plis(epsBB_tdt,sigHH,valPropreEps,V_Pr_H,ex,valPropreSig,false);
|
|
// on met à jour save_resul.cas_cal_plis en fonction
|
|
save_resul.cas_cal_plis = retour;
|
|
};
|
|
break; // fin du cas où on utilise la direction des plis sauvegardée à l'itération prec
|
|
}
|
|
default:
|
|
cout << "\n erreur non prevu : calcul_dir_plis= "<< save_resul.cas_cal_plis
|
|
<< "\n LoiCritere::Critere_plis_membrane(...";
|
|
Sortie(1);
|
|
};
|
|
|
|
// if ((ParaGlob::NiveauImpression() > 8) || (Permet_affichage() > 5))
|
|
// cout << "\n test3: LoiCritere::Critere_plis_membrane(... " << flush;
|
|
|
|
// maintenant on fait un traitement en fonction de save_resul.cas_cal_plis
|
|
// = -1 : pas de calcul de valeur propre possible en contrainte
|
|
// = 1 : pas de plis (pas de calcul de nouvelle direction )
|
|
// = -2 : pas de calcul de valeur propre de déformation
|
|
// = -3 : plis dans les deux sens, mais pas de calcul de direction propre valide
|
|
// = 2 : plis dans les deux sens, calcul des directions de plis
|
|
// = -4 : pas de calcul de vecteurs propres possible pour les contraintes
|
|
// = 3 : plis dans une seule directions, calcul ok
|
|
retour = save_resul.cas_cal_plis; // init
|
|
switch (save_resul.cas_cal_plis)
|
|
{ case -1: // pas de calcul de valeur propre possible en contrainte
|
|
{save_resul.eps_pli.Zero(); // pas de pli répertorié
|
|
break;
|
|
}
|
|
case 1: // pas de plis
|
|
{break;
|
|
}
|
|
case -2: // pas de calcul de valeur propre de déformation
|
|
{save_resul.eps_pli.Zero(); // pas de pli répertorié
|
|
break;
|
|
}
|
|
// on traite -3 et 2 en même temps, car la fin est commune
|
|
case -3: // plis dans les deux sens, mais pas de calcul de direction propre valide
|
|
case 2: // plis dans les deux sens, calcul des directions de plis
|
|
{if (save_resul.cas_cal_plis == -3)
|
|
{save_resul.eps_pli.Zero(); // pas de pli répertorié
|
|
}
|
|
else // -> cas: save_resul.cas_cal_plis == 2
|
|
{if (calcul_dir_plis != 2)
|
|
{ // calcul normal relatif à la zone relâchée
|
|
// les coordonnées des vecteurs propres sont exprimées dans l'ancienne base
|
|
// on a donc
|
|
// Vi_B(i) = gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
|
|
Mat_pleine beta(2,2);
|
|
for (int i=1;i<3;i++) // < 3 donc de 1 à 2
|
|
for (int j=1;j<3;j++)
|
|
beta(i,j) = V_Pr_H(i)(j);
|
|
Mat_pleine gamma(beta.Inverse().Transpose()); // on calcule la matrice inverse
|
|
// la matrice gamma correspond à la matrice de changement de base des g^i au gpH
|
|
// gamma est l'inverse transposée
|
|
// gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
|
|
// c-a-d= gp^i = gamma^i_j * g^j
|
|
bool deux_plis = true;
|
|
// on calcul les def de plis : save_resul.eps_pli
|
|
Passage_3ou2_vers_1(gamma,sigHH_t,beta,DepsBB,epsBB_tdt
|
|
,delta_epsBB,deux_plis,save_resul.eps_pli);
|
|
}
|
|
else // sinon calcul_dir_plis == 2, cela veut dire que l'on n'a pas recalculé
|
|
// dans la zône relâchée donc on utilise les déformations déjà calculées à t
|
|
{ save_resul.eps_pli = save_resul.eps_pli_t;
|
|
// on retient les directions des plis précendentes
|
|
(*save_resul.V_P_sig) = (*save_resul.V_P_sig_t);
|
|
};
|
|
};
|
|
// traitement des épaisseurs et des largeurs
|
|
// LoiContraintesPlanes::SaveResul_LoiContraintesPlanes * save_result_plan_inter // pour simplifier
|
|
// = (LoiContraintesPlanes::SaveResul_LoiContraintesPlanes *)
|
|
// (*(save_resul.liste_des_SaveResul.begin()));
|
|
LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble * save_result_plan_inter // pour simplifier
|
|
= (LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble *)
|
|
(*(save_resul.liste_des_SaveResul.begin()));
|
|
switch (choix_calcul_epaisseur_si_relachement_complet)
|
|
{ case 1: // cas où on impose un retour à l'épaisseur et à la largeur initiales
|
|
save_result_plan_inter->hsurh0 = 1.;
|
|
save_result_plan_inter->bsurb0 = 1.;
|
|
break;
|
|
case 2: // cas où on impose l'épaisseur et la largeur du pas précédent
|
|
{if (save_result_plan_inter->h_tsurh0 != 0.)
|
|
{save_result_plan_inter->hsurh0 = save_result_plan_inter->h_tsurh0;}
|
|
else // sinon on remet à 1.
|
|
{save_result_plan_inter->hsurh0 = 1.;};
|
|
if (save_result_plan_inter->b_tsurb0 != 0.)
|
|
{save_result_plan_inter->bsurb0 = save_result_plan_inter->b_tsurb0;}
|
|
else // sinon on remet à 1.
|
|
{save_result_plan_inter->bsurb0 = 1.;};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur *** la variable choix_calcul_epaisseur_si_relachement_complet= "
|
|
<< choix_calcul_epaisseur_si_relachement_complet
|
|
<< " est differente de 1= le seul choix actuellement implante " << endl;
|
|
Sortie(1);
|
|
break;
|
|
};
|
|
// on annule la contrainte et sa contribution éventuelle à la raideur
|
|
sigHH.Inita(0.);
|
|
//d_sigma_deps_inter.Inita(0.);
|
|
// on annule également les énergies
|
|
energ.Inita(0.);
|
|
break;
|
|
}
|
|
case -4: // pas de calcul de vecteurs propres possible pour les contraintes
|
|
{save_resul.eps_pli.Zero(); // pas de pli répertorié
|
|
break;
|
|
}
|
|
case 3: // plis dans une seule directions
|
|
{// va être traité en dehors du switch car c'est le gros morceau
|
|
break;
|
|
}
|
|
|
|
default: // normalement n'arrive jamais car Calcul_directions_plis a toujours une solution diff de 0
|
|
{ cout << "\n *** cas d'erreur inconnue dans le calcul de la direction des plis "
|
|
<< " save_resul.cas_cal_plis= "<< save_resul.cas_cal_plis ;
|
|
Sortie(1);
|
|
break;
|
|
}
|
|
};
|
|
|
|
// if ((ParaGlob::NiveauImpression() > 8) || (Permet_affichage() > 5))
|
|
// cout << "\n test4: LoiCritere::Critere_plis_membrane(... " << flush;
|
|
// le gros morceau: cas où on a une seule direction
|
|
if (save_resul.cas_cal_plis == 3)
|
|
{// choix entre la méthode historique et la nouvelle méthode
|
|
switch (choix_methode_cal_plis_memb)
|
|
{case 1:
|
|
{Premie_type_calcul_en_un_pli(epsBB_tdt,delta_epsBB
|
|
,sigHH_t,jacobien_0,jacobien,energ,energ_t
|
|
,DepsBB,module_compressibilite,module_cisaillement
|
|
,V_Pr_H,sigHH,d_sigma_deps_inter
|
|
,valPropreSig,ex);
|
|
break;
|
|
}
|
|
case 2:
|
|
{Deuxieme_type_calcul_en_un_pli(epsBB_tdt,delta_epsBB
|
|
,sigHH_t,jacobien_0,jacobien,energ,energ_t
|
|
,DepsBB,module_compressibilite,module_cisaillement
|
|
,V_Pr_H,sigHH,d_sigma_deps_inter
|
|
,valPropreSig,ex);
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n cas non prevu: choix_methode_cal_plis_memb= "
|
|
<< choix_methode_cal_plis_memb ;
|
|
Sortie(1);
|
|
break;
|
|
};
|
|
|
|
|
|
};
|
|
|
|
return retour;
|
|
};
|
|
|
|
// calcul d'une nouvelle direction de plis
|
|
// en entrée: force: = true par défaut
|
|
// si == false -> pas de calcul de direction demandée, et mise en place
|
|
// des indicateurs en cohérence avec le fait de ne pas faire de calcul
|
|
// en retour:
|
|
// = -1 : pas de calcul de valeur propre possible en contrainte
|
|
// = 1 : pas de plis (pas de calcul de nouvelle direction )
|
|
// = -2 : pas de calcul de valeur propre de déformation
|
|
// = -3 : plis dans les deux sens, mais pas de calcul de direction propre valide
|
|
// = 2 : plis dans les deux sens, calcul des directions de plis
|
|
// = -4 : pas de calcul de vecteurs propres possible pour les contraintes
|
|
// = 3 : plis dans une seule directions, calcul ok
|
|
// = 0 : erreur inconnue ??
|
|
|
|
int LoiCritere::Calcul_directions_plis(const TenseurBB & epsBB_tdt,const TenseurHH& sigHH
|
|
,Coordonnee2& valPropreEps
|
|
,Tableau <Coordonnee2H>& V_Pr_H
|
|
,const Met_abstraite::Umat_cont& ex
|
|
,Coordonnee2& valPropreSig
|
|
, bool force
|
|
)
|
|
{// dim et vérification
|
|
int dim_tens = abs(sigHH.Dimension());
|
|
V_Pr_H.Change_taille(2);
|
|
// normalement le critère ne s'applique que pour des tenseurs 2D
|
|
#ifdef MISE_AU_POINT
|
|
if (dim_tens != 2)
|
|
{ cout << "\n erreur le critere de plissement de membrane ne peut s'appliquer que pour des lois 2D"
|
|
<< " \n LoiCritere::Calcul_directions_plis(..";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
SaveResul_LoiCritere & save_resul = *((SaveResul_LoiCritere*) saveResul);
|
|
int retour = 0; // init par défaut
|
|
|
|
if (!force)
|
|
// cas particulier ou on ne veut pas de calcul, mais par contre on veut mettre les indicateurs
|
|
// cohérence avec le fait de ne pas faire de calcul
|
|
{ retour = 1; //on signale
|
|
save_resul.eps_pli.Zero(); // pas de pli répertorié
|
|
// comme il n'y a pas de plis, donc par de variation de largeur
|
|
LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble * save_result_1DCP2_inter // pour simplifier
|
|
= (LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble *)
|
|
(*(save_resul.liste_des_SaveResul.begin()));
|
|
save_result_1DCP2_inter->Zero_var_largeur(); // mise à 0 de b_surb0
|
|
}
|
|
else
|
|
{// on va calculer les valeurs propres relatives au tenseur des contraintes
|
|
|
|
TenseurHB* sigHB = (NevezTenseurHB(dim_tens,0.));
|
|
(*sigHB) = sigHH * (* ex.gijBB_tdt);
|
|
int casSig = 0; // init par défaut
|
|
// calcul des valeurs propres correspondantes à des vecteurs propres CoordonneeH
|
|
valPropreSig = sigHB->ValPropre(casSig);
|
|
// si le calcul n'est pas possible on ne fait rien mais on poste un message
|
|
if (casSig == -1)
|
|
// ==== cas d'une erreur de valeur propre sur la contrainte ====
|
|
{ if (Permet_affichage() > 4)
|
|
cout << "\n warning : le calcul des valeurs propres de contrainte n'a pas ete possible,"
|
|
<< " l'application du critere de plissement n'est pas possible, et donc non effectuee " << endl;
|
|
retour = -1;
|
|
save_resul.eps_pli.Zero(); // pas de pli répertorié
|
|
}
|
|
else if (valPropreSig(dim_tens) >= save_resul.niveau_declenchement_critere) // 0.)
|
|
{ // on se trouve dans le cas où rien n'est en compression, donc cas normal sans plissement
|
|
retour = 1; //on signale
|
|
save_resul.eps_pli.Zero(); // pas de pli répertorié
|
|
// comme il n'y a pas de plis, donc par de variation de largeur
|
|
LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble * save_result_1DCP2_inter // pour simplifier
|
|
= (LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble *)
|
|
(*(save_resul.liste_des_SaveResul.begin()));
|
|
save_result_1DCP2_inter->Zero_var_largeur(); // mise à 0 de b_surb0
|
|
}
|
|
else
|
|
{ // ==== on se trouve dans le cas où potentiellement on peut avoir des plis ====
|
|
// on regarde les valeurs propres en déformation
|
|
|
|
TenseurHB* epsHB = (NevezTenseurHB(dim_tens,0.));
|
|
(*epsHB) = (* ex.gijHH_tdt) * epsBB_tdt;
|
|
int casEps = 0; // init par défaut
|
|
// calcul des valeurs propres correspondantes à des vecteurs propres CoordonneeH
|
|
valPropreEps = epsHB->ValPropre(casEps);
|
|
if (casEps == -1)
|
|
{ if (Permet_affichage() > 4)
|
|
cout << "\n warning : le calcul des valeurs propres de deformation n'a pas ete possible,"
|
|
<< " l'application du critere de plissement n'est pas possible, et donc non effectuee " << endl;
|
|
retour = -2;
|
|
save_resul.eps_pli.Zero(); // pas de pli répertorié
|
|
}
|
|
else if (valPropreEps(1) < 0.)
|
|
// ==== cas de plis dans les deux sens ====
|
|
{// si la plus grande des valeurs propres de déformation est négative,
|
|
// dans ce cas il y a des plis dans les deux sens il faut que l'on récupère
|
|
// les def des plis
|
|
//*** il faudra grouper ce if avec le cas d'un pli de manière à éviter les erreurs de duplication ***
|
|
// -- on doit définir un repère ortho dont le premier vecteur est celui selon lequel on fait du 1D
|
|
// calcul des vecteurs propres
|
|
int cas = 0; // init
|
|
|
|
Tableau <Coordonnee> V_inter;
|
|
sigHB->VecteursPropres(valPropreSig,cas, V_inter);
|
|
// si cas = -1, la construction des vecteurs propres est impossible, on arrête le traitement
|
|
if (cas == -1)
|
|
{save_resul.eps_pli.Zero(); // on ne peut pas donner une intensité au def de pli !
|
|
// les vecteurs propres sauvegardées sont nulles
|
|
Coordonnee3 zero3; // un vecteur nul
|
|
save_resul.V_P_sig->Change_taille(3,zero3); // en 3D qui est obligatoire ici
|
|
retour = -3;
|
|
}
|
|
else
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
if (Dabs(V_inter(1) * V_inter(2)) > ConstMath::petit)
|
|
{cout << "\n== LoiCritere::Critere_plis_membrane: plis dans les deux sens ";
|
|
//cout << "\n sigma: vecteurs propres: \n "<< V_inter;
|
|
cout << "\n sigma: valeurs propres: "; valPropreSig.Affiche_1(cout);
|
|
// on va vérifier les vecteurs propres
|
|
Coordonnee2H VpH1; VpH1.Change_val(V_inter(1));
|
|
cout << "\n vecteur propre 1 VpH1= " << VpH1;
|
|
Coordonnee2H VpropreH1 = (*sigHB) * VpH1;
|
|
//cout << "\n sigHH * VH1 "<< VpropreH1;
|
|
// comme les valeurs peuvent être très petites on régularise
|
|
VpropreH1 /= (valPropreSig(1)+ConstMath::petit);
|
|
cout << "\n verif: (sigHB * VH1)/(lambda1+epsilon) = "<< VpropreH1;
|
|
CoordonneeB aB; aB.Change_val(V_inter(2));
|
|
Coordonnee2H VpH2 = aB * (* ex.gijHH_tdt);
|
|
cout << "\n vecteur propre 2 VpH2 "<< VpH2;
|
|
Coordonnee2H VpropreH2 = (*sigHB) * VpH2;
|
|
//cout << "\n sigHH * VH2 "<< VpropreH2;
|
|
VpropreH2 /= (valPropreSig(2)+ConstMath::petit);
|
|
cout << "\n verif: (sigHB * VH2)/(lambda2+epsilon) = "<< VpropreH2;
|
|
cout << endl;
|
|
};
|
|
#endif
|
|
|
|
// Le premier vecteur propre s'exprime dans la base naturelle, par contre le second s'exprime
|
|
// dans la base duale donc : V1^H, et V2_B
|
|
// donc si on veut un changement de base de naturelle en naturelle il faut calculer V2^H
|
|
V_Pr_H(1).Change_val(V_inter(1)); // le premier vecteur ok
|
|
CoordonneeB aB; aB.Change_val(V_inter(2));
|
|
V_Pr_H(2)= aB * (* ex.gijHH_tdt);
|
|
// pour l'instant ces vecteurs propres ne sont pas normés car les vecteurs
|
|
// de base ne sont pas normés, on va donc les normer
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{cout << "\n normalisation des vecteurs propres:";
|
|
};
|
|
#endif
|
|
// on calcul les vecteurs propres en orthonormee
|
|
// on définit l'espace de stockage des vecteurs propres si nécessaire
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{cout << "\n vecteurs propres en orthonormee:";
|
|
};
|
|
#endif
|
|
if (save_resul.V_P_sig == NULL)
|
|
save_resul.V_P_sig = new Tableau <Coordonnee>;
|
|
Coordonnee3 zero3; // un vecteur nul
|
|
save_resul.V_P_sig->Change_taille(3,zero3); // on va calculer en 3D qui est obligatoire ici
|
|
|
|
for (int be=1;be<3;be++) // be < 3 donc de 1 à 2
|
|
for (int al=1;al<3;al++)
|
|
(*save_resul.V_P_sig)(be) += V_Pr_H(be)(al) * ex.giB_tdt->Coordo(al);
|
|
// on norme les vecteurs après le passage en absolue (modif prise en compte)
|
|
for (int al=1;al<3;al++) // al < 3 erreur
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{cout << "\n vecteur propre ("<<al<<") local avant normalisation " << V_Pr_H(al);
|
|
cout << "\n vecteur propre ("<<al<<") global avant normalisation " << (*save_resul.V_P_sig)(al);
|
|
double norme2 = sqrt(V_Pr_H(al) * ((* ex.gijBB_tdt) * V_Pr_H(al)));
|
|
cout << "\n norme locale pour verif (doit etre = norme globale): "<< norme2;
|
|
};
|
|
#endif
|
|
double norme2 = (*save_resul.V_P_sig)(al).Norme();
|
|
(*save_resul.V_P_sig)(al) /= (norme2+ConstMath::petit);
|
|
// on norme également les vecteurs V_Pr_H qui vont être utilisés par le programme appelant
|
|
V_Pr_H(al) /= (norme2+ConstMath::petit);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{cout << "\n vecteur propre local apres normalisation " << V_Pr_H(al);
|
|
cout << "\n norme VPH("<<al<<") en globale avant normalisation "<< norme2;
|
|
cout << "\n vecteur propre ("<<al<<") global apres normalisation " << (*save_resul.V_P_sig)(al);
|
|
};
|
|
#endif
|
|
};
|
|
(*save_resul.V_P_sig)(3) = Util::ProdVec_coor((*save_resul.V_P_sig)(1),(*save_resul.V_P_sig)(2));
|
|
(*save_resul.V_P_sig)(3).Normer();// on norme le troisième vecteur *** normalement ne devrait pas être nécessaire
|
|
retour = 2;
|
|
};
|
|
}
|
|
else
|
|
{// cas où on passe d'un comportement 2D à 1D
|
|
Loi_comp_abstraite* loi = *(lois_internes.begin()); // pour simplifier
|
|
// -- donc maintenant la suite concerne le cas où seule une direction est en traction
|
|
// on va se mettre en 1D dans le repère principal de contrainte
|
|
|
|
// -- on doit définir un repère ortho dont le premier vecteur est celui selon lequel on fait du 1D
|
|
// calcul des vecteurs propres
|
|
int cas = 0; // init
|
|
Tableau <Coordonnee> V_inter;
|
|
sigHB->VecteursPropres(valPropreSig,cas, V_inter);
|
|
// si cas = -1, la construction des vecteurs propres est impossible, on arrête le traitement
|
|
if (cas == -1)
|
|
{ if (Permet_affichage() > 4)
|
|
cout << "\n warning : le calcul des valeurs propres de contrainte n'a pas ete possible,"
|
|
<< " l'application du critere de plissement n'est pas possible, et donc non effectuee " << endl;
|
|
retour = -4;
|
|
save_resul.eps_pli.Zero(); // pas de pli répertorié
|
|
// les vecteurs propres sauvegardées sont nulles
|
|
Coordonnee3 zero3; // un vecteur nul
|
|
save_resul.V_P_sig->Change_taille(3,zero3); // en 3D qui est obligatoire ici
|
|
}
|
|
else
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{cout << "\n== LoiCritere::Critere_plis_membrane plis dans un seul sens ";
|
|
//cout << "\n sigma: vecteurs propres: \n "<< V_inter;
|
|
cout << "\n sigma: valeurs propres: "; valPropreSig.Affiche_1(cout);
|
|
// on va vérifier les vecteurs propres
|
|
Coordonnee2H VpH1; VpH1.Change_val(V_inter(1));
|
|
cout << "\n vecteur propre 1 VpH1= " << VpH1;
|
|
Coordonnee2H VpropreH1 = (*sigHB) * VpH1;
|
|
//cout << "\n sigHH * VH1 "<< VpropreH1;
|
|
// comme les valeurs peuvent être très petites on régularise
|
|
VpropreH1 /= (valPropreSig(1)+ConstMath::petit);
|
|
cout << "\n verif: (sigHB * VH1)/(lambda1+epsilon) = "<< VpropreH1;
|
|
CoordonneeB aB; aB.Change_val(V_inter(2));
|
|
Coordonnee2H VpH2 = aB * (* ex.gijHH_tdt);
|
|
cout << "\n vecteur propre 2 VpH2 "<< VpH2;
|
|
Coordonnee2H VpropreH2 = (*sigHB) * VpH2;
|
|
//cout << "\n sigHH * VH2 "<< VpropreH2;
|
|
VpropreH2 /= (valPropreSig(2)+ConstMath::petit);
|
|
cout << "\n verif: (sigHB * VH2)/(lambda2+epsilon) = "<< VpropreH2;
|
|
cout << endl;
|
|
};
|
|
#endif
|
|
|
|
// Le premier vecteur propre s'exprime dans la base naturelle, par contre le second s'exprime
|
|
// dans la base duale donc : V1^H, et V2_B
|
|
// donc si on veut un changement de base de naturelle en naturelle il faut calculer V2^H
|
|
V_Pr_H(1).Change_val(V_inter(1)); // le premier vecteur ok
|
|
CoordonneeB aB; aB.Change_val(V_inter(2));
|
|
V_Pr_H(2)= aB * (* ex.gijHH_tdt);
|
|
// pour l'instant ces vecteurs propres ne sont pas normés car les vecteurs
|
|
// de base ne sont pas normés, on va donc les normer
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{cout << "\n normalisation des vecteurs propres:";
|
|
};
|
|
#endif
|
|
// on calcul les vecteurs propres en orthonormee
|
|
// on définit l'espace de stockage des vecteurs propres si nécessaire
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() >5)
|
|
{cout << "\n vecteurs propres en orthonormee:";
|
|
};
|
|
#endif
|
|
if (save_resul.V_P_sig == NULL) save_resul.V_P_sig = new Tableau <Coordonnee>;
|
|
Coordonnee3 zero3; // un vecteur nul
|
|
save_resul.V_P_sig->Change_taille(3,zero3); // on va calculer en 3D
|
|
|
|
for (int be=1;be<3;be++) // be < 3 donc de 1 à 2
|
|
for (int al=1;al<3;al++)
|
|
(*save_resul.V_P_sig)(be) += V_Pr_H(be)(al) * ex.giB_tdt->Coordo(al);
|
|
// on norme les vecteurs après le passage en absolue
|
|
for (int al=1;al<3;al++) // al < 3 erreur
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{cout << "\n vecteur propre local avant normalisation " << V_Pr_H(al);
|
|
cout << "\n vecteur propre ("<<al<<") global avant normalisation " << (*save_resul.V_P_sig)(al);
|
|
double norme2 = sqrt(V_Pr_H(al) * ((* ex.gijBB_tdt) * V_Pr_H(al)));
|
|
cout << "\n norme locale pour verif (doit etre = norme globale): "<< norme2;
|
|
};
|
|
#endif
|
|
double norme2 = (*save_resul.V_P_sig)(al).Norme();
|
|
(*save_resul.V_P_sig)(al) /= ((norme2)+ConstMath::petit);
|
|
// on norme également les vecteurs V_Pr_H qui vont être utilisés par le programme appelant
|
|
V_Pr_H(al) /= (norme2+ConstMath::petit);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{cout << "\n vecteur propre local apres normalisation " << V_Pr_H(al);
|
|
cout << "\n norme VPH("<<al<<") en globale avant normalisation "<< norme2;
|
|
cout << "\n vecteur propre ("<<al<<") global apres normalisation " << (*save_resul.V_P_sig)(al);
|
|
};
|
|
#endif
|
|
};
|
|
(*save_resul.V_P_sig)(3) = Util::ProdVec_coor((*save_resul.V_P_sig)(1),(*save_resul.V_P_sig)(2));
|
|
(*save_resul.V_P_sig)(3).Normer();// on norme le troisième vecteur *** normalement ne devrait pas être nécessaire
|
|
retour = 3;
|
|
};
|
|
};
|
|
// on supprime les tenseurs provisoires
|
|
delete epsHB;
|
|
};
|
|
// on supprime les tenseurs provisoires
|
|
delete sigHB;
|
|
};
|
|
|
|
// mise à jour de l'indicateur dans save_resul
|
|
save_resul.cas_cal_plis=retour;
|
|
// affichage conditionnelle
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 4)
|
|
{bool erreur = false;
|
|
switch (retour)
|
|
{ case -1 : cout << "\n pas de calcul de valeur propre possible en contrainte ";erreur=true;break;
|
|
case -2 : cout << "\n pas de calcul de valeur propre possible de deformation "; erreur=true;break;
|
|
case -3 : cout << "\n plis dans les deux sens, mais pas de calcul de direction propre valide "; erreur=true;break;
|
|
case -4 : cout << "\n pas de calcul de vecteurs propres possible pour les contraintes "; erreur=true;break;
|
|
case 0 : cout << "\n erreur inconnue dans le calcul directions plis ?? "; erreur=true;break;
|
|
default: break; // cas normal, aucun affichage
|
|
};
|
|
if (erreur)
|
|
if (ptintmeca_en_cours != NULL)
|
|
cout << " mail: "<< ptintmeca_en_cours->Nb_mail() << " ele: "<< ptintmeca_en_cours->Nb_ele()
|
|
<< " npti: "<< ptintmeca_en_cours->Nb_pti() ;
|
|
};
|
|
#endif
|
|
|
|
// retour
|
|
return retour;
|
|
};
|
|
|
|
|
|
// création du conteneur UMAT a partir des vecteurs propres
|
|
void LoiCritere::Creation_metrique_a_partir_vecteurs_propres_pour_Umat1D
|
|
(const Met_abstraite::Umat_cont& ex,const Mat_pleine& beta )
|
|
{
|
|
// on va se mettre en 1D dans le repère principal de contrainte
|
|
SaveResul_LoiCritere & save_resul = *((SaveResul_LoiCritere*) saveResul);
|
|
|
|
// la partie dépendant des vitesses: entre accolades pour pouvoir fermer
|
|
// ici il est nécessaire de recalculer les gradients dans la nouvelle base
|
|
// dans le cas où ils existent
|
|
{if (ex.gradVmoyBB_t != NULL)
|
|
{umat_cont_1D->gradVmoyBB_t= gradVmoyBB_t_1D_P = &gradVmoyBB_t_1D;
|
|
// on change de base
|
|
TenseurBB* t_interBB = NevezTenseurBB(*ex.gradVmoyBB_t);
|
|
t_interBB->ChBase(beta);
|
|
// puis on transfert en 1D
|
|
umat_cont_1D->gradVmoyBB_t->Affectation_trans_dimension(*t_interBB,true);
|
|
delete t_interBB;
|
|
};
|
|
if (ex.gradVmoyBB_tdt != NULL)
|
|
{umat_cont_1D->gradVmoyBB_tdt=gradVmoyBB_tdt_1D_P = &gradVmoyBB_tdt_1D;
|
|
// on change de base
|
|
TenseurBB* t_interBB = NevezTenseurBB(*ex.gradVmoyBB_tdt);
|
|
t_interBB->ChBase(beta);
|
|
// puis on transfert en 1D
|
|
umat_cont_1D->gradVmoyBB_tdt->Affectation_trans_dimension(*t_interBB,true);
|
|
delete t_interBB;
|
|
};
|
|
if (ex.gradVBB_tdt != NULL)
|
|
{umat_cont_1D->gradVBB_tdt=gradVBB_tdt_1D_P = &gradVBB_tdt_1D;
|
|
// on change de base
|
|
TenseurBB* t_interBB = NevezTenseurBB(*ex.gradVBB_tdt);
|
|
t_interBB->ChBase(beta);
|
|
// puis on transfert en 1D
|
|
umat_cont_1D->gradVBB_tdt->Affectation_trans_dimension(*t_interBB,true);
|
|
delete t_interBB;
|
|
};
|
|
}; // fin de la partie dédiée à la vitesse
|
|
|
|
// -- maintenant on construit les éléments de l'Umat pour le repère ortho avec 1 seul
|
|
// vecteur de base, et en orthonormé
|
|
// en fait on va construire un repère orthonormé à partir des vecteurs propres
|
|
// **** comme on considère que la base locale est orthonormé, normalement elle ne change pas
|
|
// pour l'instant on laisse sa construction, mais en fait une seule est suffisante
|
|
// ***** 2 oct 2018: je commente la ligne qui suit car en fait j'ai l'impression
|
|
// que la base change puisqu'elle vaut les directions propres et que ceux-ci dépendent
|
|
// du calcul .... là je ne comprends pas pourquoi j'ai mis le test ....
|
|
// if ((*(umat_cont_1D->gijBB_0))(1,1) == 0.) // pour savoir si la construction a été faite on non
|
|
{ // les vecteurs propres sont considérés normés: cf. voir leurs calculs
|
|
BaseB& giB_tdt = *(umat_cont_1D->giB_tdt); // pour simplifier
|
|
const Coordonnee& V_sig = (*save_resul.V_P_sig)(1); // pour simplifier
|
|
giB_tdt.CoordoB(1).Change_val(V_sig); // on utilise un vecteur 3D, ici dans la base orthonormee
|
|
// giB_tdt(1) = V_Pr_B(1); //.Change_val(V_sig);
|
|
// giB_tdt(1)(1) = 1.;
|
|
// idem pour les temps 0 et t
|
|
*(umat_cont_1D->giB_0)=giB_tdt;
|
|
*(umat_cont_1D->giB_t)=giB_tdt;
|
|
// idem pour le contravariant
|
|
BaseH& giH_tdt = *(umat_cont_1D->giH_tdt); // pour simplifier
|
|
// comme la base est orthonormée, H et B sont identiques
|
|
giH_tdt.CoordoH(1).Change_val(V_sig);
|
|
// giH_tdt(1).Change_val(V_Pr_B(1).Coor());
|
|
// giH_tdt(1)(1) = 1.;
|
|
*(umat_cont_1D->giH_0)=giH_tdt;
|
|
*(umat_cont_1D->giH_t)=giH_tdt;
|
|
// maintenant les métriques : ici identité
|
|
(*(umat_cont_1D->gijBB_0)).Coor(1,1)=(*(umat_cont_1D->gijHH_0)).Coor(1,1)=1.;
|
|
(*(umat_cont_1D->gijBB_t)).Coor(1,1)=(*(umat_cont_1D->gijHH_t)).Coor(1,1)=1.;
|
|
(*(umat_cont_1D->gijBB_tdt)).Coor(1,1)=(*(umat_cont_1D->gijHH_tdt)).Coor(1,1)=1.;
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{cout << "\n LoiCritere::Creation_metrique_a_partir_vecteurs_propres_pour_Umat1D( ";
|
|
cout << "\n base gipB(1) "<< giB_tdt(1)
|
|
<< "\n base gipH(1) "<< giH_tdt(1);
|
|
};
|
|
#endif
|
|
};
|
|
};
|
|
|
|
// passage des grandeurs métriques de l'ordre 3 ou 2 à 1: cas implicite
|
|
void LoiCritere::Passage_metrique_ordre_3_2_vers_1(const Met_abstraite::Umat_cont& ex,const Mat_pleine& gamma)
|
|
{};
|
|
|
|
|
|
// chgt de repère et de dimension pour les variables de passage pour le calcul final de la loi en 1D
|
|
void LoiCritere::Passage_3ou2_vers_1(const Mat_pleine& gamma,const TenseurHH & sigHH_t
|
|
,const Mat_pleine& beta
|
|
,const TenseurBB& DepsBB
|
|
,const TenseurBB & epsBB_tdt,const TenseurBB & delta_epsBB
|
|
,const bool& deux_plis,Coordonnee2& eps_pli)
|
|
{ // on s'occupe du changement de base
|
|
int dim_tens = abs(sigHH_t.Dimension());
|
|
// on crée des nouveaux tenseurs au lieu d'utiliser les tenseurs initiaux, car on fait un changement de base
|
|
// qui a pour effet de modifier le tenseur. Donc si ensuite on veut utiliser les tenseurs initiaux, ce ne sera
|
|
// plus possible comme initialement
|
|
if (dim_tens == 2)
|
|
{ Tenseur2HH sigHH_2_inter(sigHH_t); sigHH_2_inter.ChBase(gamma);
|
|
sig_HH_t_1D.Affectation_trans_dimension(sigHH_2_inter,false);// passage des tenseurs résultats à l'ordre 1
|
|
Tenseur2BB DepsBB_2_inter(DepsBB); DepsBB_2_inter.ChBase(beta);
|
|
Deps_BB_1D.Affectation_trans_dimension(DepsBB_2_inter,false);
|
|
Tenseur2BB epsBB_tdt_2_inter(epsBB_tdt);epsBB_tdt_2_inter.ChBase(beta);
|
|
/*
|
|
|
|
Lorsqu'il y a un pli dans une seule direction, le premier vecteur indique la direction selon laquelle la membrane est en traction. Son intensité est égale à la déformation mécanique dans cette même direction. Le second vecteur est normal à la première direction, son intensité est égale à la déformation cinématique selon cette direction, c'est-à-dire calculée à l'aide des déplacements des noeuds du maillage. Cette déformation est différente de la déformation mécanique calculée par la loi 3D, c'est-à-dire issue de la contraction de la zone, due à la traction suivant la première direction. Cette déformation mécanique peut-être récupéré via la deformation "DEF\_LARGEUR" associée à la loi de comportement.
|
|
|
|
Lorsqu'il y a des plis dans deux directions, les deux vecteurs indiquent les directions principales d'apparition de plis. L'intensité des vecteurs est égale à la déformation cinématique dans la direction associée.
|
|
|
|
Lorsqu'il n'y a pas de plis, l'intensité des vecteurs est nulle.
|
|
|
|
*/
|
|
|
|
eps_pli(1) = epsBB_tdt_2_inter(2,2); // on sauvegarde les def de plis
|
|
if (deux_plis) // cas où il y aurait des plis dans les deux sens
|
|
eps_pli(2) = epsBB_tdt_2_inter(1,1);
|
|
else eps_pli(2) = 0.; // cas d'un seul plis
|
|
eps_BB_1D.Affectation_trans_dimension(epsBB_tdt_2_inter,false);
|
|
Tenseur2BB delta_epsBB_2_inter(delta_epsBB); delta_epsBB_2_inter.ChBase(beta);
|
|
delta_eps_BB_1D.Affectation_trans_dimension(delta_epsBB_2_inter,false);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{cout << "\n debug LoiCritere::Passage_3ou2_vers_1( ";
|
|
cout << "\n sigHH_t= en g_i "; sigHH_t.Ecriture(cout);
|
|
cout << "\n sigHH_t= en vecteurs propres "<< sigHH_2_inter;
|
|
cout << "\n epsBB_tdt en gi ";epsBB_tdt.Ecriture(cout);
|
|
cout << "\n epsBB_tdt en vecteurs propres "<< epsBB_tdt_2_inter;
|
|
cout << "\n epsBB 1D en local vecteurs propres " << eps_BB_1D;
|
|
// en orthonormee
|
|
Tenseur3BB epsBB_3_inter;
|
|
eps_BB_1D.BaseAbsolue(epsBB_3_inter,(giH_tdt_1D));
|
|
cout << "\n base giH_tdt_1D: "<<giH_tdt_1D;
|
|
cout << "\n epsBB 1D =en orthonormee (contient que l'action de la composante de traction) \n "
|
|
<<epsBB_3_inter;
|
|
cout << "\n delta_epsBB en vecteurs propres= "<< delta_epsBB_2_inter;
|
|
cout << endl;
|
|
};
|
|
#endif
|
|
}
|
|
else // cas 3D
|
|
{ Tenseur3HH sigHH_3_inter(sigHH_t); sigHH_3_inter.ChBase(gamma);
|
|
sig_HH_t_1D.Affectation_trans_dimension(sigHH_3_inter,false);// passage des tenseurs résultats à l'ordre 1
|
|
Tenseur3BB DepsBB_3_inter(DepsBB); DepsBB_3_inter.ChBase(beta);
|
|
Deps_BB_1D.Affectation_trans_dimension(DepsBB_3_inter,false);
|
|
Tenseur3BB epsBB_tdt_3_inter(epsBB_tdt);epsBB_tdt_3_inter.ChBase(beta);
|
|
eps_BB_1D.Affectation_trans_dimension(epsBB_tdt_3_inter,false);
|
|
Tenseur3BB delta_epsBB_3_inter(delta_epsBB); delta_epsBB_3_inter.ChBase(beta);
|
|
delta_eps_BB_1D.Affectation_trans_dimension(delta_epsBB_3_inter,false);
|
|
};
|
|
// cas des jacobiens: il nous faut définir des jacobiens cohérents en 1D au niveau de leur
|
|
// variation entre 0 et tdt
|
|
// on se sert de la déformation et d'un jacobien final mis à 1. arbitrairement
|
|
jacobien_tdt_1D = 1.; // initialisation arbitraire
|
|
// -- maintenant on tient compte du type de déformation, et de sa valeur 11
|
|
// sachant que le repère actuel est normé, donc le repère initial ne peut pas être normé
|
|
// le vecteur de base est le vecteur propre
|
|
// en 1D la longueur est le jacobien
|
|
switch (type_de_deformation)
|
|
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
|
|
// cas d'une déformation d'Almansi
|
|
{ // epsBB11 = 1/2 * (1. - (l0/l)^2), dans un repère normé, ce qui est le cas ici
|
|
// -> l0 = l * sqrt(1. - 2.* eps11)
|
|
jacobien_0_1D = sqrt(1.-eps_BB_1D(1,1));
|
|
};
|
|
break;
|
|
case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
|
|
// cas d'une def logarithmique ou une approximation
|
|
// epsBB11 = log (l/l0) --> l0 = l / exp(epsBB11)
|
|
{ jacobien_0_1D = jacobien_tdt_1D / exp(eps_BB_1D(1,1));
|
|
};
|
|
break;
|
|
default :
|
|
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
|
|
<< Nom_type_deformation(type_de_deformation);
|
|
cout << "\n LoiCritere::Passage_3ou2_vers_1(.. \n";
|
|
Sortie(1);
|
|
};
|
|
|
|
};
|
|
|
|
// passage inverse: chgt de repère et de dimension pour les variables de passage
|
|
// et stockage des infos de doublement plane pour le post-traitement pour le prochain appel
|
|
// en entrée: d_sigma_deps_1D : l'opérateur qui a été calculé
|
|
// en sortie: d_sigma_deps_inter
|
|
void LoiCritere::Passage_1_vers_3ou2(const Mat_pleine& gamma,TenseurHH & sigHH
|
|
,const TenseurHHHH& d_sigma_deps_1D
|
|
,const Mat_pleine& beta
|
|
,TenseurHHHH& d_sigma_deps_inter,const Met_abstraite::Umat_cont& ex
|
|
,const Tableau <Coordonnee2H>& V_Pr_H
|
|
,const Tableau <Coordonnee>& V_P_sig)
|
|
{ // on s'occupe du changement de base
|
|
int dim_tens = abs(sigHH.Dimension());
|
|
// if (dim_tens == 2)
|
|
// { // normalement les contraintes sont nulles dans les directions autres que 11
|
|
// sigHH.Affectation_trans_dimension(sig_HH_t_1D,true);
|
|
// sigHH.ChBase(beta); // on revient dans le repère naturel
|
|
// }
|
|
// else // cas 3D
|
|
// { Tenseur3HH sigHH_3_inter(sigHH_t); sigHH_3_inter.ChBase(gamma);
|
|
// sig_HH_t_1D.Affectation_trans_dimension(sigHH_3_inter,false);// passage des tenseurs résultats à l'ordre 1
|
|
// Tenseur3BB DepsBB_3_inter(DepsBB); DepsBB_3_inter.ChBase(gamma);
|
|
// Deps_BB_1D.Affectation_trans_dimension(DepsBB_3_inter,false);
|
|
// Tenseur3BB epsBB_tdt_3_inter(epsBB_tdt);epsBB_tdt_3_inter.ChBase(gamma);
|
|
// eps_BB_1D.Affectation_trans_dimension(epsBB_tdt_3_inter,false);
|
|
// Tenseur3BB delta_epsBB_3_inter(delta_epsBB); delta_epsBB_3_inter.ChBase(gamma);
|
|
// delta_eps_BB_1D.Affectation_trans_dimension(delta_epsBB_3_inter,false);
|
|
// };
|
|
|
|
// normalement les contraintes sont nulles dans les directions autres que 11
|
|
sigHH.Affectation_trans_dimension(sig_HH_1D,true);
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{Tenseur2HH sigHH_2_inter(sigHH);
|
|
cout << "\n LoiCritere::Passage_1_vers_3ou2(";
|
|
cout << "\n sigHH (en VP) affecté du resultat 1D = "; sigHH.Ecriture(cout);
|
|
// en base orthonormee
|
|
Tenseur3HH sigHH_3_inter;
|
|
sig_HH_1D.BaseAbsolue(sigHH_3_inter,(giB_tdt_1D));
|
|
cout << "\n sigHH 1D =en orthonormee "<<sigHH_3_inter;
|
|
// cout << endl;
|
|
};
|
|
#endif
|
|
|
|
// if ((ParaGlob::NiveauImpression() > 8) || (Permet_affichage() > 5))
|
|
// cout << "\n test2: LoiCritere::Passage_1_vers_3ou2(... " << flush;
|
|
// changement de base pour la contrainte
|
|
// ici il s'agit de revenir à la base initiale il faut donc utiliser
|
|
// l'inverse gamma c-a-d finalement beta transposée
|
|
sigHH.ChBase(beta.Transpose()); // on revient dans le repère naturel
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{Tenseur2HH sigHH_2_inter(sigHH);
|
|
// cout << "\n debug LoiCritere::Passage_1_vers_3ou2(";
|
|
cout << "\n sigHH en g_i après changement de base "; sigHH.Ecriture(cout);
|
|
// en base orthonormee
|
|
Tenseur3HH sigHH_3_inter;
|
|
sigHH.BaseAbsolue(sigHH_3_inter,(* ex.giB_tdt));
|
|
cout << "\n sigHH = en orthonormee "<<sigHH_3_inter;
|
|
// cout << endl;
|
|
};
|
|
#endif
|
|
|
|
// concernant l'opérateur tangent: pour l'instant on a un opérateur en 33, mais une seule
|
|
// valeur est utilisée ! peut-être à changer par la suite, en tenseur 11
|
|
// Tenseur3HHHH& d_sig_deps_3HHHH = *((Tenseur3HHHH *) &d_sigma_deps);
|
|
// d_sig_deps_1HHHH.ChangeBase(<#TenseurHHHH &A#>, const BaseB &gi)
|
|
if (dim_tens == 2)
|
|
{ // un tenseur de travail qui contient l'opérateur tangent, mais transporté en 2D
|
|
// if ((ParaGlob::NiveauImpression() > 8) || (Permet_affichage() > 5))
|
|
// cout << "\n test3: LoiCritere::Passage_1_vers_3ou2(... " << flush;
|
|
Tenseur2HHHH tens_travail2HHHH;
|
|
bool plusZero=true;
|
|
tens_travail2HHHH.Affectation_trans_dimension(*d_sigma_deps_1D_P,plusZero);
|
|
// construction de la base qui permet de passer des vecteurs propres au gi
|
|
// donc on veut: V_al = beta_al^be g_be
|
|
// il se trouve que c'est égal à la matrice passée en paramètre : beta(al,be)
|
|
// if ((ParaGlob::NiveauImpression() > 8) || (Permet_affichage() > 5))
|
|
// cout << "\n test4: LoiCritere::Passage_1_vers_3ou2(... " << flush;
|
|
BaseB V_al_B(2,2,0.);
|
|
V_al_B.CoordoB(1)(1) = beta(1,1);V_al_B.CoordoB(1)(2) = beta(1,2);
|
|
V_al_B.CoordoB(2)(1) = beta(2,1);V_al_B.CoordoB(2)(2) = beta(2,2);
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{//cout << "\n debug LoiCritere::Passage_1_vers_3ou2(";
|
|
// vérif du calcul de la base
|
|
// cout << "\n BaseB des gi \n "<< (* ex.giB_tdt);
|
|
cout << "\n construction d'une BaseB V_al_B pour cgt base tenseurHHHH \n "<< V_al_B;
|
|
// calcul long
|
|
BaseB V_inter_B(2,2,0.);
|
|
for (int al = 1;al < 3; al++)
|
|
for (int be = 1; be < 3; be++)
|
|
V_inter_B.CoordoB(al)(be) = (* ex.giH_tdt).Coordo(be) * V_P_sig(al);
|
|
cout << "\n BaseB V_inter_B "<< V_inter_B;
|
|
|
|
// pour voir si la base est correcte on la teste sur le tenseur des contraintes
|
|
Tenseur2HH sigHH_2_inter;
|
|
// on exprime sigHH dans la base des V_P qui ici est considéré comme absolue pour la vérif
|
|
sig_HH_1D.BaseAbsolue(sigHH_2_inter,V_al_B);
|
|
cout << "\n verif via sig_HH_1D = dans la base des gi de nouveau \n "<<sigHH_2_inter;
|
|
// cout << endl;
|
|
};
|
|
#endif
|
|
|
|
// on fait le changement de base
|
|
Tenseur2HHHH* d_sig_deps_2HHHH = ((Tenseur2HHHH *) &d_sigma_deps_inter);
|
|
// d_sig_deps_2HHHH->ChangeBase(tens_travail2HHHH,gi_al_B);
|
|
tens_travail2HHHH.ChangeBase(*d_sig_deps_2HHHH,V_al_B);
|
|
}
|
|
else // cas 3D
|
|
{ //if (d_sigma_deps_inter == NULL)
|
|
// d_sigma_deps_inter = NevezTenseurHHHH(3);
|
|
// la suite dépendra de l'application donc du critère, donc on met un message
|
|
// pour l'instant
|
|
cout << "\n pour l'instant la suite n'est pas codee !! se plaindre !! "
|
|
<< "\n LoiCritere::Passage_1_vers_3ou2(..." << endl;
|
|
Sortie(2);
|
|
};
|
|
// if ((ParaGlob::NiveauImpression() > 8) || (Permet_affichage() > 5))
|
|
// cout << "\n test7: LoiCritere::Passage_1_vers_3ou2(... " << flush;
|
|
|
|
|
|
};
|
|
|
|
// fonction interne utilisée par les classes dérivées de Loi_comp_abstraite
|
|
// pour répercuter les modifications de la température
|
|
// ici utiliser pour modifier la température des lois élémentaires
|
|
// l'Enum_dure: indique quel est la température courante : 0 t ou tdt
|
|
void LoiCritere::RepercuteChangeTemperature(Enum_dure temps)
|
|
{ list <Loi_comp_abstraite *>::const_iterator ili,ilifin=lois_internes.end();
|
|
for (ili=lois_internes.begin();ili!=ilifin;ili++)
|
|
{(*ili)->temperature_0 = this->temperature_0;
|
|
(*ili)->temperature_t = this->temperature_t;
|
|
(*ili)->temperature_tdt = this->temperature_tdt;
|
|
(*ili)->dilatation=dilatation;
|
|
|
|
// on répercute également les déformations thermiques, qui ne sont utilisées
|
|
// telles quelles que pour certaines lois: ex: loi hyper-élastique
|
|
if (dilatation)
|
|
{// a- dimensionnement des tenseurs intermédiaires
|
|
int dim_tens = epsBB_therm->Dimension();
|
|
// -- cas de la déformation
|
|
if ((*ili)->epsBB_therm == NULL) { (*ili)->epsBB_therm = NevezTenseurBB(dim_tens);}
|
|
else if ((*ili)->epsBB_therm->Dimension() != dim_tens)
|
|
{ delete (*ili)->epsBB_therm;(*ili)->epsBB_therm = NevezTenseurBB(dim_tens);};
|
|
// -- cas de la vitesse de déformation
|
|
if ((*ili)->DepsBB_therm == NULL) { (*ili)->DepsBB_therm = NevezTenseurBB(dim_tens);}
|
|
else if ((*ili)->DepsBB_therm->Dimension() != dim_tens)
|
|
{ delete (*ili)->DepsBB_therm;(*ili)->DepsBB_totale = NevezTenseurBB(dim_tens);};
|
|
// b- affectation des tenseurs
|
|
(*(*ili)->epsBB_therm)=(*epsBB_therm);
|
|
(*(*ili)->DepsBB_therm)=(*DepsBB_therm);
|
|
};
|
|
|
|
(*ili)->RepercuteChangeTemperature(temps);
|
|
|
|
switch (temps)
|
|
{ case TEMPS_0: {(*ili)->temperature = &(*ili)->temperature_0; break;}
|
|
case TEMPS_t: {(*ili)->temperature = &(*ili)->temperature_t; break;}
|
|
case TEMPS_tdt: {(*ili)->temperature = &(*ili)->temperature_tdt; break;}
|
|
default:
|
|
{ cout << "\n erreur, cas de temps non prevu !! "
|
|
<< "\n LoiCritere::RepercuteChangeTemperature(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
};
|
|
};
|
|
|
|
// fonction surchargée dans les classes dérivée si besoin est
|
|
void LoiCritere::CalculGrandeurTravail
|
|
(const PtIntegMecaInterne& ptintmeca
|
|
,const Deformation & def,Enum_dure temps,const ThermoDonnee& dTP
|
|
,const Met_abstraite::Impli* ex_impli
|
|
,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
|
|
,const Met_abstraite::Umat_cont* ex_umat
|
|
,const List_io<Ddl_etendu>* exclure_dd_etend
|
|
,const List_io<const TypeQuelconque *>* exclure_Q
|
|
)
|
|
{ if (avec_ponderation)
|
|
{ // on passe en revue les grandeurs servant au calcul des fonctions de proportion
|
|
list <Ponderation >::iterator il,ilfin=list_ponderation.end();
|
|
list <double>::iterator ipfonc = fonc_ponder.begin();
|
|
for (il=list_ponderation.begin();il != ilfin;il++,ipfonc++)
|
|
{ Ponderation& ponder = (*il); // pour simplifier
|
|
(*ipfonc)=1.; // init
|
|
int taille = ponder.Type_grandeur().Taille();
|
|
for (int i=1;i<= taille; i++)
|
|
{ // calcul des pondérations
|
|
if (ponder.Valeur_aux_noeuds()(i)) // pour l'instant que des ddl patentés !
|
|
// cas d'une proportion provenant d'une interpolation aux noeuds
|
|
{ double grand = def.DonneeInterpoleeScalaire(ponder.Type_grandeur()(i).Enum(),temps);
|
|
double fonc = ponder.C_proport()(i)->Valeur(grand);
|
|
(*ipfonc) *= fonc; // on accumule multiplicativement
|
|
}
|
|
else
|
|
// sinon il s'agit d'une grandeur directement accessible au point d'intégration
|
|
// pour l'instant il n'y a pas de procédure générale de récupération, seulement des cas particuliers
|
|
{
|
|
// deux cas suivant que l'on a affaire à un ddl de base ou à un vrai ddl étendu
|
|
if (ponder.Type_grandeur()(i).Nom_vide())
|
|
{switch (ponder.Type_grandeur()(i).Enum())
|
|
{ case PROP_CRISTA:
|
|
{ const double* taux_crita = dTP.TauxCrista();
|
|
if (taux_crita != NULL)
|
|
{ double fonc = ponder.C_proport()(i)->Valeur(*taux_crita);
|
|
(*ipfonc) *= fonc; // on accumule multiplicativement
|
|
}
|
|
else
|
|
{cout << "\n erreur, le taux de cristalinite n'est pas disponible au point d'integration "
|
|
<< " il n'est pas possible de calculer la proportion pour la loi des melanges "
|
|
<< "\n LoiCritere::CalculGrandeurTravail(.... ";
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur, le type de proportion " << ponder.Type_grandeur()(i) << " n'est pas disponible "
|
|
<< " pour l'instant au point d'integration ! "
|
|
<< "\n LoiCritere::CalculGrandeurTravail(.... ";
|
|
};
|
|
}
|
|
else
|
|
{// cas d'un vrai ddl étendue
|
|
switch (ponder.Type_grandeur()(i).Position()-NbEnum_ddl())
|
|
{case 87: // cas de "def_equivalente"
|
|
{const double def_equivalente = ptintmeca.Deformation_equi_const()(1); // recup de la def equi
|
|
double fonc = ponder.C_proport()(i)->Valeur(def_equivalente);
|
|
(*ipfonc) *= fonc; // on accumule multiplicativement
|
|
//cout << "\n defEqui= " << def_equivalente << " fonct= " << fonc << " ";
|
|
break;
|
|
}
|
|
case 88: // cas de "def_duale_mises_maxi"
|
|
{const double def_duale_mises_maxi = ptintmeca.Deformation_equi_const()(3); // recup de la def equi
|
|
double fonc = ponder.C_proport()(i)->Valeur(def_duale_mises_maxi);
|
|
(*ipfonc) *= fonc; // on accumule multiplicativement
|
|
break;
|
|
}
|
|
case 89: // cas de "vitesse_def_equivalente"
|
|
{const double delta_def_equivalente = ptintmeca.Deformation_equi_const()(4); // recup du delta def equi
|
|
// recup de l'incrément de temps
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
|
|
double unSurDeltat=0;
|
|
if (Abs(deltat) >= ConstMath::trespetit)
|
|
{unSurDeltat = 1./deltat;}
|
|
else // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
|
|
{ // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
if (unSurDeltat < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "; };
|
|
unSurDeltat = ConstMath::tresgrand;
|
|
};
|
|
double vitesse_def_equi = delta_def_equivalente * unSurDeltat;
|
|
double fonc = ponder.C_proport()(i)->Valeur(vitesse_def_equi);
|
|
(*ipfonc) *= fonc; // on accumule multiplicativement
|
|
//cout << "\n vit= " << vitesse_def_equi << " f= " << fonc << endl;
|
|
break;
|
|
}
|
|
case 77: // cas de "def_duale_mises"
|
|
{const double def_duale_mises = ptintmeca.Deformation_equi_const()(2); // recup de la def equi
|
|
double fonc = ponder.C_proport()(i)->Valeur(def_duale_mises);
|
|
(*ipfonc) *= fonc; // on accumule multiplicativement
|
|
break;
|
|
}
|
|
case 78: // cas de "Spherique_eps"
|
|
{const Vecteur & epsInvar = ptintmeca.EpsInvar_const(); // recup des invariants
|
|
double spherique_eps = epsInvar(1);
|
|
double fonc = ponder.C_proport()(i)->Valeur(spherique_eps);
|
|
(*ipfonc) *= fonc; // on accumule multiplicativement
|
|
break;
|
|
}
|
|
case 81: // cas de "Spherique_sig"
|
|
{const Vecteur & sigInvar = ptintmeca.SigInvar_const(); // recup des invariants
|
|
double spherique_sig = sigInvar(1);
|
|
double fonc = ponder.C_proport()(i)->Valeur(spherique_sig);
|
|
(*ipfonc) *= fonc; // on accumule multiplicativement
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur, le type de proportion " << ponder.Type_grandeur()(i)
|
|
<< " n'est pas disponible "
|
|
<< " pour l'instant au point d'integration ! "
|
|
<< "\n LoiCritere::CalculGrandeurTravail(.... ";
|
|
break;
|
|
};
|
|
};
|
|
};
|
|
};
|
|
};
|
|
};
|
|
// enregistrement dans les variables spécifiques au point calculé
|
|
SaveResul_LoiCritere & save_resul = *((SaveResul_LoiCritere*) saveResul);
|
|
save_resul.f_ponder = fonc_ponder;
|
|
// on initialise le niveau de déclanchement avec les grandeurs de la lois
|
|
// si la valeur est différente de 0, sinon ce sera 1.
|
|
// if (niveau_declenchement_critere != 0.)
|
|
// {save_resul.niveau_declenchement_critere = niveau_declenchement_critere;}
|
|
// else
|
|
// {save_resul.niveau_declenchement_critere = 1.;};
|
|
// init par défaut
|
|
save_resul.niveau_declenchement_critere = niveau_declenchement_critere;
|
|
|
|
if (avecNiveauSigmaI_mini_pour_plis)
|
|
{// on initialise le niveau de déclanchement avec les grandeurs de la lois
|
|
// si la valeur est différente de 0, sinon ce sera 1.
|
|
// car en multiplication, si on par de 0 cela donnera toujours 0
|
|
// if (niveau_declenchement_critere == 0.)
|
|
// {save_resul.niveau_declenchement_critere = 1.;};
|
|
// non on ne fait pas cela sinon c'est compliqué à comprendre
|
|
// donc si on est dans un niveau réglable, systématiquement on met 1.
|
|
// ce qui fait que l'on a directement la valeur de la fonction
|
|
save_resul.niveau_declenchement_critere = 1.;
|
|
|
|
if (niveauF_fct_nD != NULL)
|
|
{ // on appel la fonction générique
|
|
list <SaveResul*> list_save; // inter pour l'appel de la fonction
|
|
list_save.push_back(saveResul);
|
|
|
|
Tableau <double> & tab = Loi_comp_Valeur_FnD_Evoluee
|
|
(niveauF_fct_nD->C_proport(),1
|
|
,ex_impli,ex_expli_tdt,ex_umat
|
|
,exclure_dd_etend,exclure_Q
|
|
,&list_save
|
|
);
|
|
|
|
// // comme il s'agit d'une grandeur globale, on n'a pas besoin de la récupérer
|
|
// double fonc = niveauF_fct_nD->C_proport()->Valeur_pour_variables_globales()(1);
|
|
// save_resul.niveau_declenchement_critere *= fonc;
|
|
|
|
save_resul.niveau_declenchement_critere *= tab(1);
|
|
|
|
};
|
|
if (niveauF_temps != NULL) // cas d'une dépendance au temps
|
|
{ const VariablesTemps& var_temps = ParaGlob::Variables_de_temps();
|
|
double fonc = niveauF_temps->Valeur(var_temps.TempsCourant());
|
|
save_resul.niveau_declenchement_critere *= fonc;
|
|
};
|
|
// avec une dépendance via éventuellement un ddl étendu
|
|
if (niveauF_ddlEtendu != NULL)
|
|
{ Ponderation& ponder = *niveauF_ddlEtendu; // pour simplifier
|
|
double fonc_finale = 1.; // init
|
|
int taille = ponder.Type_grandeur().Taille();
|
|
for (int i=1;i<= taille; i++)
|
|
{ // calcul des pondérations
|
|
if (ponder.Valeur_aux_noeuds()(i)) // pour l'instant que des ddl patentés !
|
|
// cas d'une proportion provenant d'une interpolation aux noeuds
|
|
{ double grand = def.DonneeInterpoleeScalaire(ponder.Type_grandeur()(i).Enum(),temps);
|
|
double fonc = ponder.C_proport()(i)->Valeur(grand);
|
|
fonc_finale *= fonc; // on accumule multiplicativement
|
|
}
|
|
else
|
|
// sinon il s'agit d'une grandeur directement accessible au point d'intégration
|
|
// pour l'instant il n'y a pas de procédure générale de récupération, seulement des cas particuliers
|
|
{
|
|
// deux cas suivant que l'on a affaire à un ddl de base ou à un vrai ddl étendu
|
|
if (ponder.Type_grandeur()(i).Nom_vide())
|
|
{switch (ponder.Type_grandeur()(i).Enum())
|
|
{ case PROP_CRISTA:
|
|
{ const double* taux_crita = dTP.TauxCrista();
|
|
if (taux_crita != NULL)
|
|
{ double fonc = ponder.C_proport()(i)->Valeur(*taux_crita);
|
|
fonc_finale *= fonc; // on accumule multiplicativement
|
|
}
|
|
else
|
|
{cout << "\n erreur, le taux de cristalinite n'est pas disponible au point d'integration "
|
|
<< " il n'est pas possible de calculer le niveau_declenchement_critere "
|
|
<< "\n LoiCritere::CalculGrandeurTravail(.... ";
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur, le type " << ponder.Type_grandeur()(i) << " n'est pas disponible "
|
|
<< " pour l'instant au point d'integration ! "
|
|
<< "\n LoiCritere::CalculGrandeurTravail(.... ";
|
|
};
|
|
}
|
|
else
|
|
{// cas d'un vrai ddl étendue
|
|
switch (ponder.Type_grandeur()(i).Position()-NbEnum_ddl())
|
|
{case 87: // cas de "def_equivalente"
|
|
{const double def_equivalente = ptintmeca.Deformation_equi_const()(1); // recup de la def equi
|
|
double fonc = ponder.C_proport()(i)->Valeur(def_equivalente);
|
|
fonc_finale *= fonc; // on accumule multiplicativement
|
|
break;
|
|
}
|
|
case 88: // cas de "def_duale_mises_maxi"
|
|
{const double def_duale_mises_maxi = ptintmeca.Deformation_equi_const()(3); // recup de la def equi
|
|
double fonc = ponder.C_proport()(i)->Valeur(def_duale_mises_maxi);
|
|
fonc_finale *= fonc; // on accumule multiplicativement
|
|
break;
|
|
}
|
|
case 89: // cas de "vitesse_def_equivalente"
|
|
{const double delta_def_equivalente = ptintmeca.Deformation_equi_const()(4); // recup du delta def equi
|
|
// recup de l'incrément de temps
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
|
|
double unSurDeltat=0;
|
|
if (Abs(deltat) >= ConstMath::trespetit)
|
|
{unSurDeltat = 1./deltat;}
|
|
else // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
|
|
{ // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
if (unSurDeltat < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "; };
|
|
unSurDeltat = ConstMath::tresgrand;
|
|
};
|
|
double vitesse_def_equi = delta_def_equivalente * unSurDeltat;
|
|
double fonc = ponder.C_proport()(i)->Valeur(vitesse_def_equi);
|
|
fonc_finale *= fonc; // on accumule multiplicativement
|
|
break;
|
|
}
|
|
case 77: // cas de "def_duale_mises"
|
|
{const double def_duale_mises = ptintmeca.Deformation_equi_const()(2); // recup de la def equi
|
|
double fonc = ponder.C_proport()(i)->Valeur(def_duale_mises);
|
|
fonc_finale *= fonc; // on accumule multiplicativement
|
|
break;
|
|
}
|
|
case 78: // cas de "Spherique_eps"
|
|
{const Vecteur & epsInvar = ptintmeca.EpsInvar_const(); // recup des invariants
|
|
double spherique_eps = epsInvar(1);
|
|
double fonc = ponder.C_proport()(i)->Valeur(spherique_eps);
|
|
fonc_finale *= fonc; // on accumule multiplicativement
|
|
break;
|
|
}
|
|
case 81: // cas de "Spherique_sig"
|
|
{const Vecteur & sigInvar = ptintmeca.SigInvar_const(); // recup des invariants
|
|
double spherique_sig = sigInvar(1);
|
|
double fonc = ponder.C_proport()(i)->Valeur(spherique_sig);
|
|
fonc_finale *= fonc; // on accumule multiplicativement
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur, le type " << ponder.Type_grandeur()(i)
|
|
<< " n'est pas disponible "
|
|
<< " pour l'instant au point d'integration ! "
|
|
<< "\n LoiCritere::CalculGrandeurTravail(.... ";
|
|
break;
|
|
};
|
|
};
|
|
};
|
|
};
|
|
// on met à jour le niveau
|
|
save_resul.niveau_declenchement_critere *= fonc_finale;
|
|
};
|
|
};
|
|
|
|
// répercution sur les classes dérivées si besoin est
|
|
list <SaveResul*>::iterator isave=save_resul.liste_des_SaveResul.begin(); // pour les saveResul des lois
|
|
list <Loi_comp_abstraite *>::const_iterator ili,ilifin=lois_internes.end();
|
|
for (ili=lois_internes.begin();ili!=ilifin;ili++,isave++)
|
|
{// passage des informations spécifique à la loi liste_des_SaveResul
|
|
(*ili)->IndiqueSaveResult(*isave);
|
|
(*ili)->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
|
|
(*ili)->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours
|
|
(*ili)->CalculGrandeurTravail(ptintmeca,def,temps,dTP
|
|
,ex_impli,ex_expli_tdt,ex_umat,exclure_dd_etend,exclure_Q);
|
|
};
|
|
};
|
|
|
|
// fonction critère de plissement de biel
|
|
// en entrée:
|
|
// implicit : si oui on est en implicite, sinon en explicite
|
|
// retour:
|
|
// 0 : il y a eu un pb que l'on n'a pas peu résoudre, rien n'a été modifié
|
|
// 1 : le critère n'a rien modifié
|
|
// 2 : l'application du critère conduit à une contrainte et un opérateur tangent nul
|
|
int LoiCritere::Critere_plis_biel(bool en_base_orthonormee,TenseurHH & sigHH_t,TenseurBB& DepsBB
|
|
,TenseurBB & epsBB_tdt,TenseurBB & delta_epsBB,double& jacobien_0,double& jacobien
|
|
,TenseurHH& sigHH,TenseurHHHH& d_sigma_deps_inter
|
|
,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite
|
|
,double& module_cisaillement
|
|
,bool implicit,const Met_abstraite::Umat_cont& ex)
|
|
{// dim et vérification
|
|
int dim_tens = abs(sigHH.Dimension());
|
|
// normalement le critère ne s'applique que pour des tenseurs 1D
|
|
#ifdef MISE_AU_POINT
|
|
if (dim_tens != 1)
|
|
{ cout << "\n erreur le critere de plissement de biel ne peut s'appliquer que pour des lois 1D"
|
|
<< " \n LoiCritere::Critere_plis_biel(..";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
SaveResul_LoiCritere & save_resul = *((SaveResul_LoiCritere*) saveResul);
|
|
|
|
// on va utiliser la valeur de la contrainte, par contre on la ramène dans un repère normé
|
|
// pour son intensité ce qui revient à calculer la trace du tenseur
|
|
double sig_I = sigHH && (* ex.gijBB_tdt);
|
|
int retour = 0; // init par défaut
|
|
if (sig_I >= niveau_declenchement_critere) // pas de plis
|
|
{ // on se trouve dans le cas où rien n'est en compression, donc cas normal sans plissement
|
|
retour = 1; //on signale
|
|
save_resul.eps_pli.Zero(); // pas de pli répertorié
|
|
}
|
|
else
|
|
{ // on se trouve dans le cas où on a des plis
|
|
|
|
// l'idée est de ce positionner au niveau du seuil de déclenchement, pour cela on va utiliser
|
|
// un facteur multiplicatif de manière à pouvoir conserver un opérateur tangent cohérent avec le seuil
|
|
// (si on impose dans une boucle de Newton, la contrainte on va récupérer une def méca et un opérateur
|
|
// tangent par rapport à cette def méca. Du coup on ne pourra pas en déduire un opérateur tangent par
|
|
// rapport à la def cinématique ...)
|
|
// double rapport = niveau_declenchement_critere / (Dabs(sig_I)+ConstMath::pasmalpetit);
|
|
|
|
|
|
|
|
|
|
// on annule la contrainte et sa contribution éventuelle à la raideur
|
|
sigHH.Inita(0.);
|
|
d_sigma_deps_inter.Inita(0.); // dans le cas explicite, ne sert à rien mais c'est rapide
|
|
// comme marqueur de l'intensité des plis on utilise la def de compression
|
|
save_resul.eps_pli(1) = epsBB_tdt && (* ex.gijHH_tdt);
|
|
// d'autre part on définit le vecteur propre normé qui n'est autre que gi_1 normé
|
|
int dim = ParaGlob::Dimension();
|
|
if (save_resul.V_P_sig == NULL)
|
|
save_resul.V_P_sig = new Tableau <Coordonnee>;
|
|
Coordonnee zero(dim); // un vecteur nul
|
|
save_resul.V_P_sig->Change_taille(2); // on n'a besoin que d'un vecteur
|
|
// **** mais pour l'instant j'en mets 2 de manière à être cohérent avec les membranes
|
|
// sinon pb en sortie des grandeurs
|
|
Coordonnee inter = ex.giB_tdt->Coordo(1);
|
|
(*save_resul.V_P_sig)(1)= inter.Normer();
|
|
retour = 2;
|
|
};
|
|
|
|
return retour;
|
|
};
|
|
|
|
|
|
// récupération de la variation relative d'épaisseur calculée: h/h0
|
|
// cette variation n'est utile que pour des lois en contraintes planes
|
|
// - pour les lois 3D : retour d'un nombre très grand, indiquant que cette fonction est invalide
|
|
// - pour les lois 2D def planes: retour de 0
|
|
// les infos nécessaires à la récupération , sont stockées dans saveResul
|
|
// qui est le conteneur spécifique au point où a été calculé la loi
|
|
double LoiCritere::HsurH0(SaveResul * saveResul) const
|
|
{ SaveResul_LoiCritere & save_resul = *((SaveResul_LoiCritere*) saveResul);
|
|
switch (type_critere)
|
|
{ case PLISSEMENT_MEMBRANE:
|
|
{ // traitement des épaisseurs
|
|
LoiContraintesPlanes::SaveResul_LoiContraintesPlanes * save_result_plan_inter // pour simplifier
|
|
= (LoiContraintesPlanes::SaveResul_LoiContraintesPlanes *)
|
|
(*(save_resul.liste_des_SaveResul.begin()));
|
|
return save_result_plan_inter->hsurh0;
|
|
break;
|
|
}
|
|
case PLISSEMENT_BIEL:
|
|
{ LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble * save_result_1DCP2_inter // pour simplifier
|
|
= (LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble *)
|
|
(*(save_resul.liste_des_SaveResul.begin()));
|
|
return save_result_1DCP2_inter->hsurh0;
|
|
break;
|
|
}
|
|
default :
|
|
cout << "\nErreur : LoiCritere::HsurH0(.. non implante pour le critere " << Nom_Critere_Loi(type_critere)
|
|
<< " !\n";
|
|
Sortie(1);
|
|
};
|
|
|
|
|
|
|
|
};
|
|
|
|
|
|
// premier type de calcul dans le cas d'un pli dans une seule direction
|
|
void LoiCritere::Premie_type_calcul_en_un_pli(const TenseurBB & epsBB_tdt,const TenseurBB & delta_epsBB
|
|
,const TenseurHH & sigHH_t
|
|
,const double& jacobien_0,const double& jacobien
|
|
,EnergieMeca & energ,const EnergieMeca & energ_t
|
|
,const TenseurBB& DepsBB
|
|
,double& module_compressibilite,double& module_cisaillement
|
|
,const Tableau <Coordonnee2H>& V_Pr_H
|
|
,TenseurHH& sigHH,TenseurHHHH& d_sigma_deps_inter
|
|
,const Coordonnee2& valPropreSig
|
|
,const Met_abstraite::Umat_cont& ex)
|
|
{
|
|
SaveResul_LoiCritere & save_resul = *((SaveResul_LoiCritere*) saveResul);
|
|
// les coordonnées des vecteurs propres sont exprimées dans l'ancienne base
|
|
// on a donc
|
|
// Vi_B(i) = gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
|
|
|
|
Mat_pleine beta(2,2);
|
|
for (int i=1;i<3;i++) // < 3 donc de 1 à 2
|
|
for (int j=1;j<3;j++)
|
|
beta(i,j) = V_Pr_H(i)(j);
|
|
Mat_pleine gamma(beta.Inverse().Transpose()); // on calcule la matrice inverse
|
|
// la matrice gamma correspond à la matrice de changement de base des g^i au gpH
|
|
// gamma est l'inverse transposée
|
|
// gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
|
|
// c-a-d= gp^i = gamma^i_j * g^j
|
|
|
|
// on va également définir une matrice de changement pour les tenseurs 3D
|
|
Tableau <Coordonnee3H> V_Pr3D_H(3);
|
|
// comme V_Pr est de dim 2, il faut transférer les coordonnées
|
|
V_Pr3D_H(1)(1) = V_Pr_H(1)(1);V_Pr3D_H(1)(2) = V_Pr_H(1)(2);
|
|
V_Pr3D_H(2)(1) = V_Pr_H(2)(1);V_Pr3D_H(2)(2) = V_Pr_H(2)(2);
|
|
// le troisième est en local 0,0,1
|
|
V_Pr3D_H(3)(3) = 1.;
|
|
Mat_pleine beta3D(3,3);// initialisée à 0
|
|
// d'abord les 2 vecteurs propres
|
|
for (int i=1;i<3;i++)
|
|
for (int j=1;j<4;j++)
|
|
beta3D(i,j) = V_Pr3D_H(i)(j);
|
|
// puis la normale
|
|
beta3D(3,3) = 1.;
|
|
Mat_pleine gamma3D(beta3D.Inverse().Transpose());
|
|
|
|
|
|
#ifdef MISE_AU_POINT
|
|
// pour tester, on passe les contraintes dans le repère principal
|
|
if (Permet_affichage() > 5)
|
|
{Tenseur2HH sigHH_2_inter(sigHH);
|
|
cout << "\n utilisation du repere ortho des vecteurs propres: V_P ";
|
|
cout << "\n sigHH= en giB : "; sigHH.Ecriture(cout);
|
|
cout << "\n base giB(1) "<< (* ex.giB_tdt)(1)
|
|
<< "\n base giB(2) "<< (* ex.giB_tdt)(2);
|
|
// en base orthonormee
|
|
Tenseur3HH sigHH_3_inter;
|
|
sigHH.BaseAbsolue(sigHH_3_inter,(* ex.giB_tdt));
|
|
// TenseurHH* ptHH = NevezTenseurHH((*(*isig)));
|
|
// (*(*isig)).BaseAbsolue(*ptHH,giB);
|
|
cout << "\n sigHH= en absolue 3D : "<<sigHH_3_inter;
|
|
cout << "\n chgt de base: Vi_B(i) = gpB(i) = beta(i,j) * gB(j) --> beta= ";beta.Affiche();
|
|
cout << "\n sigHH_2_inter avant chgt= "<< sigHH_2_inter;
|
|
sigHH_2_inter.ChBase(gamma,false);
|
|
cout << "\n sigHH= en V_P : "<< sigHH_2_inter;
|
|
// on revient à la base de départ pour vérifier que l'opération inverse ne pose
|
|
// pas de pb
|
|
sigHH_2_inter.ChBase(beta.Transpose(),false);
|
|
cout << "\n verif: sigHH après retour au repère initial giB: ";
|
|
sigHH_2_inter.Ecriture(cout);
|
|
//cout << "\n gamma * beta = "; (gamma * beta).Affiche();
|
|
cout << "\n vecteurs propres en 3D ortho: "; //<< (*save_resul.V_P_sig);
|
|
cout << "\n VP1 "<< (*save_resul.V_P_sig)(1);
|
|
cout << "\n VP2 "<< (*save_resul.V_P_sig)(2);
|
|
cout << "\n VP3 "<< (*save_resul.V_P_sig)(3);
|
|
// vérif des vecteurs propres
|
|
Coordonnee3B VpB1; VpB1.Change_val((*save_resul.V_P_sig)(1));
|
|
Coordonnee3H VpropreH1 = sigHH_3_inter * VpB1;
|
|
//cout << "\n sigHH_3_inter * VP1 "<< VpropreH1;
|
|
VpropreH1 /= (valPropreSig(1)+ConstMath::petit);
|
|
cout << "\n (sigHH_3_ortho * VP1) / (lambda1+epsilon)= "<< VpropreH1;
|
|
Coordonnee3B VpB2; VpB2.Change_val((*save_resul.V_P_sig)(2));
|
|
Coordonnee3H VpropreH2 = sigHH_3_inter * VpB2;
|
|
//cout << "\n sigHH_3_inter * VP2 "<< VpropreH2;
|
|
VpropreH2 /= (valPropreSig(2)+ConstMath::petit);
|
|
cout << "\n (sigHH_3_ortho * VP2) / (lambda2+epsilon)= "<< VpropreH2;
|
|
};
|
|
#endif
|
|
|
|
{ // --il faut créer un save_result pour l'appel des contraintes doublement planes
|
|
// on récupère le conteneur qui bien qu'associé à une loi CP, est également
|
|
// en réalité un conteneur CP2 car c'est comme cela qu'il a été créé via
|
|
// LoiCritere::LoiCritere::New_et_Initialise()
|
|
//*** là il va falloir revoir cette partie, car en fait ce n'est pas sûr qui faut utiliser un conteneur intermédiaire
|
|
//*** cela pose pas mal de pb, car il faut de nouveau mettre à jour les infos du conteneur initial, dans le cas des lois
|
|
//*** incrémentales, or pour l'instant ce n'est pas fait, ça c'est une erreur
|
|
/* *** l'idée initiale était que l'utilisation du 1D était transitoire en particulier d'une itération à l'autre.
|
|
donc supposons qu'à une itération on ait un plis dans une direction, on change toutes les grandeurs tensorielles en fonction
|
|
de cette direction (par exemple les tenseurs de ref) et donc à l'itération suivante on n'aura plus le même état de démarrage.
|
|
Mais ... en fait on remet les choses dans le repère initial à la fin du traitement:
|
|
save_resul.save_result_1DCP2->ChBase_des_grandeurs(gamma3D);
|
|
donc dans ce cas là, je ne comprends pas pourquoi on utilise pas directement le conteneur initial
|
|
Pour l'instant, je laisse en l'état (25 nov 2015) et après la réunion avec le CNES, je vais faire les tests en incrémental et modifier
|
|
car c'est possible qu'il y a quelque chose qui m'échappe.
|
|
Par contre à minima je remets les def de largeur, dans le conteneur initial
|
|
*/
|
|
}
|
|
|
|
|
|
LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble * save_result_1DCP2_inter // pour simplifier
|
|
= (LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble *)
|
|
(*(save_resul.liste_des_SaveResul.begin()));
|
|
// on affecte les grandeurs à un conteneur de travail
|
|
if (save_resul.save_result_1DCP2 == NULL)
|
|
{ SaveResul * inter = save_result_1DCP2_inter->Nevez_SaveResul();
|
|
save_resul.save_result_1DCP2 = ((LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble*) inter);
|
|
}
|
|
else // sinon on peut affecter directement
|
|
*(save_resul.save_result_1DCP2) = *((Loi_comp_abstraite::SaveResul *) save_result_1DCP2_inter);
|
|
|
|
// . -- changement de repère pour venir dans le nouveau repère ortho
|
|
// on change toutes les infos tensorielles contenues pour qu'elles soient exprimées dans le nouveau repère
|
|
// ceci ne concerne pas les grandeurs scalaires et les grandeurs exprimées dans le repère global
|
|
|
|
// là il faut utiliser une matrice 3x3 car on a systématiquement des grandeurs 3D, car la loi interne est 3D
|
|
save_resul.save_result_1DCP2->ChBase_des_grandeurs(beta3D,gamma3D);
|
|
|
|
// il faut initialiser les métriques en 1D à null, pour qu'elles soient recalculer systématiquement car par exemple
|
|
// la métrique à 0 change d'une itération à l'autre et d'un incrément à l'autre
|
|
{Deformation::Stmet met_ref_00,met_ref_t,met_ref_tdt; // def de métrique vide, qui vont être construite pendant le calcul
|
|
save_resul.save_result_1DCP2->Affectation_metriques_locale(met_ref_00,met_ref_t,met_ref_tdt);
|
|
};
|
|
|
|
loi_1DCP2_de_3D->IndiquePtIntegMecaInterne(ptintmeca_en_cours); // on transmet le ptintmeca global
|
|
loi_1DCP2_de_3D->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours
|
|
loi_1DCP2_de_3D->IndiqueSaveResult(save_resul.save_result_1DCP2); // on indique à la loi le conteneur qu'il faut utiliser
|
|
////-- debug --
|
|
{//#ifdef MISE_AU_POINT
|
|
//if ( save_resul.save_result_1DCP2->meti_ref_00.giB_ != NULL)
|
|
// cout << "\n debug Critere_plis_membrane ";
|
|
//#endif
|
|
////-- fin debug --
|
|
// -- appel du programme de contraintes doublement planes -> récup des résultats dans le nouveau repère
|
|
// 7 nov 2018: ?? normalement on est en base orthonormée ! bool base_orthonormee = false;
|
|
}
|
|
bool base_orthonormee = true;
|
|
// conditionnnement de la nouvelle métrique
|
|
Creation_metrique_a_partir_vecteurs_propres_pour_Umat1D(ex,beta);
|
|
// chgt de variable pour les autres variables de passage: contrainte, def, D, delta eps ...
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{
|
|
//cout << "\n debug LoiCritere::Critere_plis_membrane(";
|
|
// for (int a=1; a<4;a++)
|
|
// {cout << "\n giH_0("<<a<<"): "; (*(umat_cont_1D->giB_0))(a).Affiche();
|
|
// cout << "\n giH_0("<<a<<"): "; (*(umat_cont_1D->giB_0))(a).Affiche();
|
|
// };
|
|
cout << "\n epsBB initial en local 2D "; epsBB_tdt.Ecriture(cout);
|
|
// en orthonormee
|
|
Tenseur3BB epsBB_3_inter;
|
|
epsBB_tdt.BaseAbsolue(epsBB_3_inter,*(ex.giH_tdt));
|
|
cout << "\n epsBB d'entrée =en orthonormee "<<epsBB_3_inter;
|
|
cout << endl;
|
|
}
|
|
#endif
|
|
|
|
bool deux_plis = false;
|
|
Passage_3ou2_vers_1(gamma,sigHH_t,beta,DepsBB,epsBB_tdt
|
|
,delta_epsBB,deux_plis,save_resul.eps_pli);
|
|
|
|
// avant le calcul on initialise l'énergie
|
|
energ.Inita(0.); // car sinon on cumule l'énergie sans plis et avec plis !!
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
cout << "\n LoiCritere::Critere_plis_membrane ==> calcul en 1D CP2";
|
|
#endif
|
|
loi_1DCP2_de_3D->Calcul_dsigma_deps(base_orthonormee,sig_HH_t_1D,Deps_BB_1D
|
|
,eps_BB_1D,delta_eps_BB_1D,jacobien_0_1D,jacobien_tdt_1D
|
|
,sig_HH_1D,d_sigma_deps_1D,energ,energ_t,module_compressibilite, module_cisaillement
|
|
,*umat_cont_1D);
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 4)
|
|
{ if (Permet_affichage() > 5)
|
|
{cout << "\n\n ... LoiCritere::Critere_plis_membrane apres CP2: ";
|
|
cout << "\n sig_HH_1D "<< sig_HH_1D;
|
|
cout << "\n d_sigma_deps_1D "<< d_sigma_deps_1D;
|
|
};
|
|
if ((!isfinite(d_sigma_deps_1D(1,1,1,1))) || (isnan(d_sigma_deps_1D(1,1,1,1))) )
|
|
{cout << "\n *** attention d_sigma_deps_3D(1,1,1,1)= "<< d_sigma_deps_1D(1,1,1,1);
|
|
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
|
|
cout << "\n LoiCritere::Critere_plis_membrane ==> calcul en 1D CP2";
|
|
cout << flush;
|
|
//// --- pour le débug
|
|
// loi_1DCP2_de_3D->Calcul_dsigma_deps(base_orthonormee,sig_HH_t_1D,Deps_BB_1D
|
|
// ,eps_BB_1D,delta_eps_BB_1D,jacobien_0_1D,jacobien_tdt_1D
|
|
// ,sig_HH_1D,d_sigma_deps_1D,energ,energ_t,module_compressibilite, module_cisaillement
|
|
// ,*umat_cont_1D);
|
|
//// ---fin debug
|
|
};
|
|
Tenseur1HH deltasigHH = (d_sigma_deps_1D && delta_eps_BB_1D);
|
|
cout << "\n deltasigHH "<< deltasigHH;
|
|
cout << endl;
|
|
};
|
|
#endif
|
|
if (Permet_affichage() > 4)
|
|
if (//(sig_HH_1D(1,1) <= -1.)||
|
|
(std::isinf(sig_HH_1D(1,1)))||( std::isnan(sig_HH_1D(1,1))))
|
|
{cout << "\n *** attention sig_HH_1D(1,1) negatif ou infinie (debug LoiCritere::Critere_plis_membrane(.. \n";
|
|
cout << sig_HH_1D(1,1)
|
|
<< endl;
|
|
cout << "\n ***** erreur detectee : on continue quand meme ... ";
|
|
// Sortie(2);
|
|
};
|
|
|
|
// -- passage des résultats: sig, dsig, module, énergies dans l'ancien repère
|
|
// passage inverse: chgt de repère et de dimension pour les variables de passage
|
|
// et stockage des infos de doublement plane pour le post-traitement pour le prochain appel
|
|
Passage_1_vers_3ou2(gamma,sigHH,d_sigma_deps_1D
|
|
,beta,d_sigma_deps_inter,ex,V_Pr_H
|
|
,*(save_resul.V_P_sig));
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{ cout << "\n apres Passage_1_vers_3ou2, suite LoiCritere::Critere_plis_membrane";
|
|
cout << "\n sigHH en giB: "; sigHH.Ecriture(cout);
|
|
cout << "\n d_sigma_deps_inter en giB: ";d_sigma_deps_inter.Ecriture(cout);
|
|
Tenseur2HH deltasigHH = (d_sigma_deps_inter && delta_epsBB);
|
|
cout << "\n deltasig2HH en giB: "<< deltasigHH;
|
|
cout << endl;
|
|
};
|
|
#endif
|
|
|
|
if (Permet_affichage() > 4)
|
|
for (int al = 1; al<3;al++)
|
|
for (int be=1;be<3;be++)
|
|
if (//(sigHH(al,be) < 0.)||
|
|
(std::isinf(sigHH(al,be))||( std::isnan(sigHH(al,be)))))
|
|
{cout << "\n *** attention (sigHH("<<al<<","<<be
|
|
<< ")="<< sigHH(al,be) << " nan debug LoiCritere::Critere_plis_membrane( \n";
|
|
cout << "\n ***** erreur detectee : on continue quand meme ... ";
|
|
};
|
|
// -- mise à jour du save_result de contraintes planes avec celui des doublements planes en tenant
|
|
// compte du chgt de repère
|
|
//*** ça ce sera à revoir: cf. la remarque plus haut ***
|
|
save_resul.save_result_1DCP2->ChBase_des_grandeurs(beta3D,gamma3D);
|
|
// on met à jour les def de largeur et d'épaisseur sur le conteneur initial
|
|
save_result_1DCP2_inter->Transfert_var_hetb(*save_resul.save_result_1DCP2);
|
|
|
|
};
|
|
|
|
// deuxième type de calcul dans le cas d'un pli dans une seule direction
|
|
void LoiCritere::Deuxieme_type_calcul_en_un_pli
|
|
(const TenseurBB & epsBB_tdt,const TenseurBB & delta_epsBB
|
|
,const TenseurHH & sigHH_t
|
|
,double& jacobien_0,double& jacobien
|
|
,EnergieMeca & energ,const EnergieMeca & energ_t
|
|
,const TenseurBB& DepsBB
|
|
,double& module_compressibilite,double& module_cisaillement
|
|
,const Tableau <Coordonnee2H>& V_Pr_H
|
|
,TenseurHH& sigHH,TenseurHHHH& d_sigma_deps_inter
|
|
,const Coordonnee2& valPropreSig
|
|
,const Met_abstraite::Umat_cont& ex
|
|
)
|
|
{
|
|
SaveResul_LoiCritere & save_resul = *((SaveResul_LoiCritere*) saveResul);
|
|
// le conteneur Met_abstraite::Umat_cont& ex correspond à une métrique 2D
|
|
// il nous faut l'équivalent en 3D, on se sert de la loi cp et de la méthode ad hoc
|
|
const Met_abstraite::Umat_cont* ex3D = loi_2DCP_de_3D->Passage_metrique_ordre2_vers_3(ex);
|
|
// en fait le repère 2D a été complété par la normale aux 2 premiers vecteurs
|
|
|
|
// -- en entrée on connait la direction des vecteurs propres de contrainte
|
|
// correspondant au calcul sans prise en compte de pli
|
|
// on va calculer les bases associées: ViB et ViH: ici 3 vecteurs en 3D
|
|
// a) ViH = les vecteurs en contravariant: c'est en fait = V_Pr_H pour
|
|
// les 2 premiers vecteurs
|
|
// mais on doit passer de 2D en 3D, manuellement
|
|
ViH.CoordoH(1)(1) = V_Pr_H(1)(1);ViH.CoordoH(1)(2) = V_Pr_H(1)(2);
|
|
ViH.CoordoH(1)(3) = 0.; // normalement pas de composante suivant la normale
|
|
ViH.CoordoH(2)(1) = V_Pr_H(2)(1);ViH.CoordoH(2)(2) = V_Pr_H(2)(2);
|
|
ViH.CoordoH(2)(3) = 0.; // normalement pas de composante suivant la normale
|
|
// pour le 3 ième vecteur on se sert des grandeurs sauvegardées
|
|
for (int i=1;i<4;i++) // < 4 donc de 1 à 3
|
|
ViH.CoordoH(3)(i) = (*save_resul.V_P_sig)(3) * ex3D->giH_tdt->Coordo(i);
|
|
// b) ViB -> là il faut les calculer
|
|
// on choisit d'utiliser les grandeurs stockées qui représentent les
|
|
// les vecteurs orthonormées
|
|
for (int j=1;j<4;j++)
|
|
for (int i=1;i<4;i++) // < 4 donc de 1 à 3
|
|
ViB.CoordoB(j)(i) = (*save_resul.V_P_sig)(j) * ex3D->giB_tdt->Coordo(i);
|
|
|
|
// on calcul les def de plis : save_resul.eps_pli
|
|
Mat_pleine beta(2,2);
|
|
for (int i=1;i<3;i++) // < 3 donc de 1 à 2
|
|
for (int j=1;j<3;j++)
|
|
beta(i,j) = V_Pr_H(i)(j);
|
|
Mat_pleine gamma(beta.Inverse().Transpose()); // on calcule la matrice inverse
|
|
Tenseur2BB epsBB_tdt_2_inter(epsBB_tdt);epsBB_tdt_2_inter.ChBase(beta);
|
|
save_resul.eps_pli(1) = epsBB_tdt_2_inter(2,2); // on sauvegarde les def de plis
|
|
save_resul.eps_pli(2) = 0.; // cas d'un seul plis
|
|
|
|
|
|
#ifdef MISE_AU_POINT
|
|
// pour tester et afficher: calcul des matrices de changement de base
|
|
if (Permet_affichage() > 5)
|
|
{// les coordonnées des vecteurs propres sont exprimées dans l'ancienne base
|
|
// on a donc
|
|
// Vi_B(i) = gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
|
|
|
|
// Mat_pleine beta(2,2);
|
|
// for (int i=1;i<3;i++) // < 3 donc de 1 à 2
|
|
// for (int j=1;j<3;j++)
|
|
// beta(i,j) = V_Pr_H(i)(j);
|
|
Mat_pleine gamma(beta.Inverse().Transpose()); // on calcule la matrice inverse
|
|
// la matrice gamma correspond à la matrice de changement de base des g^i au gpH
|
|
// gamma est l'inverse transposée
|
|
// gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
|
|
// c-a-d= gp^i = gamma^i_j * g^j
|
|
|
|
// on va également définir une matrice de changement pour les tenseurs 3D
|
|
Tableau <Coordonnee3H> V_Pr3D_H(3);
|
|
// comme V_Pr est de dim 2, il faut transférer les coordonnées
|
|
V_Pr3D_H(1)(1) = V_Pr_H(1)(1);V_Pr3D_H(1)(2) = V_Pr_H(1)(2);
|
|
V_Pr3D_H(2)(1) = V_Pr_H(2)(1);V_Pr3D_H(2)(2) = V_Pr_H(2)(2);
|
|
// le troisième est en local 0,0,1
|
|
V_Pr3D_H(3)(3) = 1.;
|
|
Mat_pleine beta3D(3,3);// initialisée à 0
|
|
// d'abord les 2 vecteurs propres
|
|
for (int i=1;i<3;i++)
|
|
for (int j=1;j<4;j++)
|
|
beta3D(i,j) = V_Pr3D_H(i)(j);
|
|
// puis la normale
|
|
beta3D(3,3) = 1.;
|
|
Mat_pleine gamma3D(beta3D.Inverse().Transpose());
|
|
Tenseur2HH sigHH_2_inter(sigHH);
|
|
cout << "\n utilisation du repere ortho des vecteurs propres: V_P ";
|
|
cout << "\n sigHH= en giB : "; sigHH.Ecriture(cout);
|
|
cout << "\n base giB(1) "<< (* ex3D->giB_tdt)(1)
|
|
<< "\n base giB(2) "<< (* ex3D->giB_tdt)(2);
|
|
// en base orthonormee
|
|
Tenseur3HH sigHH_3_inter;
|
|
sigHH.BaseAbsolue(sigHH_3_inter,(* ex3D->giB_tdt));
|
|
cout << "\n sigHH= en absolue 3D : "<<sigHH_3_inter;
|
|
cout << "\n chgt de base: Vi_B(i) = gpB(i) = beta(i,j) * gB(j) --> beta= ";beta.Affiche();
|
|
cout << "\n sigHH_2_inter avant chgt= "<< sigHH_2_inter;
|
|
sigHH_2_inter.ChBase(gamma,false);
|
|
cout << "\n sigHH= en V_P : "<< sigHH_2_inter;
|
|
// on revient à la base de départ pour vérifier que l'opération inverse ne pose
|
|
// pas de pb
|
|
sigHH_2_inter.ChBase(beta.Transpose(),false);
|
|
cout << "\n verif: sigHH après retour au repère initial giB: ";
|
|
sigHH_2_inter.Ecriture(cout);
|
|
//cout << "\n gamma * beta = "; (gamma * beta).Affiche();
|
|
cout << "\n vecteurs propres en 3D ortho: "; //<< (*save_resul.V_P_sig);
|
|
cout << "\n VP1 "<< (*save_resul.V_P_sig)(1);
|
|
cout << "\n VP2 "<< (*save_resul.V_P_sig)(2);
|
|
cout << "\n VP3 "<< (*save_resul.V_P_sig)(3);
|
|
// vérif des vecteurs propres
|
|
Coordonnee3B VpB1; VpB1.Change_val((*save_resul.V_P_sig)(1));
|
|
Coordonnee3H VpropreH1 = sigHH_3_inter * VpB1;
|
|
//cout << "\n sigHH_3_inter * VP1 "<< VpropreH1;
|
|
VpropreH1 /= (valPropreSig(1)+ConstMath::petit);
|
|
cout << "\n (sigHH_3_ortho * VP1) / (lambda1+epsilon)= "<< VpropreH1;
|
|
Coordonnee3B VpB2; VpB2.Change_val((*save_resul.V_P_sig)(2));
|
|
Coordonnee3H VpropreH2 = sigHH_3_inter * VpB2;
|
|
//cout << "\n sigHH_3_inter * VP2 "<< VpropreH2;
|
|
VpropreH2 /= (valPropreSig(2)+ConstMath::petit);
|
|
cout << "\n (sigHH_3_ortho * VP2) / (lambda2+epsilon)= "<< VpropreH2;
|
|
};
|
|
#endif
|
|
|
|
LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble * save_result_1DCP2_inter // pour simplifier
|
|
= (LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble *)
|
|
(*(save_resul.liste_des_SaveResul.begin()));
|
|
|
|
loi_1DCP2_de_3D->IndiqueSaveResult(save_result_1DCP2_inter); // on indique à la loi le conteneur qu'il faut utiliser
|
|
loi_1DCP2_de_3D->IndiquePtIntegMecaInterne(ptintmeca_en_cours); // on transmet le ptintmeca global
|
|
loi_1DCP2_de_3D->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours
|
|
bool base_orthonormee = false; // ici on va travailler avec les gi
|
|
// on utilise donc la métrique passée en paramètre
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{
|
|
//cout << "\n debug LoiCritere::Critere_plis_membrane(";
|
|
// for (int a=1; a<4;a++)
|
|
// {cout << "\n giH_0("<<a<<"): "; (*(umat_cont_1D->giB_0))(a).Affiche();
|
|
// cout << "\n giH_0("<<a<<"): "; (*(umat_cont_1D->giB_0))(a).Affiche();
|
|
// };
|
|
cout << "\n epsBB initial en local 2D "; epsBB_tdt.Ecriture(cout);
|
|
// en orthonormee
|
|
Tenseur3BB epsBB_3_inter;
|
|
epsBB_tdt.BaseAbsolue(epsBB_3_inter,*(ex.giH_tdt));
|
|
cout << "\n epsBB d'entrée =en orthonormee "<<epsBB_3_inter;
|
|
cout << endl;
|
|
}
|
|
#endif
|
|
|
|
// passage des informations liées à la déformation de 2 vers 3
|
|
Passage_deformation_contrainte_ordre2_vers_3(DepsBB,epsBB_tdt,delta_epsBB,sigHH_t);
|
|
// avant le calcul on initialise l'énergie
|
|
energ.Inita(0.); // car sinon on cumule l'énergie sans plis et avec plis !!
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{
|
|
cout << "\n LoiCritere::Critere_plis_membrane ==> calcul en 3D CP2";
|
|
}
|
|
#endif
|
|
loi_1DCP2_de_3D->Calcul_dsigma_deps
|
|
(base_orthonormee,&ViB,&ViH,sig_HH_t_3D,Deps_BB_3D
|
|
,eps_BB_3D,delta_eps_BB_3D,jacobien_0,jacobien
|
|
,sig_HH_3D,d_sig_deps_3D_HHHH
|
|
,energ,energ_t
|
|
,module_compressibilite,module_cisaillement
|
|
,*ex3D);
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{ cout << "\n\n ... LoiCritere::Deuxieme_type_calcul_en_un_pli apres CP2: ";
|
|
cout << "\n sigHH_3D en giB: "; sig_HH_3D.Ecriture(cout);
|
|
if (Permet_affichage() > 6)
|
|
{cout << "\n d_sig_deps_3D_HHHH en giB: ";d_sig_deps_3D_HHHH.Ecriture(cout);};
|
|
Tenseur3HH deltasigHH = (d_sig_deps_3D_HHHH && delta_eps_BB_3D);
|
|
cout << "\n deltasig3HH en giB: "<< deltasigHH;
|
|
cout << endl;
|
|
};
|
|
#endif
|
|
|
|
if (Permet_affichage() > 4)
|
|
{bool erreur = false;
|
|
for (int al = 1; al<4;al++)
|
|
for (int be=1;be<4;be++)
|
|
if (//(sigHH(al,be) < 0.)||
|
|
(std::isinf(sig_HH_3D(al,be))||( std::isnan(sig_HH_3D(al,be)))))
|
|
{cout << "\n *** attention (sigHH_3D("<<al<<","<<be
|
|
<< ")="<< sig_HH_3D(al,be) << " valeur nan LoiCritere::Critere_plis_membrane( \n"
|
|
<< "\n LoiCritere::Deuxieme_type_calcul_en_un_pli apres CP2:.. \n";
|
|
erreur = true;
|
|
};
|
|
if (erreur)
|
|
{cout << "\n sigHH= ";sig_HH_3D.Ecriture(cout);
|
|
cout << "\n **** erreur nan sur loi critere !! **** on continue quand meme ... ";
|
|
// Sortie(2);
|
|
};
|
|
};
|
|
|
|
// passage des informations liées à la nouvelle contrainte de 2 vers 3
|
|
// et à l'opérateur tangent : méthode 2
|
|
// non !! pas de mise à jour --> // mise à jour de save_resul.eps_pli
|
|
// *** a supprimer le passage de save_resul.eps_pli
|
|
// car c'est une erreur
|
|
|
|
|
|
bool deux_plis = false;
|
|
Passage_contrainte_et_operateur_tangent_ordre2_vers_3
|
|
(sigHH,d_sigma_deps_inter,deux_plis,save_resul.eps_pli);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
{ cout << "\n retour en dim 2 : ";
|
|
cout << "\n sigHH en giB: "; sigHH.Ecriture(cout);
|
|
cout << "\n d_sig_deps_HHHH en giB: ";d_sigma_deps_inter.Ecriture(cout);
|
|
Tenseur2HH deltasigHH = (d_sigma_deps_inter && delta_epsBB);
|
|
cout << "\n deltasig2HH en giB: "<< deltasigHH;
|
|
cout << endl;
|
|
};
|
|
#endif
|
|
|
|
};
|
|
|
|
|
|
// fonction pre_critère de plissement de membrane
|
|
// dans le cas d'un critère pli (plutôt seconde méthode), l'incrément de déformation
|
|
// dépend de la situation précédente: avec pli ou pas
|
|
// du coup on utilise un delta_eps intermédiaire au lieu du delta cinématique
|
|
// ce deltat_eps doit être précalculé
|
|
// en entrée:
|
|
// implicit : si oui on est en implicite, sinon en explicite
|
|
// retour:
|
|
void LoiCritere::Pre_Critere_plis_membrane
|
|
(TenseurBB & epsBB_tdt_,TenseurBB & delta_epsBB
|
|
,const Met_abstraite::Expli_t_tdt* ex_expli
|
|
,const Met_abstraite::Impli* ex_impli
|
|
,const Met_abstraite::Umat_cont* ex_umat
|
|
)
|
|
{ // on commence par récupérer les métriques en fonctions du conteneur ex valide
|
|
const Tenseur2HH * gijHH_; // pour simplifier
|
|
const Tenseur2BB * gijBB_; // " " " "
|
|
const Tenseur2HH * gijHH_t_; // " " " "
|
|
const Tenseur2BB * gijBB_t_; // " " " "
|
|
if (ex_expli != NULL)
|
|
{ gijHH_ = ((Tenseur2HH*) ex_expli->gijHH_tdt); // pour simplifier
|
|
gijBB_ = ((Tenseur2BB*) ex_expli->gijBB_tdt); // " " " "
|
|
gijHH_t_ = ((Tenseur2HH*) ex_expli->gijHH_t); // " " " "
|
|
gijBB_t_ = ((Tenseur2BB*) ex_expli->gijBB_t); // " " " "
|
|
}
|
|
else if (ex_impli != NULL)
|
|
{ gijHH_ = ((Tenseur2HH*) ex_impli->gijHH_tdt); // pour simplifier
|
|
gijBB_ = ((Tenseur2BB*) ex_impli->gijBB_tdt); // " " " "
|
|
gijHH_t_ = ((Tenseur2HH*) ex_impli->gijHH_t); // " " " "
|
|
gijBB_t_ = ((Tenseur2BB*) ex_impli->gijBB_t); // " " " "
|
|
}
|
|
else if (ex_umat != NULL)
|
|
{ gijHH_ = ((Tenseur2HH*) ex_umat->gijHH_tdt); // pour simplifier
|
|
gijBB_ = ((Tenseur2BB*) ex_umat->gijBB_tdt); // " " " "
|
|
gijHH_t_ = ((Tenseur2HH*) ex_umat->gijHH_t); // " " " "
|
|
gijBB_t_ = ((Tenseur2BB*) ex_umat->gijBB_t); // " " " "
|
|
}
|
|
else
|
|
{cout << "\n **** erreur: pas de conteneur de metrique disponible ..."
|
|
<< "\n LoiCritere::Pre_Critere_plis_membrane(...";
|
|
Sortie(1);
|
|
};
|
|
|
|
const Tenseur2HH & gijHH = *gijHH_; // pour simplifier
|
|
const Tenseur2BB & gijBB = *gijBB_; // " " " "
|
|
const Tenseur2HH & gijHH_t = *gijHH_t_; // " " " "
|
|
const Tenseur2BB & gijBB_t = *gijBB_t_; // " " " "
|
|
|
|
// et la déformation
|
|
const Tenseur2BB & epsBB_tdt = *((Tenseur2BB*) &epsBB_tdt_); // " " " "
|
|
|
|
SaveResul_LoiCritere & save_resul = *((SaveResul_LoiCritere*) saveResul);
|
|
int retour = 0; // init par défaut
|
|
LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble * save_result_1DCP2_inter // pour simplifier
|
|
= (LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble *)
|
|
(*(save_resul.liste_des_SaveResul.begin()));
|
|
|
|
// l'objectif est de mettre à jour : delta_eps_BB_2D
|
|
delta_eps_BB_2D = delta_epsBB; // initialisation par défaut
|
|
|
|
// on regarde l'état précédent: sans plis, avec un pli, avec 2 plis ?
|
|
int cas_cal_plis_t = save_resul.cas_cal_plis_t;
|
|
// en fait seule les cas avec un deux plis, sans erreur, sont exploitables
|
|
// ou le cas sans pli, pour tous les autres cas, on considèrera qu'il n'y avait pas de plis
|
|
// car les infos sauvegardées (def méca) ne sont pas correctes
|
|
switch (save_resul.cas_cal_plis_t)
|
|
{ case -1: case 1: case -2: case -3: case -4:// pas de plis
|
|
{// a priori rien à faire, delta_eps_BB_2D a déjà été initialisée
|
|
break;
|
|
}
|
|
case 2: case 3: // plis dans une seule directions ou dans deux directions
|
|
{ // L'incrément de déformation, entre t et t+dt est obtenu par la différence entre les
|
|
// déformations t+dt et les déformations mécaniques à t, transportée de manière cohérente
|
|
// avec la dérivée de Jauman à t+dt: 1/2 deux fois covariant + deux fois contra
|
|
// a) on commence par récupérer les def méca à l'instant t
|
|
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_result_1DCP2_inter->eps_mecaBB_t));
|
|
// la déformation 33 est directement gérée par la loi 2D CP via les variations d'épaisseurs
|
|
// sachant que l'axe 3 est indépendant des deux autres axes donc on ne s'en occupe pas
|
|
eps_BB_2D_t.Affectation_3D_a_2D(eps_mecaBB_t);
|
|
// b) on fait le transport cohérent avec Jauman
|
|
eps_HH_2D_t = gijHH_t * eps_BB_2D_t * gijHH_t;
|
|
delta_epsBB = epsBB_tdt - 0.5 * (gijBB * eps_HH_2D_t * gijBB + eps_BB_2D_t);
|
|
break;
|
|
}
|
|
case 0:
|
|
// cas où aucun calcul n'a encore été effectué on n'a rien n'a faire
|
|
break;
|
|
default: // normalement n'arrive jamais car Calcul_directions_plis a toujours une solution diff de 0
|
|
{ cout << "\n *** cas d'erreur inconnue dans le calcul de la direction des plis "
|
|
<< " save_resul.cas_cal_plis= "<< save_resul.cas_cal_plis
|
|
<< "\n LoiCritere::Pre_Critere_plis_membrane(... ";
|
|
Sortie(1);
|
|
break;
|
|
}
|
|
};
|
|
};
|
|
|
|
// passage des informations liées à la déformation de 2 vers 3
|
|
void LoiCritere::Passage_deformation_contrainte_ordre2_vers_3
|
|
(const TenseurBB& DepsBB,const TenseurBB & epsBB_tdt
|
|
,const TenseurBB & delta_epsBB,const TenseurHH& sig_HH_t)
|
|
|
|
{ // on complète avec des 0
|
|
bool plusZero = true;
|
|
Deps_BB_3D.Affectation_2D_a_3D(*((Tenseur2BB*) &DepsBB),plusZero);
|
|
eps_BB_3D.Affectation_2D_a_3D(*((Tenseur2BB*) &epsBB_tdt),plusZero);
|
|
delta_eps_BB_3D.Affectation_2D_a_3D(*((Tenseur2BB*) &delta_epsBB),plusZero);
|
|
sig_HH_t_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sig_HH_t),plusZero);
|
|
};
|
|
|
|
// passage des informations liées à la nouvelle contrainte de 2 vers 3
|
|
// et à l'opérateur tangent : méthode 2
|
|
void LoiCritere::Passage_contrainte_et_operateur_tangent_ordre2_vers_3
|
|
(TenseurHH& sig_HH_tdt_,TenseurHHHH& d_sig_deps_HHHH_
|
|
,const bool& deux_plis,Coordonnee2& eps_pli)
|
|
{ // normalement les directions 1 et 2 sont identiques en 2D et 3D
|
|
// on recopie donc simplement les grandeurs
|
|
Tenseur2HH & sigHH = *((Tenseur2HH*) &sig_HH_tdt_); // " " " "
|
|
sigHH.Affectation_3D_a_2D(sig_HH_3D);
|
|
// idem pour l'opérateur tangent
|
|
Tenseur2HHHH & d_sig_deps_HHHH = *((Tenseur2HHHH*) &d_sig_deps_HHHH_); // pour simplifier
|
|
d_sig_deps_HHHH.TransfertDunTenseurGeneral(d_sig_deps_3D_HHHH);
|
|
|
|
/* Lorsqu'il y a un pli dans une seule direction, le premier vecteur indique la direction selon laquelle la membrane est en traction. Son intensité est égale à la déformation mécanique dans cette même direction. Le second vecteur est normal à la première direction, son intensité est égale à la déformation cinématique selon cette direction, c'est-à-dire calculée à l'aide des déplacements des noeuds du maillage. Cette déformation est différente de la déformation mécanique calculée par la loi 3D, c'est-à-dire issue de la contraction de la zone, due à la traction suivant la première direction. Cette déformation mécanique peut-être récupéré via la deformation "DEF\_LARGEUR" associée à la loi de comportement.
|
|
|
|
Lorsqu'il y a des plis dans deux directions, les deux vecteurs indiquent les directions principales d'apparition de plis. L'intensité des vecteurs est égale à la déformation cinématique dans la direction associée.
|
|
|
|
Lorsqu'il n'y a pas de plis, l'intensité des vecteurs est nulle.
|
|
|
|
*/
|
|
|
|
//--- modif 24 sept 2020:
|
|
// les def pour les plis ont déjà été calculées dans LoiCritere::Passage_3ou2_vers_1(..
|
|
// et il s'agit uniquement de def cinématique donc la suite n'est pas nécessaire et même incorrecte
|
|
|
|
// SaveResul_LoiCritere & save_resul = *((SaveResul_LoiCritere*) saveResul);
|
|
// LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble * save_result_1DCP2_inter // pour simplifier
|
|
// = (LoiContraintesPlanesDouble::SaveResul_LoiContraintesPlanesDouble *)
|
|
// (*(save_resul.liste_des_SaveResul.begin()));
|
|
// Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_result_1DCP2_inter->eps_mecaBB));
|
|
//
|
|
// eps_pli(1) = eps_mecaBB(2,2); // on sauvegarde les def de plis
|
|
// if (deux_plis) // cas où il y aurait des plis dans les deux sens
|
|
// eps_pli(2) = eps_mecaBB(1,1);
|
|
// else eps_pli(2) = 0.; // cas d'un seul plis
|
|
// Tenseur3BB& eps_P_mecaBB = *((Tenseur3BB*) (save_result_1DCP2_inter->eps_P_mecaBB));
|
|
//
|
|
// eps_pli(1) = eps_P_mecaBB(2,2); // on sauvegarde les def de plis
|
|
// if (deux_plis) // cas où il y aurait des plis dans les deux sens
|
|
// eps_pli(2) = eps_P_mecaBB(1,1);
|
|
// else eps_pli(2) = 0.; // cas d'un seul plis
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|