Commit 1f23f102 authored by Sebastien Timmermans's avatar Sebastien Timmermans 🎹
Browse files

[New] error of NR loop closing close geo in mbs_part

parent 847d1c75
......@@ -44,7 +44,7 @@ MbsAux * initMbsAux(MbsData *s)
Ncons = s->Ncons;
mbs_aux->NRerr = s->NRerr;
mbs_aux->MAX_NR_ITER = 100;
mbs_aux->MAX_NR_ITER = 97;
if (Ncons>0)
{
......
......@@ -26,19 +26,16 @@ int mbs_close_geo(MbsData *s, MbsAux *mbs_aux)
int i,j;
int iter=0;
int nL,nC;
int err = 0;
double d;
iter = 0;
mbs_aux->norm_h=1.0;
while((mbs_aux->norm_h > mbs_aux->NRerr) && (iter++ <= mbs_aux->MAX_NR_ITER))
{
// Calcul des contraintes et de la Jacobienne
mbs_calc_hJ(s,mbs_aux);
// Norme des contraintes (en supposant que toutes les contraintes independantes sont au debut ???)
mbs_aux->norm_h = norminf_vector(mbs_aux->h,s->nhu);
......@@ -54,8 +51,14 @@ int mbs_close_geo(MbsData *s, MbsAux *mbs_aux)
}
// Decomposition LU de la matrice -Jv
ludcmp(mbs_aux->mJv,s->nqv,mbs_aux->ind_mJv,&d);
err = ludcmp(mbs_aux->mJv,s->nqv,mbs_aux->ind_mJv,&d);
if (err < 0)
{
mbs_msg("\t >CLOSE GEO> Error in mbs_close_geo du LU decomposition of Jv matrix \n");
return err;
}
if(mbs_aux->norm_h > mbs_aux->NRerr)
{
// err
......
......@@ -610,8 +610,6 @@ void mbs_run_part(MbsPart* mbs_part, MbsData* mbs_data)
err = PART_coord_part(mbs_data, mbs_part, ind_u_des, mbs_data->nqu, ind_c, mbs_data->nqc, 1); // d'abord on cherche un choix valable
if (err < 0)
{
mbs_msg(">>PART>>\n");
......@@ -636,17 +634,19 @@ void mbs_run_part(MbsPart* mbs_part, MbsData* mbs_data)
{
mbs_msg(">>PART>>\n");
mbs_msg(">>PART>> ***** mbs_run_part : impossible to close the MBS : *****\n");
mbs_msg(">>PART>> after %d iteration of Newton Raphson \n", nbiter);
mbs_msg(">>PART>> -> try with other desired independent variables u\n");
mbs_msg(">>PART>> and/or\n");
mbs_msg(">>PART>> -> try with other initial values for q (u, v)\n");
mbs_msg(">>PART>> -> ... but have a look at the proposed {u,v} partition\n");
mbs_msg(">>PART>> Loop closing Error: Newton-Raphson iteration overrun during partitioning !\n");
mbs_msg(">>PART>>\n");
mbs_error_msg("Loop closing Error: Newton-Raphson iteration overrun during partitioning !\n");
}
mbs_error_msg(" [%d] in mbs_run_part !! \n", 100 + nbiter); }
else if (nbiter < 0)
{
mbs_msg(">>PART>>\n");
mbs_msg(">>PART>> Error during mbs_part \n");
mbs_error_msg(" [%d] in mbs_run_part !! \n", -100 + nbiter);
}
else
{
......@@ -659,7 +659,7 @@ void mbs_run_part(MbsPart* mbs_part, MbsData* mbs_data)
// Second partitioning on the closed structure (possible refinment still on the basis of ind_u_des)
err = PART_coord_part(mbs_data, mbs_part, ind_u_des, mbs_data->nqu, ind_c, mbs_data->nqc, 0);
if (err == -1) // not very likely but should be considered ...
if (err < 0) // not very likely but should be considered ...
{
mbs_msg(">>PART>>\n");
mbs_msg(">>PART>> ***** mbs_run_part : Second coordinate partitioning interrupted : *****\n");
......
......@@ -15,7 +15,7 @@
#define NRANSI
#define TINY 1.0e-20
void ludcmp(double **a, int n, int *indx, double *d)
int ludcmp(double **a, int n, int *indx, double *d)
{
int i,imax,j,k;
double big,dum,sum,temp;
......@@ -25,12 +25,15 @@ void ludcmp(double **a, int n, int *indx, double *d)
vv=get_dvec_1(n);
*d=1.0;
for (i=1;i<=n;i++) {
big=0.0;
for (j=1;j<=n;j++)
if ((temp=fabs(a[i][j])) > big) big=temp;
for (i = 1; i <= n; i++) {
big = 0.0;
for (j = 1; j <= n; j++)
if ((temp = fabs(a[i][j])) > big) big = temp;
if (big == 0.0)
mbs_error_msg("Singular matrix in routine ludcmp");
{
mbs_msg("\t >LU DCMP> Singular matrix in routine ludcmp");
return -2;
}
vv[i]=1.0/big;
}
for (j=1;j<=n;j++) {
......@@ -69,7 +72,7 @@ void ludcmp(double **a, int n, int *indx, double *d)
free_dvec_1(vv);
}
void ludcmp_0(double **a, int n, int *indx, double *d)
int ludcmp_0(double **a, int n, int *indx, double *d)
{
int i, imax, j, k;
double big, dum, sum, temp;
......@@ -85,7 +88,8 @@ void ludcmp_0(double **a, int n, int *indx, double *d)
if ((temp = fabs(a[i][j])) > big) big = temp;
if (big == 0.0)
{
mbs_error_msg("Singular matrix in routine ludcmp\n");
mbs_msg("\t >LU DCMP> Singular matrix in routine ludcmp_0");
return -2;
}
vv[i] = 1.0 / big;
}
......
......@@ -11,8 +11,8 @@
#define nr_fct_h
/*--------------------*/
void ludcmp(double **a, int n, int *indx, double *d);
void ludcmp_0(double **a, int n, int *indx, double *d);
int ludcmp(double **a, int n, int *indx, double *d);
int ludcmp_0(double **a, int n, int *indx, double *d);
void lubksb(double **a, int n, int *indx, double b[]);
void lubksb_0(double **a, int n, int *indx, double b[]);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment