// 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 "Algori.h"
#include "MatBand.h"
#include "Assemblage.h"
#include "ElemMeca.h"
// ------ fait partie de la classe Algori -----------
// initialisation du calcul de la remontée des contraintes aux noeuds
// initialisation valable également dans le cas où aucun calcul n'a précédé ce calcul
void Algori::InitRemontSigma(LesMaillages * lesMail,LesReferences* lesRef,
DiversStockage* ,Charge* ,
LesCondLim* lesCondLim,LesContacts* ,Resultats* )
{//mise en place des coordonnées à t, qui sont nécessaire pour calculer les éléments
// géométriques sur une pièce déformée, même si les ddl sont autres que les déplacements
// dans le cas où les coordonnées existent déjà, il n'y a aucune modification
lesMail->Insert_coord1();
// l'idée ici est d'initialiser la gestion d'un nouveau groupe de ddl
// ==== 1 - on commence par inactiver les ddl qui sont en cours
lesMail->Inactive_ddl();
// ==== 2 - on définit un nouveau groupe de ddl avec la gestion rapide associé
// introduction des ddl de contraintes (s'ils n'existent pas déjà sinon aucune action)
lesMail->Plus_ddl_Sigma();
////--- debug
//cout << "\n debug Algori::InitRemontSigma ";
//lesMail->Noeud_LesMaille(1,1).Affiche(1);
//cout << endl;
////--- fin debug
// mise a zero des ddl actifs et on laisse les ddl à t+dt
// mais il faut d'abord activer les ddl sigma
lesMail->Active_ddl_Sigma();
// le deuxième paramètre de ZeroDdl est false pour indiquer que les coordonnées à t
// ne sont pas à initialiser aux valeurs à t = 0, ceci permet par exemple de travailler
// sur une structure déformée si elle a été calculée auparavant
lesMail->ZeroDdl(true,false); // mise à jour de tous les ddl actifs sauf de déplacement
// ==== 3 - def de l'assemblage pour tous les ddl de type sigma
nbddl_sigeps = lesMail->NbTotalDdlActifs(); // nb total de ddl de sigma
// on définit un nouveau cas de second membre
// à travers la définition d'une instance de la classe assemblage
// si elle existe déjà on l'utilise directement
if (cas_ass_sig == NULL)
// cas où les conteneurs sont à créer
{ cas_ass_sig = new Nb_assemb(lesMail->InitNouveauCasAssemblage(1));
Ass_sig = new Assemblage(*cas_ass_sig);
}
// mise a jour des pointeurs d'assemblage
lesMail->MiseAJourPointeurAssemblage(Ass_sig->Nb_cas_assemb());
// ==== 4 - def de l'assemblage pour tous les ddl de type sigma11 uniquement (c-a-d le premier ddl)
// on rend actif le premier ddl mécanique, en fait on se servira des
//autres qui sont supposé être à la suite mais pour éviter d'avoir une
//grosse matrice de raideur le second membre est décomposé en plusieurs
//SM, chacun relatif à "un" ddl de contrainte
lesMail->Inactive_ddl_Sigma();
lesMail->Active_premier_ddl_Sigma();
// récup du nombre de composantes du tenseur de contrainte en absolu
int nbcomposante = ParaGlob::NbCompTens ();
nbddl_sigeps11 = lesMail->NbTotalDdlActifs(SIG11) / nbcomposante; // nb total de ddl de sigma11
// on définit un nouveau cas de second membre
// à travers la définition d'une instance de la classe assemblage
// si elle existe déjà on l'utilise directement
if (cas_ass_sig11 == NULL)
// cas où les conteneurs sont à créer
{ cas_ass_sig11 = new Nb_assemb(lesMail->InitNouveauCasAssemblage(1));
Ass_sig11 = new Assemblage(*cas_ass_sig11);
}
// mise a jour des pointeurs d'assemblage
lesMail->MiseAJourPointeurAssemblage(Ass_sig11->Nb_cas_assemb());
// ==== 5 - définition de la matrice pour la résolution du problème de minimisation
// dans le cas où c'est déjà fait on ne change rien
if (mat_glob_sigeps == NULL)
//on est obligé de faire une petite manipe avec un tableau, car l'appel à changé
{ Tableau < Mat_abstraite*> tab_inter(1); tab_inter(1)=NULL;
Choix_matriciel(nbddl_sigeps11,tab_inter,lesMail,lesRef
,Ass_sig11->Nb_cas_assemb(),lesCondLim);
mat_glob_sigeps = tab_inter(1);
};
// choix de la résolution
mat_glob_sigeps->Change_Choix_resolution(pa.Type_resolution(),pa.Type_preconditionnement());
// def du tableau de seconds membres globaux avec initialisation à zéro
if (vglob_sigeps == NULL)
vglob_sigeps = new Tableau <Vecteur> (nbcomposante,Vecteur(nbddl_sigeps11));
// définition du vecteur global de tous les ddl de sigma
if (vglobal_sigeps == NULL)
vglobal_sigeps = new Vecteur(nbcomposante*nbddl_sigeps11);
};
// initialisation du calcul de la remontée des déformations aux noeuds
// initialisation valable également dans le cas où aucun calcul n'a précédé ce calcul
void Algori::InitRemontEps(LesMaillages * ,LesReferences* ,
DiversStockage* ,Charge* ,
LesCondLim* ,LesContacts* ,Resultats* )
{// pour l'instant rien
};
// initialisation du calcul d'erreur aux éléments
// initialisation valable également dans le cas où aucun calcul n'a précédé ce calcul
void Algori::InitErreur(LesMaillages * lesMail,LesReferences* ,
DiversStockage* ,Charge* ,
LesCondLim* ,LesContacts* ,Resultats* )
{// pour la remontée des valeurs aux noeuds
// l'idée ici est d'initialiser la gestion d'un nouveau groupe de ddl
// ==== 1 - on commence par inactiver les ddl qui sont en cours
lesMail->Inactive_ddl();
// ==== 2 - on définit un nouveau groupe de ddl avec la gestion rapide associé
// introduction des ddl d'erreur (s'ils n'existent pas déjà sinon aucune action)
lesMail->Plus_ddl_Erreur();
// activation des ddl Erreur
lesMail->Active_ddl_Erreur();
// mise a zero des ddl actifs et on laisse les ddl à t+dt
// le deuxième paramètre de ZeroDdl est false pour indiquer que les coordonnées à t
// ne sont pas à initialiser aux valeurs à t = 0, ceci permet par exemple de travailler
// sur une structure déformée si elle a été calculée auparavant
lesMail->ZeroDdl(true,false); // mise à jour de tous les ddl actifs sauf de déplacement
// ==== 3 - def de l'assemblage pour tous les ddl de type erreur
nbddl_err= lesMail->NbTotalDdlActifs(); // nb total de ddl d'erreur
// on définit un nouveau cas de second membre
// à travers la définition d'une instance de la classe assemblage
// si elle existe déjà on l'utilise directement
if (cas_ass_err == NULL)
// cas où les conteneurs sont à créer
{ cas_ass_err = new Nb_assemb(lesMail->InitNouveauCasAssemblage(1));
Ass_err = new Assemblage(*cas_ass_err);
}
// mise a jour des pointeurs d'assemblage
lesMail->MiseAJourPointeurAssemblage(Ass_err->Nb_cas_assemb());
// ==== 4 - définition de la matrice pour la résolution du problème de minimisation
// dans le cas où c'est déjà fait on ne change rien
// ici contrairement au cas des contraintes où des déformations, on considère qu'il y a déjà eu
// la remonté des contraintes, qui est obligatoire pour calculer l'erreur aux éléments qui elle est
// obligatoire pour calculer l'erreur aux noeuds => on utilise la même matrice de raideur
if (mat_glob_sigeps == NULL)
{ cout << "\n erreur la remontee des contraintes n'a pas ete faite, on ne peut pas continuer"
<< "\n Algori::InitErreur(...";
Sortie(1);
};
// def du tableau du second membre global avec initialisation à zéro
// on utilise la place de celui utilisé pour sigma
if (vglob_sigeps == NULL)
{ cout << "\n erreur la remontee des contraintes n'a pas ete faite, on ne peut pas continuer"
<< "\n Algori::InitErreur(...";
Sortie(1);
};
};
// calcul de la remontée des contraintes aux noeuds
void Algori::RemontSigma(LesMaillages * lesMail)
{ // affichage
if (ParaGlob::NiveauImpression() >= 2 )
cout << "\n================================================================="
<< "\n REMONTEE DES CONTRAINTES AUX NOEUDS "
<< "\n=================================================================\n";
lesMail->Inactive_ddl(); // inactive les ddl en cours
// mise a zero des ddl actifs et on laisse les ddl à t+dt
// mais il faut d'abord activer les ddl sigma
lesMail->Active_ddl_Sigma();
// le deuxième paramètre de ZeroDdl est false pour indiquer que les coordonnées à t
// ne sont pas à initialiser aux valeurs à t = 0, ceci permet par exemple de travailler
// sur une structure déformée si elle a été calculée auparavant
lesMail->ZeroDdl(true,false);
// on rend actif le premier ddl mécanique, en fait on se servira des
//autres qui sont supposé être à la suite mais pour éviter d'avoir une
//grosse matrice de raideur le second membre est décomposé en plusieurs
//SM, chacun relatif à "un" ddl de contrainte
lesMail->Inactive_ddl_Sigma();
lesMail->Active_premier_ddl_Sigma();
// def du nombre de composantes du tenseur de contrainte en absolu
int nbcomposante = ParaGlob::NbCompTens ();
// pour simplifier
Mat_abstraite& matglob = *mat_glob_sigeps;
mat_glob_sigeps->Initialise(0.);
Tableau <Vecteur>& vglob = *vglob_sigeps;
Vecteur& vglobal = *vglobal_sigeps;
vglobal_sigeps->Zero();
// ----- Calcul de la raideur et des seconds membres
Tableau <Vecteur> * sol; // pointeur du vecteur solution
// boucle sur les elements
int nbMaillage = lesMail->NbMaillage();
for (int nbMail =1; nbMail<= nbMaillage; nbMail++)
{int nbelement = lesMail->Nombre_element(nbMail);
for (int ne=1; ne<= nbelement;ne++)
{
// calcul du résidu et de la matrice de raideur pour le calcul d'erreur
ElemMeca & el = *((ElemMeca*) &lesMail->Element_LesMaille(nbMail,ne)); // l'element
Element::Er_ResRaid resu = el.ContrainteAuNoeud_ResRaid();
Tableau<Noeud *>& taN = el.Tab_noeud(); // tableau de noeuds de l'el
// création du tableau d'un ddl SIG11, pour chaque noeud concerné des
// élément traité
const DdlElement& TabDdl = el.Tableau_de_Sig1() ;
// assemblage de la raideur
Ass_sig11->AssembMatSym (matglob,*(resu.raidErr),TabDdl,taN); // de la raideur
// assemblage des seconds membres
Ass_sig11->AssemSM (vglob,*(resu.resErr),TabDdl,taN); // du second membre
// cout << "numéro =" << ne << endl;
// (resu.res)->Affiche();
// (resu.raid)->Affiche();
// assemblage
};
};
// affichage éventuelle de la matrice de raideur et du second membre
if (ParaGlob::NiveauImpression() >= 10)
{ string entete = " affichage de la matrice de raideur pour le calcul d'erreur ";
matglob.Affichage_ecran(entete);
entete = " affichage du second membre nb= ";
int nbsm = vglob.Taille();
for (int i=1;i<=nbsm;i++)
vglob(i).Affichage_ecran(entete + ChangeEntierSTring(i) +" ");
};
//------------- resolution -----------------
// matglob.Affiche1(1,6,1,1,6,1);
// for (int i=1;i<= vglob.Taille();i++)
// vglob(i).Affiche();
// vglobal.Affiche();
// cout << "\n entree une lettre"; int toto; cin >> toto;
// ici sol en fait = vecglob qui est ecrase par la resolution
// la résolution se fait en deux temps pour que la même matrice de raideur puisse servir plusieur fois
matglob.Preparation_resol(); // préparation
sol = &(matglob.Simple_Resol_systID(vglob)); // résolution sans modification de la matrice de raideur
// affichage éventuelle de la solution
if (ParaGlob::NiveauImpression() >= 10)
{ string entete = " affichage de la solutions pour le calcul d'erreur ";
entete = " affichage de la composante nb= ";
int nbsm = vglob.Taille();
for (int i=1;i<=nbsm;i++)
(*sol)(i).Affichage_ecran(entete + ChangeEntierSTring(i) +" ");
};
// sol = &(matglob.Resol_systID(vglob));
// sol->Affiche();
////---debug
// cout << "\n debug Algori::RemontErreur ";
// (*sol)(1).Affiche();
// (*sol)(2).Affiche();
// (*sol)(3).Affiche();
//// (*sol)(4).Affiche();
//// (*sol)(5).Affiche();
//// (*sol)(6).Affiche();
cout << endl;
////---fin debug
// ------------- résultat transféré dans les noeuds ---------
// construction du vecteur avec tous les ddl sigma
int iposition = 1;
for (int j=1; j<= nbddl_sigeps11; j++)
for (int i=1;i<=nbcomposante;i++,iposition++)
vglobal(iposition) = vglob(i)(j);
// activation de tous les ddl sigma
lesMail->Active_ddl_Sigma();
// modif 6 nov 2014
// // actualisation de tous les sigma a t
lesMail->PlusDelta_t(vglobal,Ass_sig->Nb_cas_assemb());
// actualisation de tous les sigma a tdt, et on passe les info
// car c'est les grandeurs de tdt qui sont visualisées !
lesMail->ChangeDdla_tdt(vglobal,Ass_sig->Nb_cas_assemb());
// fin des calculs
};
// calcul de la remontée des contraintes aux
// noeuds, puis de l'erreur, après un calcul de mécanique.
// calcul de la remontée des déformations aux noeuds
void Algori::RemontEps(LesMaillages * )
{// pour l'instant rien
};
void Algori::RemontErreur(LesMaillages * lesMail)
{
if (ParaGlob::NiveauImpression() >= 2 )
cout << "\n=========================================================================="
<< "\n CALCUL D'ERREUR AUX ELEMENTS ET REMONTEE DES VALEURS AUX NOEUDS "
<< "\n=========================================================================\n ";
// 1) Calcul de l'erreur sur chaque élément
// type indique le type d'erreur retenue
// type = 1 : cas d'un calcul aux moindres carrés
// et retour de l'erreur relative totale sur les maillages en cours
int type = 1;
// Calcul de l'erreur sur l'ensemble des éléments
// type indique le type d'erreur retenue
// type = 1 : cas d'un calcul aux moindres carrés
// et retour un tableau de tableau de grandeurs sur les maillages en cours
// ret(i) : concerne le maillage i
// ret(i)(1) : somme des erreurs sur l'ensemble des éléments: est homogêne à
// un |delta contrainte| * domaine
// ret(i)(2) : somme de la grandeur de ref du calcul d'erreur sur l'ensemble des
// éléments: est homogêne à une |contrainte| * domaine
// ret(i)(3) : le maxi pour les tous les éléments de |delta contrainte| * domaine
// ret(i)(4) : le maxi pour les tous les éléments de |contrainte| * domaine
// ret(i)(5) : le maxi pour les tous les éléments de l'erreur relative
Tableau <Tableau <double > > ret(lesMail->ErreurSurChaqueElement(type));
if (ParaGlob::NiveauImpression() >= 2 )
{int nbmaillage = ret.Taille();
for (int i=1;i<=nbmaillage;i++)
{cout << "\n --- maillage: "<<i
<< "\n cumul des erreurs sur les elements: " << ret(i)(1)
<< " versus cumul de la |grandeur| observee = "<< ret(i)(2)
<< "\n ==> erreur cumulee relative= "<<ret(i)(1)/MaX(ret(i)(3), ConstMath::trespetit)
<< "\n |maxi| sur les elements de l'erreur= "<< ret(i)(3)
<< " versus |maxi| sur les element de la |grandeur| observee = "<< ret(i)(4)
<< "\n ==> maxi de l'erreur relative sur les element = "<<ret(i)(5)
<< '\n';
};
};
// 2) remontée des valeurs aux noeuds
// on inactive les ddl
lesMail->Inactive_ddl();
// activation des ddl Erreur
lesMail->Active_ddl_Erreur();
// mise a zero des ddl actifs et on laisse les ddl à t+dt
lesMail->ZeroDdl(true,false);
// mise a jour des pointeurs d'assemblage
lesMail->MiseAJourPointeurAssemblage(Ass_err->Nb_cas_assemb());
// pour le second membre on utilise le premier vecteur de second membre utilisé
// pour les contraintes : vglob(1) et localement resu.resErr(1)
// on ne recalcule pas la raideur, mais on utilise celle calculée lors de la résolution
// pour les contraintes. Il y a sauvegarde de la matrice de raideur lors de la
// résolution, on peut donc la réutiliser.
// ----- Calcul du second membre
// on commence par initialiser le second membre à zero
Vecteur & vgloberr = (*vglob_sigeps)(1);
vgloberr.Zero();
Mat_abstraite& matglob = *mat_glob_sigeps;
// boucle sur les elements
int nbMaillage = lesMail->NbMaillage();
for (int nbMail =1; nbMail<= nbMaillage; nbMail++)
{int nbelement = lesMail->Nombre_element(nbMail);
for (int ne=1; ne<= nbelement;ne++)
{
// calcul du résidu
ElemMeca & el = *((ElemMeca*) &lesMail->Element_LesMaille(nbMail,ne)); // l'element
Element::Er_ResRaid resu = el.ErreurAuNoeud_ResRaid();
Tableau<Noeud *>& taN = el.Tab_noeud(); // tableau de noeuds de l'el
// assemblage du second membre
Ass_err->AssemSM (vgloberr,*((*(resu.resErr))(1)),el.Tableau_de_ERREUR(),taN);
// cout << "numéro =" << ne << endl;
// (resu.res)->Affiche();
// (resu.raid)->Affiche();
// assemblage
}
}
//------------- resolution -----------------
// matglob.Affiche(1,4,1,1,4,1);
// for (int i=1;i<= vglob.Taille();i++)
// vglob(i).Affiche();
// matglob.Affiche(1,8,1,7,8,1);
// matglob.Affiche(1,18,1,13,18,1);
// Vecteur * solEr = &(matglob.Resol_systID_sans_triangul(vglob(1)));
Vecteur * solEr = &(matglob.Simple_Resol_systID(vgloberr)); // résolution sans modification de la matrice de raideur
////---debug
// cout << "\n debug Algori::RemontErreur ";
// solEr->Affiche();
// cout << endl;
////---fin debug
// sol->Affiche();
// ------------- résultat transféré dans les noeuds ---------
// actualisation des ddl actifs a t
lesMail->PlusDelta_t(*solEr,Ass_err->Nb_cas_assemb());
// inactivation des ddl d'erreur
lesMail->Inactive_ddl_Erreur();
// fin des calculs
};