432 lines
19 KiB
C++
Executable file
432 lines
19 KiB
C++
Executable file
// This file is part of the Herezh++ application.
|
|
//
|
|
// The finite element software Herezh++ is dedicated to the field
|
|
// of mechanics for large transformations of solid structures.
|
|
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
|
|
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
|
|
//
|
|
// Herezh++ is distributed under GPL 3 license ou ultérieure.
|
|
//
|
|
// Copyright (C) 1997-2021 Université Bretagne Sud (France)
|
|
// AUTHOR : Gérard Rio
|
|
// E-MAIL : gerardrio56@free.fr
|
|
//
|
|
// This program is free software: you can redistribute it and/or modify
|
|
// it under the terms of the GNU General Public License as published by
|
|
// the Free Software Foundation, either version 3 of the License,
|
|
// or (at your option) any later version.
|
|
//
|
|
// This program is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty
|
|
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
|
// See the GNU General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU General Public License
|
|
// along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
//
|
|
// For more information, please consult: <https://herezh.irdl.fr/>.
|
|
|
|
#include "Algo_zero.h"
|
|
|
|
// contiend que les méthodes qui ne sont pas des templates
|
|
|
|
#include "MathUtil.h"
|
|
#include "ConstMath.h"
|
|
|
|
|
|
// CONSTRUCTEURS :
|
|
// constructeur par défaut
|
|
Algo_zero::Algo_zero() :
|
|
prec_res_abs(1.e-3) // précision absolue par défaut sur le résidu
|
|
,prec_res_rel(1.e-3) // précision relative par défaut sur le résidu
|
|
,iter_max(20) // nombre d'itération maxi pour un incrément
|
|
,coef_mini_delta_x(ConstMath::pasmalpetit) // précision relative sur le mini de l'incrément de x
|
|
,mini_delta_x( ConstMath::trespetit) // minimum de delta x permis
|
|
,maxi_delta_x(ConstMath::pasmalgrand) // maximum de delta x permis
|
|
,nbMaxiIncre(20) // nombre maxi d'incrément permis
|
|
,maxi_test_moins1(2)
|
|
,permet_affichage(0)
|
|
{}
|
|
// constructeur de copie
|
|
Algo_zero::Algo_zero(const Algo_zero& a) :
|
|
prec_res_abs(a.prec_res_abs),prec_res_rel(a.prec_res_rel)
|
|
,iter_max(a.iter_max),coef_mini_delta_x(a.coef_mini_delta_x)
|
|
,mini_delta_x(a.mini_delta_x),maxi_delta_x(a.maxi_delta_x),nbMaxiIncre(a.nbMaxiIncre)
|
|
,maxi_test_moins1(a.maxi_test_moins1)
|
|
,permet_affichage(a.permet_affichage)
|
|
{}
|
|
|
|
// DESTRUCTEUR :
|
|
Algo_zero::~Algo_zero()
|
|
{}
|
|
|
|
// METHODES PUBLIQUES :
|
|
// initialisation de tous les paramètres à leurs valeurs par défaut
|
|
void Algo_zero::Init_param_val_defaut()
|
|
{ prec_res_abs = 1.e-3 ; // précision absolue par défaut sur le résidu
|
|
prec_res_rel = 1.e-3 ; // précision relative par défaut sur le résidu
|
|
iter_max = 20; // nombre d'itération maxi pour un incrément
|
|
coef_mini_delta_x = ConstMath::pasmalpetit; // précision relative sur le mini de l'incrément de x
|
|
mini_delta_x = ConstMath::trespetit; // minimum de delta x permis
|
|
maxi_delta_x = ConstMath::pasmalgrand ; // maximum de delta x permis
|
|
nbMaxiIncre = 20; // nombre maxi d'incrément permis
|
|
maxi_test_moins1 = 2;
|
|
permet_affichage = 0; // niveau par défaut d'affichage
|
|
}
|
|
|
|
|
|
// affichage à l'écran des infos
|
|
void Algo_zero::Affiche() const
|
|
{ cout << "\n Algo_zero:_parametres: "
|
|
<< " prec_res_abs= " << prec_res_abs
|
|
<< " prec_res_rel= " << prec_res_rel
|
|
<< " iter_max= " << iter_max
|
|
<< " coef_mini_delta_x= " << coef_mini_delta_x
|
|
<< " mini_delta_x= " << mini_delta_x
|
|
<< " maxi_delta_x= " << maxi_delta_x
|
|
<< " nbMaxiIncre= " << nbMaxiIncre
|
|
<< " maxi_test_moins1= "<< maxi_test_moins1
|
|
<< " permet_affichage " << permet_affichage
|
|
<< endl;
|
|
};
|
|
|
|
//----- lecture écriture de restart -----
|
|
// cas donne le niveau de la récupération
|
|
// = 1 : on récupère tout
|
|
// = 2 : on récupère uniquement les données variables (supposées comme telles)
|
|
void Algo_zero::Lecture_base_info (ifstream& ent,const int cas)
|
|
{ if (cas == 1)
|
|
{string nom;
|
|
ent >> nom;
|
|
if (nom != "Algo_zero:_para:")
|
|
{ cout << "\n erreur en lecture des parametres d'Algo_zero, on devait lire le mot clef "
|
|
<< "Algo_zero:_para: et on a lu "<< nom << endl;
|
|
Sortie(1);
|
|
};
|
|
ent >> nom >> prec_res_abs
|
|
>> nom >> prec_res_rel
|
|
>> nom >> iter_max
|
|
>> nom >> coef_mini_delta_x
|
|
>> nom >> mini_delta_x
|
|
>> nom >> maxi_delta_x
|
|
>> nom >> nbMaxiIncre
|
|
>> nom >> maxi_test_moins1;
|
|
ent >> nom >> permet_affichage;
|
|
};
|
|
};
|
|
|
|
|
|
// cas donne le niveau de sauvegarde
|
|
// = 1 : on sauvegarde tout
|
|
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
|
|
void Algo_zero::Ecriture_base_info(ofstream& sort,const int cas)
|
|
{ if (cas==1)
|
|
{sort << "\n Algo_zero:_para:"
|
|
<< " prec_res_abs= " << prec_res_abs
|
|
<< " prec_res_rel= " << prec_res_rel
|
|
<< " iter_max= " << iter_max
|
|
<< " coef_mini_delta_x= " << coef_mini_delta_x
|
|
<< " mini_delta_x= " << mini_delta_x
|
|
<< " maxi_delta_x= " << maxi_delta_x
|
|
<< " nbMaxiIncre= " << nbMaxiIncre
|
|
<< " maxi_test_moins1= "<< maxi_test_moins1
|
|
<< " permet_affichage= " << permet_affichage
|
|
<< endl;
|
|
};
|
|
|
|
};
|
|
|
|
|
|
// recherche du zéro d'une équation du second degré dans l'espace des réels: ax^2+bx+c=0.
|
|
// cas indique les différents cas:
|
|
// = 1 : il y a deux racines, racine2 > racine1
|
|
// = 2 : il y a une racine double: racine1=racine2
|
|
// = 0 : pas de racine
|
|
// = -1 : pas de racine car a,b,c sont tous inférieur à la précision de traitement
|
|
// = -2 : pas de racine car a,b sont inférieur à la précision de traitement tandis que c est non nul
|
|
// = -3 : une racine simple car a est inférieur à la précision de traitement, donc considéré comme nul
|
|
// b et c sont non nul -> equa du premier degré, solution: racine1
|
|
void Algo_zero::SecondDegre(double a, double b, double c1, double& racine1, double& racine2, int& cas) const
|
|
{ // 1== on regarde tout d'abord si ce n'est pas une équation du premier degré
|
|
if (Dabs(a) <= ConstMath::trespetit)
|
|
{ // on regarde si le coeff restant de x n'est pas nul
|
|
if ( Dabs(b) <= ConstMath::trespetit)
|
|
{ // on regarde si le terme constant est nul
|
|
if (Dabs(c1) <= ConstMath::trespetit)
|
|
// pas de solution possible mais l'équation est valide à ConstMath::trespetit près
|
|
{ cas = -1;}
|
|
else
|
|
// pas de solution possible et l'équation n'est pas possible
|
|
{ cas = -2;}
|
|
}
|
|
else
|
|
// b est non nul et l'équation est du premier degré
|
|
{ cas = -3; racine1= -c1/b; racine2=racine1; return;};
|
|
};
|
|
// 2== maintenant on est sur d'avoir une équa du second degré
|
|
// calcul du discriminant
|
|
double delta = b*b - 4. * a * c1;
|
|
// choix en fonction de la valeur du discriminant
|
|
if (delta < 0.)
|
|
{ // on regarde si l'on est près de la précision de la machine
|
|
if (Dabs(delta/a) <= 1000*ConstMath::trespetit)
|
|
{ // dans ce cas on met delta à 0
|
|
cas = 2;
|
|
racine1 = (-b)/(2.*a); racine2=racine1;
|
|
}
|
|
else
|
|
// pas de solution
|
|
{cas = 0;}
|
|
}
|
|
else if (delta == 0.)
|
|
{ // une solution double
|
|
cas = 2;
|
|
racine1 = (-b)/(2.*a); racine2=racine1;
|
|
}
|
|
else
|
|
{ // deux racines séparées
|
|
cas = 1; double r = sqrt(delta);
|
|
racine1 = (-b-r)/(2.*a);racine2 = (-b+r)/(2.*a);
|
|
};
|
|
};
|
|
|
|
// recherche du zéro d'une équation du troisième degré dans l'espace des réels: ax^3+bx^2+cx+d=0.
|
|
// cas indique les différents cas:
|
|
// = 1 : il y a deux racines , racine2 > racine1
|
|
// = 2 : il y a une racine double
|
|
// = 3 : il y a une seule racine réelle racine1 (et deux complexes non calculées)
|
|
// = 4 : il y a une racine simple: racine1, et une racine double: racine 2
|
|
// = 5 : il y a 3 racines réelles: racine1, racine2, racine3
|
|
// = 0 : pas de racine
|
|
// = -1 : pas de racine car a,b,c sont tous inférieur à la précision de traitement
|
|
// = -2 : pas de racine car a,b sont inférieur à la précision de traitement tandis que c est non nul
|
|
// = -3 : une racine simple car a est inférieur à la précision de traitement, donc considéré comme nul
|
|
// b et c sont non nul -> equa du premier degré, solution: racine1
|
|
void Algo_zero::TroisiemeDegre(double a, double b, double c1, double d,
|
|
double& racine1, double& racine2, double& racine3, int& cas)
|
|
{
|
|
// 1== on regarde tout d'abord si ce n'est pas une équation du second degré
|
|
if (Dabs(a) <= ConstMath::trespetit)
|
|
{ // appel de la résolution d'une équation du second degré
|
|
SecondDegre(b,c1,d,racine1,racine2,cas);
|
|
#ifdef MISE_AU_POINT
|
|
if (cas==1)
|
|
{if ((Dabs(racine1) > ConstMath::pasmalgrand) || (Dabs(racine2) > ConstMath::pasmalgrand) )
|
|
{ cout << "\n sans doute erreur dans le calcul des racines de l'equation du 3ieme degre"
|
|
<< "\n resultat trop grand " << racine1 << " " << racine2;
|
|
}
|
|
}
|
|
else
|
|
{if (Dabs(racine1) > ConstMath::pasmalgrand)
|
|
{ cout << "\n sans doute erreur dans le calcul des racines de l'equation du 3ieme degre"
|
|
<< "\n resultat trop grand " << racine1;
|
|
}
|
|
};
|
|
#endif
|
|
}
|
|
else
|
|
{ // cas d'une réelle équation du 3ième degré
|
|
// on normalise l'équation par a;
|
|
double bb = b/a; double cc1 = c1/a;double dd = d/a;
|
|
// on appelle la routine de quartic
|
|
double v3[4];// les racines
|
|
int nb_root= algo_quartic.cubic(bb,cc1,dd,v3);
|
|
// on recopie les racines
|
|
cas = nb_root+2;
|
|
switch (nb_root)
|
|
{ case 3: racine3=v3[2];
|
|
/* #ifdef MISE_AU_POINT
|
|
double x=racine3;
|
|
if (Dabs(x*x*x + bb*x*x +cc1*x+dd) > ConstMath::petit)
|
|
{ cout << "\n erreur la troisieme racine de l'equation du 3ieme degre est fausse equa= "
|
|
<< Dabs(x*x*x + bb*x*x +cc1*x+dd)
|
|
<< "\n Algo_zero::TroisiemeDegre(..";
|
|
}
|
|
#endif*/
|
|
case 2: racine2=v3[1];
|
|
/* #ifdef MISE_AU_POINT
|
|
x=racine2;
|
|
if (Dabs(x*x*x + bb*x*x +cc1*x+dd) > ConstMath::petit)
|
|
{ cout << "\n erreur la seconde racine de l'equation du 3ieme degre est fausse equa= "
|
|
<< Dabs(x*x*x + bb*x*x +cc1*x+dd)
|
|
<< "\n Algo_zero::TroisiemeDegre(..";
|
|
}
|
|
#endif*/
|
|
case 1: racine1=v3[0];
|
|
/* #ifdef MISE_AU_POINT
|
|
x=racine1;
|
|
if (Dabs(x*x*x + bb*x*x +cc1*x+dd) > ConstMath::petit)
|
|
{ cout << "\n erreur la premiere racine de l'equation du 3ieme degre est fausse equa= "
|
|
<< Dabs(x*x*x + bb*x*x +cc1*x+dd)
|
|
<< "\n Algo_zero::TroisiemeDegre(..";
|
|
}
|
|
#endif*/
|
|
break;
|
|
default:
|
|
cout << "\n erreur pas de racine pour l'equation cubique !"
|
|
<< "\n Algo_zero::TroisiemeDegre(..";
|
|
Sortie(1);
|
|
}
|
|
};
|
|
|
|
// cout << "\n résolution 3 degré cas= " << cas << " racine1= " << racine1
|
|
// << " racine2= " << racine2 << " racine3= " << racine3;
|
|
|
|
};
|
|
|
|
|
|
/*
|
|
void Algo_zero::TroisiemeDegre(double a, double b, double c1, double d,
|
|
double& racine1, double& racine2, double& racine3, int& cas) const
|
|
{
|
|
|
|
// 1== on regarde tout d'abord si ce n'est pas une équation du second degré
|
|
if (Dabs(a) <= ConstMath::trespetit)
|
|
{ // appel de la résolution d'une équation du second degré
|
|
SecondDegre(b,c1,d,racine1,racine2,cas);
|
|
#ifdef MISE_AU_POINT
|
|
if (cas==1)
|
|
{if ((Dabs(racine1) > ConstMath::pasmalgrand) || (Dabs(racine2) > ConstMath::pasmalgrand) )
|
|
{ cout << "\n sans doute erreur dans le calcul des racines de l'equation du 3ieme degre"
|
|
<< "\n resultat trop grand " << racine1 << " " << racine2;
|
|
}
|
|
}
|
|
else
|
|
{if (Dabs(racine1) > ConstMath::pasmalgrand)
|
|
{ cout << "\n sans doute erreur dans le calcul des racines de l'equation du 3ieme degre"
|
|
<< "\n resultat trop grand " << racine1;
|
|
}
|
|
};
|
|
#endif
|
|
}
|
|
else
|
|
{ // cas d'une réelle équation du 3ième degré
|
|
// on normalise l'équation par a;
|
|
a=1.; b/=a; c1/=a;d/=a;
|
|
// on commence par regarder si ce n'est pas une racine triple
|
|
// si p est nulle il n'y a qu'une seule racine
|
|
|
|
|
|
|
|
|
|
|
|
double lambda = -b/3.;
|
|
if (( Dabs(c1- 3.*lambda*lambda) <= 100.*ConstMath::trespetit) &&
|
|
( Dabs(d + lambda*lambda*lambda) <= 100.*ConstMath::trespetit))
|
|
{ // cas d'une racine tripe
|
|
cas = 3; racine1=racine2=racine3=lambda;
|
|
#ifdef MISE_AU_POINT
|
|
double x=racine1;
|
|
if (Dabs(x*x*x + b*x*x +c1*x+d) > ConstMath::petit)
|
|
{ cout << "\n erreur la racine triple de l'equation du 3ieme degre est fausse equa= "
|
|
<< Dabs(x*x*x + b*x*x +c1*x+d);
|
|
}
|
|
#endif
|
|
return;
|
|
}
|
|
// cas où on n'a pas de racine triple
|
|
double r1=-b/(3.*a); //double x= r1 + y
|
|
double b2=b*b; double a2=a*a;double b3=b2*b; double a3=a*a2;
|
|
double p = (3*a*c1-b2)/(3.*a2); double q = (2.*b3+27.*a2*d-9.*a*b*c1)/(27.*a3);
|
|
double p3 = p*p*p; double q2=q*q;
|
|
double s=4.*p3 + 27.*q2;
|
|
if (s>0.)
|
|
{ // cas d'une seule racine réelle: racine1
|
|
double dd=sqrt(s/27.);
|
|
double untier=1./3.;
|
|
// résolution des pb de petits nombres ou de nombres négatifs
|
|
double ddmoinsq=(-q+dd);
|
|
double signeddmoinsq=1.; if (ddmoinsq < 0.) signeddmoinsq=-1.;
|
|
ddmoinsq = Dabs(ddmoinsq);
|
|
double qPlusDd =(q+dd);
|
|
double signeqPlusDd=1.; if (qPlusDd < 0.) signeqPlusDd=-1.;
|
|
qPlusDd = Dabs(qPlusDd);
|
|
// si le nombre est trop petit le log tend vers -l'infini
|
|
if (ddmoinsq <= 100.*ConstMath::trespetit) {ddmoinsq=1.;signeddmoinsq=0.;};
|
|
if (qPlusDd <= 100.*ConstMath::trespetit) {qPlusDd=1.;signeqPlusDd=0.;};
|
|
double y = signeddmoinsq * pow(ddmoinsq*0.5,untier)
|
|
-signeqPlusDd*pow(qPlusDd*0.5,untier);
|
|
// tout ça à la place de l'expression qui suit !!
|
|
// double y = pow((-q+dd)*0.5,untier)+pow(-(q+dd)*0.5,untier);
|
|
racine1=racine2=racine3=r1+y;
|
|
#ifdef MISE_AU_POINT
|
|
if (Dabs(racine1) > ConstMath::pasmalgrand)
|
|
{ cout << "\n sans doute erreur dans le calcul des racines de l'equation du 3ieme degre"
|
|
<< "\n resultat trop grand " << racine1;
|
|
}
|
|
#endif
|
|
#ifdef MISE_AU_POINT
|
|
double x=racine1;
|
|
if (Dabs(x*x*x + b*x*x +c1*x+d) > ConstMath::petit)
|
|
{ cout << "\n erreur la racine simple de l'equation du 3ieme degre est fausse equa= "
|
|
<< Dabs(x*x*x + b*x*x +c1*x+d);
|
|
}
|
|
#endif
|
|
cas= 3;
|
|
}
|
|
else if (s == 0.)
|
|
{ // cas où il y a une racine simple: racine 1, et une racine double: racine2
|
|
double y1= 3.*q/p; double y2=-y1/2.;
|
|
racine1=r1+y1; racine2=r1+y2; cas = 4;
|
|
#ifdef MISE_AU_POINT
|
|
if ((Dabs(racine1) > ConstMath::pasmalgrand) || (Dabs(racine2) > ConstMath::pasmalgrand) )
|
|
{ cout << "\n sans doute erreur dans le calcul des racines de l'equation du 3ieme degre"
|
|
<< "\n resultat trop grand " << racine1 << " " << racine2;
|
|
}
|
|
#endif
|
|
#ifdef MISE_AU_POINT
|
|
double x=racine1;
|
|
if (Dabs(x*x*x + b*x*x +c1*x+d) > ConstMath::petit)
|
|
{ cout << "\n erreur la racine simple premiere de l'equation du 3ieme degre est fausse equa= "
|
|
<< Dabs(x*x*x + b*x*x +c1*x+d);
|
|
}
|
|
x=racine2;
|
|
if (Dabs(x*x*x + b*x*x +c1*x+d) > ConstMath::petit)
|
|
{ cout << "\n erreur la racine double de l'equation du 3ieme degre est fausse equa= "
|
|
<< Dabs(x*x*x + b*x*x +c1*x+d);
|
|
}
|
|
#endif
|
|
}
|
|
else // cas où s < 0.
|
|
{ // cas où il y a trois racines réelles
|
|
double dd=sqrt((-4.*p3+27.*q2)/27.); double r=2.*sqrt(-p/3.);
|
|
double teta1=ConstMath::Pi/2.;
|
|
if (Dabs(q) > 10*ConstMath::trespetit)teta1 = atan(-dd/q);
|
|
if (cos(teta1)*(-q) < 0) teta1 +=ConstMath::Pi/2;
|
|
// if (q<0.) {teta1 = atan(-dd/q);}
|
|
// else if (q>0.) {teta1 = atan(-dd/q) + ConstMath::Pi;};
|
|
// // sinon teta1=pi/2 ce qui est déjà le cas
|
|
double y1=r*cos(teta1/3.);
|
|
double y2=r*cos(teta1/3.+2.*ConstMath::Pi/3.);
|
|
double y3=r*cos(teta1/3.+4.*ConstMath::Pi/3.);
|
|
racine1=r1+y1; racine2=r1+y2; racine3=r1+y3;
|
|
cas = 5;
|
|
#ifdef MISE_AU_POINT
|
|
if ((Dabs(racine1) > ConstMath::pasmalgrand) || (Dabs(racine2) > ConstMath::pasmalgrand)
|
|
|| (racine3 > ConstMath::pasmalgrand))
|
|
{ cout << "\n sans doute erreur dans le calcul des racines de l'equation du 3ieme degre"
|
|
<< "\n resultat trop grand " << racine1 << " " << racine2 << " " << racine2;
|
|
}
|
|
#endif
|
|
#ifdef MISE_AU_POINT
|
|
double x=racine1;
|
|
if (Dabs(x*x*x + b*x*x +c1*x+d) > ConstMath::petit)
|
|
{ cout << "\n erreur la premiere racine (des 3) de l'equation du 3ieme degre est fausse equa= "
|
|
<< Dabs(x*x*x + b*x*x +c1*x+d);
|
|
}
|
|
x=racine2;
|
|
if (Dabs(x*x*x + b*x*x +c1*x+d) > ConstMath::petit)
|
|
{ cout << "\n erreur la seconde racine (des 3) de l'equation du 3ieme degre est fausse equa= "
|
|
<< Dabs(x*x*x + b*x*x +c1*x+d);
|
|
}
|
|
x=racine3;
|
|
if (Dabs(x*x*x + b*x*x +c1*x+d) > ConstMath::petit)
|
|
{ cout << "\n erreur la troisieme racine (des 3) de l'equation du 3ieme degre est fausse equa= "
|
|
<< Dabs(x*x*x + b*x*x +c1*x+d);
|
|
}
|
|
#endif
|
|
};
|
|
};
|
|
}
|
|
*/
|