// 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 "CondLim.h"
#include "ConstMath.h"
#include "ParaGlob.h"
// ---------- variables globales ------------
// déclaration d'une variable de travail, utilisée par CondlineaireCHRepere
Mat_pleine CondLim::matinter(1,1);
// ---------- fin variables globales ------------
// contructeur
CondLim::CondLim() :
VImpSM(),VImpRaid()
{ };
// destructeur
CondLim::~CondLim()
{ };
// concernant les sous class
CondLim::ImpRaid::ImpRaid() : // part defaut
ligne(),colonne()
{};
CondLim::ImpRaid::ImpRaid(int i,const Vecteur& vL,const Vecteur& vC,double val) : // entree nomale
ilicol(i), ligne(vL),colonne(vC),valeur(val)
{};
CondLim::ImpRaid::ImpRaid(const ImpRaid & a) :
ilicol(a.ilicol), ligne(a.ligne),colonne(a.colonne),valeur(a.valeur)
{};
CondLim::ImpRaid& CondLim::ImpRaid::operator=(const ImpRaid& a) // assigment
{ this->ilicol = a.ilicol; this->ligne = a.ligne; this->colonne = a.colonne;
this->valeur = a.valeur;
return *this;
};
CondLim::ImpSM& CondLim::ImpSM::operator=(const ImpSM& a) // assigment
{ this->iligne = a.iligne; this->valeur = a.valeur;
return *this;
};
// valeur imposee au second membre
// vecglob : le second membre, i : la position globale du ddl impose
// val : la valeur a imposer
// vec2 : est un second vecteur éventuel (si != NULL) sur lequel on impose les mêmes CL que vecglob
// mais sans sauvegarde (correspond par exemple à une partie de vecglob)
void CondLim::Val_imposee_Sm(Vecteur& vecglob,int i,double val,Vecteur* vec2)
{ // on sauvegarde la valeur du second membre existante et on remplace
// par le ddl impose
ImpSM sauv(i,vecglob(i));
VImpSM.push_back(sauv);
vecglob(i) = val;
if (vec2 != NULL)
(*vec2)(i) = val;
};
// valeur imposee a la matrice et au second membre si la valeur est
// differente de zero
// matglob : la matrice, i : la position globale du ddl impose
// val : la valeur a imposer
// vec2 : est un second vecteur éventuel (si != NULL) sur lequel on impose les mêmes CL que vecglob
// mais sans sauvegarde (correspond par exemple à une partie de vecglob)
void CondLim::Val_imposee_Mat(Mat_abstraite & matglob,Vecteur& vecglob,
int i,double val,Vecteur* vec2)
{// on sauvegarde la ligne et la colonne de la matrice, correspondant
// au ddl impose
Vecteur vL = matglob.LigneSpe(i);
Vecteur vC = matglob.ColonneSpe(i);
ImpRaid sauv(i,vL,vC,vecglob(i)); // on sauvegarde la réaction
VImpRaid.push_back(sauv);
// on impose sur la matrice et on met les répercussions sur le second membre, mais
// on n'impose pas sur le second membre à ce stade, car en fait la valeur imposer
// pourrait éventuellement encore changer (sans doute non, car les coeff hors diag sont
// maintenant nul, mais on laisse quand même le système à deux temps)
if (val != 0.) // cas ou il faut modifier le second membre
{ vecglob -= val*matglob.Colonne(i); // recup de la colonne et modif
if (vec2 != NULL)
*vec2 -= val*matglob.Colonne(i);
// je crois que la ligne d'après est fausse, car là on applique la cl, alors que l'objet de la routine concerne uniquement la matrice
// et les répercussions sur le second membre. Si val est diff de 0, vecglob(i) va changer et la valeur sauvegardé ensuite par
// CondLim::Val_imposee_Sm( ne sera pas bonne
// vecglob(i) = val; // mais a jour le second membre
};
matglob.MetValLigne(i,0.); // mise de la ligne
matglob.MetValColonne(i,0.); // et de la colonne a zero
matglob(i,i) = 1.; // puis l'element diagonnal = 1
};
// cas particulier de valeur imposee a une matrice
// c'a-dire val sur la diagonale
// matglob : la matrice, i : la position globale du ddl impose
// dans ce cas-ci il n'y a aucune information sauvegardée, des valeurs modifiées
// cela signifie qu'après cette fonction, les appels aux routines pour la
// remontée aux réactions, n'ont aucun sens
void CondLim::Val_imposSimple_Mat(Mat_abstraite & matglob,int i,double val)
{matglob.MetValLigne(i,0.); // mise de la ligne
matglob.MetValColonne(i,0.); // et de la colonne a zero
matglob(i,i) = val; // puis l'element diagonnal = 1
};
// remontee aux efforts apres resolution
// ceci pour la ligne i, dans le cas d'un ddl bloque
double CondLim::RemonteDdlBloqueMat(Mat_abstraite & matglob,Vecteur& solution, int iligne)
{ // on boucle sur le container pour retrouver le ddl
bool cherche = false;
list <ImpRaid>::iterator ip;
for (ip = VImpRaid.begin();ip != VImpRaid.end(); ip++)
if ((*ip).ilicol == iligne)
{ matglob.RemplaceLigneSpe(iligne,(*ip).ligne);
return matglob.Prod_Ligne_vec (iligne,solution);
cherche = true;
};
if (cherche != true)
{ cout << "\n par d effort sur ce ddl nb = " << iligne;
cout << "\n double CondLim::RemonteDdlBloqueMat ( etc ..." << endl;
};
return 0.;
};
// retourne le maxi des efforts exterieurs
double CondLim::MaxEffort(int & ili)
{ double res = 0;
ili = 1;
list <ImpSM>::iterator ip;
for (ip = VImpSM.begin();ip != VImpSM.end(); ip++)
if (Abs((*ip).valeur) > res)
{res = Abs((*ip).valeur);
ili = (*ip).iligne;
}
return res;
};
// retourne la valeur initiale au second membre avant condition limite
// en fonction du pointeur d'assemblage
double CondLim::ValReact(int & ili)
{ list <ImpSM>::iterator ip;
for (ip = VImpSM.begin();ip != VImpSM.end(); ip++)
if ((*ip).iligne == ili )
{ return (*ip).valeur;
}
// on n'a rien trouver
cout << "\n erreur la ligne demandee ne correspond pas a un ddl bloque ";
cout << "\n CondLim::ValReact(int & ili) ";
Sortie (1);
return 0.; // pour pas avoir de message a la compile
};
//----- cas de conditions limites lineaires entre plusieurs ddl-----
// class qui permet d'enregistrer le changement de repere
CondLim::lineaires::lineaires () : // par defaut
val(),pt()
{ indice = 0; valeur = 0.; };
CondLim::lineaires::lineaires (int i,const Vecteur& vL,const Tableau<int> & ptt,double vale):
indice(i),val(vL),pt(ptt), valeur(vale)
{};
CondLim::lineaires::lineaires (const lineaires & a) : // constructeur de copie
indice(a.indice),val(a.val),pt(a.pt),valeur(a.valeur)
{};
CondLim::lineaires& CondLim::lineaires::operator=(const lineaires& a) // assigment
{ this->indice = a.indice; this->val = a.val;this->pt = a.pt;this->valeur = a.valeur;
return *this; };
// premier temps : on prépare la mise en place de la condition par un changement
// de repère sur la matrice et sur le second membre
// pt : tableau des pointeurs de ddl concerne, pt(i) = la position du ddl i
// dans la matrice globale
// val : tableau des coefficients de la condition lineaire
// valeur : valeur a laquelle est egale la condition lineaire
// cond lineaire -> somme des val(i) * ddl(pt(i)) = valeur
// la procedure modifie les reperes d'expression des ddl, mais sauvegarde
// les infos permettant de reconstruire les reperes initiaux
// vec2 : est un second vecteur éventuel (si != NULL) sur lequel on impose les mêmes CL que vecglob
// mais sans sauvegarde (correspond par exemple à une partie de vecglob)
void CondLim::CondlineaireCHRepere(Mat_abstraite & matglob,Vecteur& vecglob,
const Tableau<int> & pt,const Vecteur& val, double valeur,Vecteur* vec2)
{
// //--- debug
// cout << "\n CondLim::CondlineaireCHRepere ";
// cout << "\n valeur = " << valeur << " val= " << val << " pt= " << pt << endl ;
// //--- fin debug
// on commence par normer la condition limite
int tail = pt.Taille();
double norme = val.Norme();
if ( norme <= ConstMath::trespetit )
{ cout << "\n erreur , la condition lineaire est formee de coefficients nulles !! ";
cout << "\n CondLim::Condlineaire( etc ... " << endl;
Sortie(1);
};
// norme = 1. ; // pour le débug !!!
// ==->par construction le premier indice correspond a la direction a bloquer
// --ici dans la normalisation je pense qu'il y a un pb donc on supprime momentanément
// double unsurnorme = 1. / norme ;
// valeur *= unsurnorme;
// --
//norme = 1;
// ancienne méthode : calcul du vecteur V
// val *= unsurnorme ;
// Vecteur V = -val/val(1);
// V(1) = 1./val(1);
// nouvelle méthode : calcul du vecteur V
Vecteur V(val);
V /= (-val(1));
V(1) = 1 / val(1);
// V(1) = norme / val(1);
// modif du vecteur en Vprime
V(1) -= 1.;
// calcul de la matrice de raideur modifiee: pour l'instant ne fonctionne que pour des matrices carrees ou bandes
// //--- debug
// cout << "\n CondLim::CondlineaireCHRepere ";
// cout << "\n V= " << V << endl ;
// //--- fin debug
#ifdef MISE_AU_POINT
if ((matglob.Type_matrice() != CARREE) && (matglob.Type_matrice() != CARREE_SYMETRIQUE)
&& (matglob.Type_matrice() != BANDE_SYMETRIQUE) && (matglob.Type_matrice() != BANDE_NON_SYMETRIQUE)
&& (matglob.Type_matrice() != BANDE_NON_SYMETRIQUE_LAPACK)
&& (matglob.Type_matrice() != BANDE_SYMETRIQUE_LAPACK)
&& (matglob.Type_matrice() != CARREE_LAPACK) && (matglob.Type_matrice() != CARREE_SYMETRIQUE_LAPACK))
{ cout << "\n erreur: pour l'instant la mise en place de condition lineaire par rotation n'est implante "
<< " que pour les matrices carree ou bande, symetrique ou pas , alors que la matrice ici est "
<< matglob.Type_matrice() << " ";
};
#endif
// la condition est la suivante sur la matrice de raideur:
// kp(i,j) = k(i,j) + Cp(ip) * K(a, j) + K(i,a) * Cp(jp) + Cp(ip) K(a,a) Cp(jp)
// pour tout Cp(ip) non nul, on doit balayer tous les j et idemn pour jp non nul, balayer les i
// mais!!! il faut tenir compte du type de stockage, à savoir certain K(i,j) peuvent ne pas être licite
// car ils n'existent pas (cas en dehors de la bande par exemple pour un stockage bande)
// a est l'indice où est imposé la condition limite, Kp est la nouvelle valeur de la matrice de raideur
// Cp vaut C sauf pour la première composante qui vaut : Cp(1) = C(1) -1
// on commence par une matrice carrée pour les tests pour faciliter la mise au point
// ici tous les K(i,j) existent a priori
if ((matglob.Type_matrice() == CARREE)|| (matglob.Type_matrice() == CARREE_LAPACK))
{ int nbligne = matglob.Nb_ligne();
int nbcol = matglob.Nb_colonne();
int ic = pt(1); // indice du ddl imposé
//---- pour débug
// // on symétrise la matrice
// for (int ix = 1; ix <= nbligne; ix++)
// for (int iy=1;iy <= ix; iy++)
// matglob(ix,iy) = matglob(iy,ix) = 0.5*(matglob(ix,iy)+matglob(iy,ix));
// // affiche les termes non symétriques de la matrice s'il y en a
// cout << "\n disymetrie eventuelle avant CLL ";
// matglob.Symetrie();
//---- fin débug
// on récupère la ligne et la colonne qui passe par ce point
Vecteur Bpj = matglob.Ligne(ic);
Vecteur Api = matglob.Colonne(ic);
////--- debug
// cout << "\n CondLim::CondlineaireCHRepere ";
// cout << "\n Bpj = "; Bpj.Affiche();
// cout << "\n Api= "; Api.Affiche();
// cout << endl ;
////--- fin debug
double matglobii = matglob(ic,ic); // pour simplifier
// calcul des termes modiés de la matrice
// 1) termes calculées à partir du terme diagonal matglobii: K_4
for (int ilc =1; ilc <= tail; ilc++)
for (int jlc =1; jlc <= tail; jlc++)
{int l = pt(ilc);int m = pt(jlc);
matglob(l,m) += V(ilc) * matglobii * V(jlc);
};
// 2) termes non diagonaux: K_2
for (int col=1;col <= tail;col++)
for (int ix = 1; ix <= nbligne; ix++)
{ int m = pt(col);
matglob (ix,m) += Api(ix) * V(col);
};
// 3) termes non diagonaux: K_3
for (int iy=1;iy <= nbcol;iy++)
for (int lig=1; lig <= tail; lig++)
{ int l = pt(lig);
matglob (l,iy) += V(lig) * Bpj(iy);
};
//---- pour débug
// // affiche les termes non symétriques de la matrice s'il y en a
// cout << "\n disymetrie eventuelle apres CLL ";
// matglob.Symetrie();
//---- fin débug
}
else
//--- cas d'une matrice quelconque ----
{ int nbligne = matglob.Nb_ligne();
int ic = pt(1); // indice du ddl imposé
// on récupère la ligne et la colonne qui passe par ce point
Vecteur Bpj = matglob.Ligne(ic);
Vecteur Api = matglob.Colonne(ic);
double matglobii = matglob(ic,ic); // pour simplifier
////--- debug
// cout << "\n CondLim::CondlineaireCHRepere ";
// cout << "\n matglobii = "<< matglobii ;
// cout << "\n Bpj = "; Bpj.Affiche();
// cout << "\n Api= "; Api.Affiche();
// cout << endl ;
////--- fin debug
//
// calcul des termes modifiés de la matrice
// if (!matglob.Symetrie())
if (!Symetrique_Enum_matrice(matglob.Type_matrice()))
{ //$$ cas d'une matrice non symétrique
// 1) termes calculées à partir du terme diagonal matglobii: K_4
for (int ilc =1; ilc <= tail; ilc++)
for (int jlc =1; jlc <= tail; jlc++)
{int l = pt(ilc);int m = pt(jlc);
matglob(l,m) += V(ilc) * matglobii * V(jlc);
};
// 2) deuxième et troisième terme
for (int col=1;col <= tail;col++)
for (int ix = 1; ix <= nbligne; ix++) // on sait que les termes non nuls sont ceux compris dans lb
{ int m = pt(col);
double r = Api(ix) * V(col);
// normalement il n'y aura de terme non null qu'au niveau de place existant dans la matrice !!
// ce qui permet d'éviter de tester l'appartenance ou non à la bande par exemple pour une matrice bande
// NON!! ce n'est pas vrai, il faut vérifier que l'on ne sort pas des clous
// en effet on fait varier ix de 1 à nbligne donc toute la matrice !! donc forcément on sortira de la bande
// si c'est un stockage bande
// on ne fait le calcul que si le terme de la matrice K existe
if (matglob.Existe(ix,m))
{if (Dabs(r) > ConstMath::trespetit)
matglob (ix,m) += r;
};
// on fait l'inverse: au niveau des indices, pour éviter une nouvelle double boucle
r = V(col) * Bpj(ix);
if (matglob.Existe(m,ix))
{if (Dabs(r) > ConstMath::trespetit)
matglob (m,ix) += r;
};
};
}
else //$$ cas d'une matrice symétrique
{ // 1) termes calculées à partir du terme diagonal matglobii: K_4
for (int ilc =1; ilc <= tail; ilc++)
for (int jlc =1; jlc <= tail; jlc++)
{int l = pt(ilc);int m = pt(jlc);
if (m >= l)
matglob(l,m) += V(ilc) * matglobii * V(jlc);
};
// 2) deuxième et troisième terme
for (int col=1;col <= tail;col++)
for (int ix = 1; ix <= nbligne; ix++) // on sait que les termes non nuls sont ceux compris dans lb
{ int m = pt(col);
double r = Api(ix) * V(col);
// normalement il n'y aura de terme non null qu'au niveau de place existant dans la matrice !!
// ce qui permet d'éviter de tester l'appartenance ou non à la bande par exemple pour une matrice bande
// NON!! ce n'est pas vrai, il faut vérifier que l'on ne sort pas des clous
// en effet on fait varier ix de 1 à nbligne donc toute la matrice !! donc forcément on sortira de la bande
// si c'est un stockage bande
// on ne fait le calcul que si le terme de la matrice K existe
if (matglob.Existe(ix,m))
{ if ((Dabs(r) > ConstMath::trespetit) && (m >= ix ))
matglob (ix,m) += r;
};
// on fait l'inverse: au niveau des indices, pour éviter une nouvelle double boucle
r = V(col) * Bpj(ix);
if (matglob.Existe(m,ix))
{if ((Dabs(r) > ConstMath::trespetit) && (ix >= m))
matglob (m,ix) += r;
};
};
////--- debug
//{
// Vecteur Bpj_ = matglob.Ligne(ic);
// Vecteur Api_ = matglob.Colonne(ic);
// double matglobii_ = matglob(ic,ic); // pour simplifier
// cout << "\n CondLim::CondlineaireCHRepere ";
// cout << "\n matglobii = "<< matglobii_ ;
// cout << "\n Bpj_ = "; Bpj_.Affiche();
// cout << "\n Api_= "; Api_.Affiche();
// cout << endl ;
// }
////--- fin debug
};
};
// calcul du second membre modifie, et éventuellement d'un second vecteur
int ili = pt(1);
if (vec2 == NULL)
{double ri = vecglob(ili);
for (int lig=1; lig <= tail; lig++)
{int l = pt(lig);
vecglob (l) += ri * V(lig);
};
}
else // cas où on applique également au second vecteur
{double ri = vecglob(ili);
double ri2 = (*vec2)(ili);
for (int lig=1; lig <= tail; lig++)
{ int l = pt(lig);
vecglob (l) += ri * V(lig);
(*vec2) (l) += ri2 * V(lig);
};
};
// retour du Vprime en V et sauvegarde
V(1) += 1.;
lineaires sauv(1,V,pt,valeur);
// on sauvegarde la condition lineaire au début de la liste de manière ensuite si l'on
// parcourt la liste de commencer par la dernière entrée et de remonter vers les plus anciennes
Vlineaires.push_front(sauv);
};
// second temps : imposition des blocages correspondant aux conditions lineaires
// vec2 : est un second vecteur éventuel (si != NULL) sur lequel on impose les mêmes CL que vecglob
// mais sans sauvegarde (correspond par exemple à une partie de vecglob)
void CondLim::CondlineaireImpose (Mat_abstraite & matglob,Vecteur& vecglob,Vecteur* vec2)
{ list <lineaires>::iterator ip,ipfin = Vlineaires.end();
if (Vlineaires.size() != 0)
for (ip = Vlineaires.begin();ip != ipfin; ip++)
{ lineaires& lin = (*ip);
int ptassemb = lin.pt(lin.indice); // pointeur d'assemblage globale
// on impose la condition limite a la raideur et au second membre
// NB: a priori ici la valeur imposée est toujours nulle (vue la méthode employée
// pour tenir compte des conditions linéaires, cf. théorie)
// on impose sur la matrice et on met les répercussions sur le second membre,
// mais on n'impose pas sur le second membre à ce stade, car en fait la valeur imposer
// pourrait éventuellement encore changer (sans doute non, car les coeff hors diag sont
// maintenant nul, mais on laisse quand même le système à deux temps)
Val_imposee_Mat(matglob,vecglob,ptassemb,lin.valeur,vec2);
// là on impose définitivement les valeurs au second membre
Val_imposee_Sm(vecglob,ptassemb,lin.valeur,vec2);
}
// affichage éventuelle de la matrice de raideur et du second membre
if (ParaGlob::NiveauImpression() >= 10)
{ string entete = " affichage de la matrice de raideur (puissance interne) apres blocage dus aux ";
entete += " conditions lineaires limites";
matglob.Affichage_ecran(entete);
entete = " affichage du second membre (puissance interne) apres blocage dus aux ";
entete += " conditions lineaires limites";
vecglob.Affichage_ecran(entete);
}
};
//expression du vecteur resultat dans les reperes initiaux
// sol : la solution, est modifie et retournee dans les reperes initiaux
void CondLim::RepInitiaux( Vecteur& sol)
{ // on balaie les conditions limites lineaires
// ceux-ci sont ordonnées de manière à ce que la première conditions limites est la
// dernière rentrée
list <lineaires>::iterator ip,ipfin= Vlineaires.end();
//-- debug
// Vecteur intersol(Vlineaires.size());
// int num = 1;
//-- fin debug
// on revient par rotation inverse
for (ip = Vlineaires.begin();ip != ipfin; ip++)
{ Vecteur& val = (*ip).val;
int indice = (*ip).indice;
Tableau<int> & pt = (*ip).pt;
int tail = val.Taille();
double r = 0;
//-- debug
//cout << "\n pt(indice)= " << pt(indice) ;
//-- fin debug
for (int lig=1; lig <= tail; lig++)
{ r += val(lig) * sol (pt(lig));
//-- debug
//if (val(lig) != 0.)
// cout << " val(" << lig << ")= " << val(lig) << " sol (" << pt(lig) << ")= " << sol (pt(lig)) ;
//-- fin debug
};
//-- debug
//cout << " r= " << r << endl;
//-- fin debug
sol(pt(indice)) = r; // à commenter pour le débug
// intersol(num) = r; // debug
// num++; // debug
};
//-- debug
// num = 1;
// for (ip = Vlineaires.begin();ip != ipfin; ip++)
// { int indice = (*ip).indice;
// Tableau<int> & pt = (*ip).pt;
// sol(pt(indice)) = intersol(num);
// }
//-- fin debug
// affichage éventuelle du vecteur solution : deltat_ddl
if (Vlineaires.size() != 0)
if (ParaGlob::NiveauImpression() >= 10)
{ string entete = " affichage du vecteur solution (delta_ddl) apres retour dans les reperes initiaux";
sol.Affichage_ecran(entete); };
};
// application d'une condition linéaire seule, avec en retour, la situation de la condition linéaire
// imposée, ramené dans le repère initial
void CondLim::CondiLineaireImposeComplet(Mat_abstraite & matglob,Vecteur& vecglob,
const Tableau<int> & pt,const Vecteur& val, double valeur,Vecteur* vec2)
{
// //--- debug
// cout << "\n CondLim::CondiLineaireImposeComplet ";
// // pour tester la rotation
// cout << "\n valeur = " << valeur << " val= " << val << " pt= " << pt << endl ;
// cout << "\n matrice au début: ";
// matglob.Affiche2(1,16,1,1,16,1,14);
// cout << "\n second membre au début: ";vecglob.Affiche();
//{ Vecteur vol(16);
// vol(1) = 0.1; vol(2)=0.;
// vol(3) = 0.1; vol(4)=0.015;
// vol(5) = 0.05; vol(6)=0.;
// vol(7) = 0.05; vol(8)=0.015;
// vol(9) = 0.05; vol(10)=0.;
// vol(11) = 0.05; vol(12)=0.015;
// vol(13) = 0.; vol(14)=0.;
// vol(15) = 0.; vol(16)=0.015;
// Vecteur toto = matglob.Prod_mat_vec(vol);
// cout << "\n résidu matglob*deplacement: ";toto.Affiche();
//}
//
// //--- fin debug
// I) ========= première partie : on commence par changer de repère pour l'application de la condition limite
// on commence par normer la condition limite
// int tail = pt.Taille();
double norme = val.Norme();
if ( norme <= ConstMath::trespetit )
{ cout << "\n erreur , la condition lineaire est formee de coefficients nulles !! ";
cout << "\n CondLim::Condlineaire( etc ... " << endl;
Sortie(1);
};
Vecteur V(val);
V /= (-val(1));
V(1) = 1 / val(1);
// calcul de la matrice de raideur modifiee: pour l'instant ne fonctionne que pour des matrices carrees ou bandes
Rotation(matglob,vecglob,pt,V,vec2);
// on sauvegarde la condition lineaire au début de la liste de manière ensuite si l'on
// parcourt la liste de commencer par la dernière entrée et de remonter vers les plus anciennes
lineaires sauv(1,V,pt,valeur);
Vlineaires.push_front(sauv);
// II) ========= deuxième partie : on applique la condition fixée, due à la condition limite
// //--- debug
// cout << "\n CondLim::CondiLineaireImposeComplet ";
// cout << "\n matrice après retation et avant application du ddl bloqué : ";
// matglob.Affiche2(1,16,1,1,16,1,14);
// cout << "\n second membre après retation et avant application du ddl bloqué: ";vecglob.Affiche();
// //--- fin debug
int ptassemb = pt(sauv.indice);
// on impose la condition limite a la raideur et au second membre
// NB: a priori ici la valeur imposée est toujours nulle (vue la méthode employée
// pour tenir compte des conditions linéaires, cf. théorie)
// on impose sur la matrice et on met les répercussions sur le second membre,
// mais on n'impose pas sur le second membre à ce stade, car en fait la valeur imposer
// pourrait éventuellement encore changer (sans doute non, car les coeff hors diag sont
// maintenant nul, mais on laisse quand même le système à deux temps)
Val_imposee_Mat(matglob,vecglob,ptassemb,sauv.valeur,vec2);
// là on impose définitivement les valeurs au second membre
Val_imposee_Sm(vecglob,ptassemb,sauv.valeur,vec2);
// //--- debug
// cout << "\n CondLim::CondiLineaireImposeComplet ";
// cout << "\n matrice après rotation et après application du ddl bloqué : ";
// matglob.Affiche2(1,16,1,1,16,1,14);
// cout << "\n second membre après rotation et après application du ddl bloqué: ";vecglob.Affiche();
// //--- fin debug
// III) ========= troisième partie : on revient dans le repère initiale
// calcul de la matrice de raideur modifiee: pour l'instant ne fonctionne que pour des matrices carrees ou bandes
V = val; // récup de la condition initiale
Rotation(matglob,vecglob,pt,V,vec2);
// là on impose définitivement les valeurs au second membre
// Val_imposee_Sm(vecglob,ptassemb,sauv.valeur,vec2);
// //--- debug
// cout << "\n CondLim::CondiLineaireImposeComplet ";
// cout << "\n matrice à la fin: ";
// matglob.Affiche2(1,16,1,1,16,1,14);
// cout << "\n second membre à la fin: ";vecglob.Affiche();
//{ Vecteur vol(16);
// vol(1) = 0.1; vol(2)=0.;
// vol(3) = 0.1; vol(4)=0.015;
// vol(5) = 0.05; vol(6)=0.;
// vol(7) = 0.05; vol(8)=0.015;
// vol(9) = 0.05; vol(10)=0.;
// vol(11) = 0.05; vol(12)=0.015;
// vol(13) = 0.; vol(14)=0.;
// vol(15) = 0.; vol(16)=0.015;
// Vecteur toto = matglob.Prod_mat_vec(vol);
// cout << "\n résidu matglob*deplacement: ";toto.Affiche();
//}
// Sortie(1);
//--- fin debug
};
// remise a zero des sauvegardes de second membre et de matrice
void CondLim::EffaceSauvegarde()
{ // effacement des deux containers
VImpSM.erase(VImpSM.begin(),VImpSM.end());
VImpRaid.erase(VImpRaid.begin(),VImpRaid.end());
};
// remise a zero des sauvegardes de condition lineaire
void CondLim::EffaceCoLin()
{ Vlineaires.erase(Vlineaires.begin(),Vlineaires.end());
};
//---------- méthodes internes --------
// méthode pour faire une rotation de la matrice et du second membre suivant un vecteur
// vec2 : est un second vecteur éventuel (si != NULL)
void CondLim::Rotation(Mat_abstraite & matglob,Vecteur& vecglob,
const Tableau<int> & pt,const Vecteur& val,Vecteur* vec2)
{
Vecteur V(val);
int tail = pt.Taille();
V(1) -= 1.;
// calcul de la matrice de raideur modifiee: pour l'instant ne fonctionne que pour des matrices carrees ou bandes
#ifdef MISE_AU_POINT
if ((matglob.Type_matrice() != CARREE) && (matglob.Type_matrice() != CARREE_SYMETRIQUE)
&& (matglob.Type_matrice() != BANDE_SYMETRIQUE) && (matglob.Type_matrice() != BANDE_NON_SYMETRIQUE)
&& (matglob.Type_matrice() != BANDE_NON_SYMETRIQUE_LAPACK)
&& (matglob.Type_matrice() != BANDE_SYMETRIQUE_LAPACK)
&& (matglob.Type_matrice() != CARREE_LAPACK) && (matglob.Type_matrice() != CARREE_SYMETRIQUE_LAPACK))
{ cout << "\n erreur: pour l'instant la mise en place de condition lineaire par rotation n'est implante "
<< " que pour les matrices carree ou bande, symetrique ou pas , alors que la matrice ici est "
<< matglob.Type_matrice() << " ";
};
#endif
// la condition est la suivante sur la matrice de raideur:
// kp(i,j) = k(i,j) + Cp(ip) * K(a, j) + K(i,a) * Cp(jp) + Cp(ip) K(a,a) Cp(jp)
// pour tout Cp(ip) non nul, on doit balayer tous les j et idemn pour jp non nul, balayer les i
// mais!!! il faut tenir compte du type de stockage, à savoir certain K(i,j) peuvent ne pas être licite
// car ils n'existent pas (cas en dehors de la bande par exemple pour un stockage bande)
// a est l'indice où est imposé la condition limite, Kp est la nouvelle valeur de la matrice de raideur
// Cp vaut C sauf pour la première composante qui vaut : Cp(1) = C(1) -1
// on commence par une matrice carrée pour les tests pour faciliter la mise au point
// ici tous les K(i,j) existent a priori
if ((matglob.Type_matrice() == CARREE)|| (matglob.Type_matrice() == CARREE_LAPACK))
{ int nbligne = matglob.Nb_ligne();
int nbcol = matglob.Nb_colonne();
int ic = pt(1); // indice du ddl imposé
//---- pour débug
// // on symétrise la matrice
// for (int ix = 1; ix <= nbligne; ix++)
// for (int iy=1;iy <= ix; iy++)
// matglob(ix,iy) = matglob(iy,ix) = 0.5*(matglob(ix,iy)+matglob(iy,ix));
// // affiche les termes non symétriques de la matrice s'il y en a
// cout << "\n disymetrie eventuelle avant CLL ";
// matglob.Symetrie();
//---- fin débug
// on récupère la ligne et la colonne qui passe par ce point
Vecteur Bpj = matglob.Ligne(ic);
Vecteur Api = matglob.Colonne(ic);
////--- debug
// cout << "\n CondLim::CondlineaireCHRepere ";
// cout << "\n Bpj = "; Bpj.Affiche();
// cout << "\n Api= "; Api.Affiche();
// cout << endl ;
////--- fin debug
double matglobii = matglob(ic,ic); // pour simplifier
// calcul des termes modiés de la matrice
// 1) termes calculées à partir du terme diagonal matglobii: K_4
for (int ilc =1; ilc <= tail; ilc++)
for (int jlc =1; jlc <= tail; jlc++)
{int l = pt(ilc);int m = pt(jlc);
matglob(l,m) += V(ilc) * matglobii * V(jlc);
};
// 2) termes non diagonaux: K_2
for (int col=1;col <= tail;col++)
for (int ix = 1; ix <= nbligne; ix++)
{ int m = pt(col);
matglob (ix,m) += Api(ix) * V(col);
};
// 3) termes non diagonaux: K_3
for (int iy=1;iy <= nbcol;iy++)
for (int lig=1; lig <= tail; lig++)
{ int l = pt(lig);
matglob (l,iy) += V(lig) * Bpj(iy);
};
//---- pour débug
// // affiche les termes non symétriques de la matrice s'il y en a
// cout << "\n disymetrie eventuelle apres CLL ";
// matglob.Symetrie();
//---- fin débug
}
else
//--- cas d'une matrice quelconque ----
{ int nbligne = matglob.Nb_ligne();
int ic = pt(1); // indice du ddl imposé
// on récupère la ligne et la colonne qui passe par ce point
Vecteur Bpj = matglob.Ligne(ic);
Vecteur Api = matglob.Colonne(ic);
double matglobii = matglob(ic,ic); // pour simplifier
////--- debug
// cout << "\n CondLim::CondlineaireCHRepere ";
// cout << "\n matglobii = "<< matglobii ;
// cout << "\n Bpj = "; Bpj.Affiche();
// cout << "\n Api= "; Api.Affiche();
// cout << endl ;
////--- fin debug
//
// calcul des termes modifiés de la matrice
// if (!matglob.Symetrie())
if (!Symetrique_Enum_matrice(matglob.Type_matrice()))
{ //$$ cas d'une matrice non symétrique
// 1) termes calculées à partir du terme diagonal matglobii: K_4
for (int ilc =1; ilc <= tail; ilc++)
for (int jlc =1; jlc <= tail; jlc++)
{int l = pt(ilc);int m = pt(jlc);
matglob(l,m) += V(ilc) * matglobii * V(jlc);
};
// 2) deuxième et troisième terme
for (int col=1;col <= tail;col++)
for (int ix = 1; ix <= nbligne; ix++) // on sait que les termes non nuls sont ceux compris dans lb
{ int m = pt(col);
double r = Api(ix) * V(col);
// normalement il n'y aura de terme non null qu'au niveau de place existant dans la matrice !!
// ce qui permet d'éviter de tester l'appartenance ou non à la bande par exemple pour une matrice bande
// NON!! ce n'est pas vrai, il faut vérifier que l'on ne sort pas des clous
// en effet on fait varier ix de 1 à nbligne donc toute la matrice !! donc forcément on sortira de la bande
// si c'est un stockage bande
// on ne fait le calcul que si le terme de la matrice K existe
if (matglob.Existe(ix,m))
{if (Dabs(r) > ConstMath::trespetit)
matglob (ix,m) += r;
};
// on fait l'inverse: au niveau des indices, pour éviter une nouvelle double boucle
r = V(col) * Bpj(ix);
if (matglob.Existe(m,ix))
{if (Dabs(r) > ConstMath::trespetit)
matglob (m,ix) += r;
};
};
}
else //$$ cas d'une matrice symétrique
{ // 1) termes calculées à partir du terme diagonal matglobii: K_4
for (int ilc =1; ilc <= tail; ilc++)
for (int jlc =1; jlc <= tail; jlc++)
{int l = pt(ilc);int m = pt(jlc);
if (m >= l)
matglob(l,m) += V(ilc) * matglobii * V(jlc);
};
// 2) deuxième et troisième terme
for (int col=1;col <= tail;col++)
for (int ix = 1; ix <= nbligne; ix++) // on sait que les termes non nuls sont ceux compris dans lb
{ int m = pt(col);
double r = Api(ix) * V(col);
// normalement il n'y aura de terme non null qu'au niveau de place existant dans la matrice !!
// ce qui permet d'éviter de tester l'appartenance ou non à la bande par exemple pour une matrice bande
// NON!! ce n'est pas vrai, il faut vérifier que l'on ne sort pas des clous
// en effet on fait varier ix de 1 à nbligne donc toute la matrice !! donc forcément on sortira de la bande
// si c'est un stockage bande
// on ne fait le calcul que si le terme de la matrice K existe
if (matglob.Existe(ix,m))
{ if ((Dabs(r) > ConstMath::trespetit) && (m >= ix ))
matglob (ix,m) += r;
};
// on fait l'inverse: au niveau des indices, pour éviter une nouvelle double boucle
r = V(col) * Bpj(ix);
if (matglob.Existe(m,ix))
{if ((Dabs(r) > ConstMath::trespetit) && (ix >= m))
matglob (m,ix) += r;
};
};
////--- debug
//{
// Vecteur Bpj_ = matglob.Ligne(ic);
// Vecteur Api_ = matglob.Colonne(ic);
// double matglobii_ = matglob(ic,ic); // pour simplifier
// cout << "\n CondLim::CondlineaireCHRepere ";
// cout << "\n matglobii = "<< matglobii_ ;
// cout << "\n Bpj_ = "; Bpj_.Affiche();
// cout << "\n Api_= "; Api_.Affiche();
// cout << endl ;
// }
////--- fin debug
};
};
// calcul du second membre modifie, et éventuellement d'un second vecteur
int ili = pt(1);
if (vec2 == NULL)
{double ri = vecglob(ili);
for (int lig=1; lig <= tail; lig++)
{int l = pt(lig);
vecglob (l) += ri * V(lig);
};
}
else // cas où on applique également au second vecteur
{double ri = vecglob(ili);
double ri2 = (*vec2)(ili);
for (int lig=1; lig <= tail; lig++)
{ int l = pt(lig);
vecglob (l) += ri * V(lig);
(*vec2) (l) += ri2 * V(lig);
};
};
// retour à la situation initiale
V(1) += 1.;
};