Commit 292586cf authored by Olivier Lantsoght's avatar Olivier Lantsoght 🏁

Merge branch 'release_1.11.2' into 'dev'

Release 1.11.2

See merge request !246
parents ab296c4b 2f06848b
Pipeline #5024 passed with stages
in 19 minutes and 41 seconds
## Changes for next release
* [New] error management in mbs load, part and dirdyn (not yet in equil and modal)
* ? Coupling VPM branch merge : to check with DG
## Robotran v1.11.2
* [New] error management in all module
* See the wiki in gitlab for the corresponding .pptx
* The results files are written when an error is raised during dirdyn process.
* [New] CI automatic tests added for equil module + file tree of example project remodeled.
* The results files are written when an error is raised during dirdyn and solvekin process.
* [New][CI]
* automatic tests added for equil module + file tree of example project remodeled.
* Ubuntu 18.04 is used for testing
* [New] usefull function`scalar_product_0` and `matrix_product_0` created (previously exists for fixed size vector)
* [New] module called [Solvekin] to study kinematics. Refer to MBsysC documentation for instruction or example in `/ExampleProjects/TutorialProjects/analysisModules/`
* [Python] interface:
* [New] mbs_gensensor is callable from Python interface
* [Update] interface for W10 working, see instructions in the "readme_aux/interface_python.md"
* [Fix] The user model include the "lut.h" header file when needed.
Not yet done :
* ? Coupling VPM branch merge : to check with DG
* ? error management in equil and modal
* [DEV] The arguments of some core function of Modal and Equilibrium module has changed. No impact for user.
## Robotran v1.11.1
* Fixing the server call for symbolic generation
......
......@@ -107,7 +107,7 @@ pre.writeBodies(bdys,chemin=datbox_path)
pre.writeDrvDof(bdys,chemin=datbox_path)
pre.writeDofIni(bdys,chemin=datbox_path)
pre.writeModels(mods,chemin=datbox_path)
pre.writeBulkBehav(mats,chemin=datbox_path)
pre.writeBulkBehav(mats,chemin=datbox_path,dim=dim, gravy=[0.,-9.81,0.])
pre.writeTactBehav(tacs,sees,chemin=datbox_path)
pre.writeVlocRlocIni(chemin=datbox_path)
pre.writePostpro(commands=post, parts=bdys, path=datbox_path)
......
......@@ -133,6 +133,7 @@ int main(int argc, char const *argv[])
mbs_data->q[R2_pendulum_id] = 0.2;
mbs_equil = mbs_new_equil(mbs_data);
// equil options (see documentations for additional options)
mbs_equil->options->grad_lpk = 1; // using another method of computation
mbs_equil->options->senstol = 1e-06;
mbs_equil->options->verbose = 1;
mbs_equil->options->save2file = 1;
......@@ -181,6 +182,7 @@ int main(int argc, char const *argv[])
mbs_equil->options->verbose = 0;
// --- Variable addition, xe (see user_equil_fxe)
mbs_equil->options->nxe = 1;
mbs_equil->options->grad_lpk = 1; // using another method of computation
mbs_equil->options->compute_uxd = 0;
mbs_equil_addition(mbs_equil->options);
mbs_equil->options->xe_ptr[1] = &(mbs_data->Qa[R2_pendulum_id]);
......@@ -202,6 +204,7 @@ int main(int argc, char const *argv[])
// equil options (see documentations for additional options)
mbs_equil->options->mode = 2;
mbs_equil->options->senstol = 1e-06;
mbs_equil->options->grad_lpk = 1; // using another method of computation
mbs_equil->options->verbose = 1;
mbs_equil->options->compute_uxd = 0;
// set a desired speed for the cart (i.e. the speed of the all system)
......@@ -246,6 +249,7 @@ int main(int argc, char const *argv[])
mbs_equil = mbs_new_equil(mbs_data);
// equil options (see documentations for additional options)
mbs_equil->options->senstol = 1e-06;
mbs_equil->options->grad_lpk = 1; // using another method of computation
mbs_equil->options->verbose = 1;
mbs_equil->options->save2file = 1;
mbs_equil->options->compute_uxd = 0;
......
......@@ -551,6 +551,10 @@ MbsData* MDS_create_MBSdataStruct(MDS_gen_strct* mds_gen_strct, MbsData* s)
s->lambda[i] = 0.0;
}
}
else
{
s->lambda = NULL;
}
// Link Data
......@@ -743,7 +747,7 @@ void mbs_delete_data(MbsData *s)
if (s->lrod != NULL)
free(s->lrod);
if (s->nqv)
if (s->lambda != NULL)
free(s->lambda);
// Links
......
......@@ -51,7 +51,7 @@ MbsEquil* mbs_new_equil_aux(MbsData* mbs_data, MbsAux* mbs_aux)
equil->nx = 0;
equil->x = NULL;
equil->nxs;
equil->nxs = 0;
equil->xs = NULL;
equil->nxns = 0;
equil->xns = NULL;
......@@ -74,9 +74,7 @@ MbsEquil* mbs_new_equil_aux(MbsData* mbs_data, MbsAux* mbs_aux)
opts->relax = 0;
opts->relaxcoeff = 0.9;
opts->soft = 0;
opts->senstol = 1e-06;
opts->equitol = 1e-06;
opts->devjac = 1e-06;
opts->itermax = 30;
opts->mode = 1; // static equilibrium
opts->verbose = 1;
......@@ -104,6 +102,17 @@ MbsEquil* mbs_new_equil_aux(MbsData* mbs_data, MbsAux* mbs_aux)
opts->buffersize = -1;
opts->max_save_user = 12;
// gradient computation through dev
opts->senstol = 1e-06;
opts->devjac = 1e-06;
// gradient computation through lpk
opts->grad_lpk = 0;
opts->lpk_itermax = 10;
opts->lpk_relincr = 1e-2;
opts->lpk_absincr = 1e-3;
opts->lpk_equitol = 1e-6;
opts->lpk_lintol = 1e-3;
equil->options = opts;
return equil;
......@@ -126,8 +135,10 @@ void mbs_delete_equil(MbsEquil* eq, MbsData* mbs_data)
if (eq->options->compute_uxd==1 && mbs_data->Nux > 0) { free_dmat_0(eq->grad_uxd); }
if (eq->nxe > 0) { free_dmat_0(eq->grad_fxe); }
free_dmat_0(eq->grad);
mbs_delete_lpk(eq->lpk);
if (eq->options->grad_lpk == 1)
{
mbs_delete_lpk(eq->lpk);
}
free(eq->qd_saved);
free(eq->qdd_saved);
//-------
......@@ -364,7 +375,17 @@ void mbs_equil_init(MbsEquil* eq, MbsData* s)
if (s->Nux > 0) { eq->grad_uxd = get_dmat_0(s->Nux, eq->nx); }
if (eq->nxe > 0) { eq->grad_fxe = get_dmat_0(eq->nxe, eq->nx); }
eq->grad = get_dmat_0(eq->nx, eq->nx);
eq->lpk = mbs_new_lpk(s, eq->nx, MAX(MAX(s->nqu, s->Nux),eq->nxe));
if (eq->options->grad_lpk == 1)
{
eq->lpk = mbs_new_lpk(s, eq->nx, MAX(MAX(s->nqu, s->Nux), eq->nxe));
eq->lpk->itermax = eq->options->lpk_itermax;
eq->lpk->relincr = eq->options->lpk_relincr;
eq->lpk->absincr = eq->options->lpk_absincr;
eq->lpk->equitol = eq->options->lpk_equitol;
eq->lpk->lintol = eq->options->lpk_lintol;
}
// 6. Equilibrium saving to file (anim, results)
if (eq->options->save2file)
......@@ -711,7 +732,14 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
// 1. Gradient test - search for sensitive variables
// 1.1 Gradient estimation for f(x) :
err = mbs_equil_grad(eq, s);
if (eq->options->grad_lpk == 1)
{
err = mbs_equil_grad_lpk(eq, s);
}
else
{
err = mbs_equil_grad_dev(eq, s);
}
if (err<0) // Error management
{
mbs_msg(">>EQUIL>>\n");
......@@ -828,16 +856,20 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
}
else // numerical estimation (similar to 1.1)
{
eq->lpk->verbose = 0;
eq->options->verbose = eq->options->verbose - 1;
err=mbs_equil_grad(eq, s);
if (eq->options->grad_lpk == 1)
{
err = mbs_equil_grad_lpk(eq, s);
}
else
{
err = mbs_equil_grad_dev(eq, s);
}
if (err<0) // Error management
{
mbs_msg(">>EQUIL>>\n");
mbs_msg(">>EQUIL>> Error during equilibrium process at iteration %d !\n", eq->iter);
mbs_error_msg("[%d] in mbs_equil_fsolvepk !! \n", -200 + err);
}
eq->options->verbose = eq->options->verbose + 1;
}
// Gradient Resizing (sensitive variables only) : [nfs * nx]
for (i = 0; i < eq->nxs; i++)
......@@ -1040,12 +1072,16 @@ void mbs_equil_fsolvepk(int(*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*,
free_dvec_0(xsR);
}
int mbs_equil_grad(MbsEquil *eq,MbsData *s)
int mbs_equil_grad_lpk(MbsEquil *eq,MbsData *s)
{
int err = 0;
int i,j;
int EqChoice;
for (i = 0; i < eq->lpk->nx; i++) { eq->lpk->x_ptr[i] = eq->x_ptr[i]; }
for (i = 0; i < eq->lpk->nx; i++)
{
eq->lpk->x_ptr[i] = eq->x_ptr[i];
}
if (0) // disactivated for visibility...
{
......@@ -1111,6 +1147,99 @@ int mbs_equil_grad(MbsEquil *eq,MbsData *s)
return err;
}
int mbs_equil_grad_dev(MbsEquil *eq, MbsData *s)
{
int i, j;
int err = 0;
int isSensitive;
int nx, nf;
double *xd, *f, *fd; // xd is used for different calculations.
int *xs, *xns;
int nxs, nxns;
// 0. initialization (and allocation)
nx = eq->nx;
nf = nx;
xd = get_dvec_0(nx);
f = get_dvec_0(nx);
fd = get_dvec_0(nx);
xs = get_ivec_0(nx);
xns = get_ivec_0(nx);
// 0.1 Solver might either use a analythical evaluation or a numerical estimation for the Jacobian
// in case of analythical evaluation a additional test might be necessary
// 1. Gradient test - search for sensitive variables
err = mbs_equil_fct(eq->x, nx, f, eq, eq->aux, s);
if (err < 0) // Error management
{
mbs_msg("\t in EQUIL GRAD, computation of grad_Rr = delta Rr/delta x \n");
return err;
}
copy_dvec_0(eq->x, xd, nx);
for (i = 0; i < nx; i++) //xd is, here, a small deviation from x
{
if (fabs(eq->x[i]) > 1e-2) // if x[i] !=0
{
xd[i] = eq->x[i] * (1.0 + eq->options->devjac);
}
else // if x[i] ~=0
{
xd[i] = eq->x[i] + eq->options->devjac;
}
err = mbs_equil_fct(xd, nx, fd, eq, eq->aux, s);
if (err < 0) // Error management
{
mbs_msg("\t in EQUIL GRAD, computation of grad_Rr = delta Rr/delta x \n");
return err;
}
for (j = 0; j < nf; j++) // gradient estimation
{
eq->grad[j][i] = (fd[j] - f[j]) / (xd[i] - eq->x[i]);
}
xd[i] = eq->x[i];
}
// 1.2. Gradient testing (removal of non-sensitive variables)
nxs = 0;
nxns = 0;
for (i = 0; i < nx; i++) // test all equilibrium variables
{
isSensitive = 0;
for (j = 0; j < nf; j++)
{
if (fabs(eq->grad[j][i]) > fabs(eq->options->senstol)) // sensitive if there is a non-zero element among the ith column
{
isSensitive = 1;
break;
}
}
if (isSensitive) // the ith variable is sensitive
{
xs[nxs] = i; // write the vector, eq->xs, with the index of sensitive variables
(nxs)++; // increment the number of sensitive variables
}
else // the ith variable is non sensitive (all the column is 0)
{
xns[nxns] = i; // write the vector, eq->xs, with the index of non sensitive variables
(nxns)++; // increment the number of non sensitive variables
for (j = 0; j < nf; j++)
{
eq->grad[j][i] = 0.0; // set gradient line for the non sensitive variables to 0
}
}
}
free_dvec_0(xd);
free_dvec_0(f);
free_dvec_0(fd);
free_ivec_0(xs);
free_ivec_0(xns);
return err;
}
void mbs_print_equil(MbsEquil* eq)
......
......@@ -134,12 +134,19 @@ void mbs_print_equil(MbsEquil* eq);
void mbs_equil_fsolvepk(int (*fun_ptr)(double*, int, double*, MbsEquil*, MbsAux*, MbsData*), MbsEquil *eq, MbsAux *aux, MbsData *s );
/**
* Compute the gradient of f(x) for a given configuration x
* Compute the gradient of f(x) for a given configuration x based on the lpk parabolic fitting method
*
* \p eq ?Q to change for a more general use of function of fsolvepk
* \p s ?Q to change for a more general use of function of fsolvepk
*/
int mbs_equil_grad(MbsEquil *eq, MbsData *s);
int mbs_equil_grad_lpk(MbsEquil *eq, MbsData *s);
/**
* Compute the gradient of f(x) for a given configuration x based on a deviation computation
*
* \p eq ?Q to change for a more general use of function of fsolvepk
* \p s ?Q to change for a more general use of function of fsolvepk
*/
int mbs_equil_grad_dev(MbsEquil *eq, MbsData *s);
#endif
......@@ -233,7 +233,7 @@ int mbs_linearipk(double **GK, MbsLpk *lpk, MbsAux *aux, MbsData *s, int EqChoi
// waiting column k is still not acepted
else // Error management
{
mbs_msg(">>LINPK>> Column %d/%d no convergence after %d iterations\n", k + 1, lpk->nx, lpk->itermax);
mbs_msg(">>LINPK>> Column %d/%d no convergence after %d iterations\n", lpk->diverge_ind[k] + 1, lpk->nx, lpk->itermax);
mbs_msg(">>LINPK>> Restart with other parameters !\n");
err = -10;
return err;
......
......@@ -32,7 +32,7 @@ MbsModal* mbs_new_modal(MbsData* mbs_data)
return mbs_new_modal_aux(mbs_data, mbs_aux);
}
MbsModal* mbs_new_modal_aux(MbsData* mbs_data, MbsAux* mbs_aux)
MbsModal* mbs_new_modal_aux(MbsData* s, MbsAux* mbs_aux)
{
MbsModal *mo; // short for modal general structure : mbs_modal
MbsModalOptions *opts;
......@@ -60,8 +60,15 @@ MbsModal* mbs_new_modal_aux(MbsData* mbs_data, MbsAux* mbs_aux)
opts->compute_JS = 0; // as it is quite time-consuming
opts->norm_evec = 1;
mo->options = opts;
// gradient computation
opts->lpk_itermax = 10;
opts->lpk_relincr = 1e-2;
opts->lpk_absincr = 1e-3;
opts->lpk_equitol = 1e-6;
opts->lpk_lintol = 1e-3;
mo->options = opts;
return mo;
}
......@@ -95,6 +102,12 @@ void mbs_modal_init(MbsModal *mo, MbsData *s)
}
mo->lpk = mbs_new_lpk(s, MAX(s->nqu, s->Nux), MAX(s->nqu, s->Nux)); // voluntarily oversized (in order to avoid extra lpk structure)
mo->lpk->itermax=mo->options->lpk_itermax;
mo->lpk->relincr = mo->options->lpk_relincr;
mo->lpk->absincr = mo->options->lpk_absincr;
mo->lpk->equitol = mo->options->lpk_equitol;
mo->lpk->lintol = mo->options->lpk_lintol;
mo->ss = mbs_new_statespace(s->nqu, s->Nux, 0, 0, mo->options->compute_uxd, mo->options->compute_JS);
if (mo->options->verbose == 1)
......@@ -915,6 +928,8 @@ int mbs_modal_save_anim(MbsModal* mo, MbsData* s, char* filepath)
}
free_dvec_0(q_print);
free_dvec_0(q_mode_r);
free_dvec_0(q_mode_phi);
free(filename);
return err; // signature should be changed if no return is needed
......
......@@ -32,7 +32,7 @@ MbsModal* mbs_new_modal(MbsData* mbs_data);
* \p mbs_data the data structure of the model for which
* the modal analysis will be computed
*/
MbsModal* mbs_new_modal_aux(MbsData* mbs_data,MbsAux* mbs_aux);
MbsModal* mbs_new_modal_aux(MbsData* s, MbsAux* mbs_aux);
/**
*
......
......@@ -607,9 +607,9 @@ int mbs_fct_solvekin(MbsData *s, MbsSolvekin *sk)
err3 = mbs_lut_res_interp(sk->qdd, s->tsim, lutqdd);
for(i=1; i<=s->nqu; i++){
index = s->qu[i];
s->q[index] = lutq[index];
s->qd[index] = lutqd[index];
s->qdd[index] = lutqdd[index];
s->q[index] = lutq[i];
s->qd[index] = lutqd[i];
s->qdd[index] = lutqdd[i];
}
}
if (err1<0){ // Error management
......
......@@ -16,7 +16,7 @@
#include <stdio.h>
#include <string.h>
#include "glm/gtx/transform.hpp"
#include "gtx/transform.hpp"
#define SEPARATOR " ,\n\r\t"
......
......@@ -17,6 +17,13 @@ typedef struct MbsEquilOptions
int save2file; ///< 1: results saved 0: not saved [default = 1]
int compute_uxd; ///< no = 0, yes = 1, flag to compute the extra constitutive differential equations as equil equations [default = 1]
int grad_lpk; ///< options to compute the gradient through the lpk parabolic fitting method. If not activated, the gradient is computed by small deviation with options devJac [default = 0]
int lpk_itermax; ///< options for the lpk gradient [default = 10]
double lpk_relincr; ///< options for the lpk gradient [default = 1e-2]
double lpk_absincr; ///< options for the lpk gradient [default = 1e-3]
double lpk_equitol; ///< options for the lpk gradient [default = 1e-6]
double lpk_lintol; ///< options for the lpk gradient [default = 1e-3]
char* resfilename; ///< The keyword used for determining the name of result files [default:equil]
const char* respath; ///< Path in which result file are saved. [default:resultsR folder]
const char* animpath; ///< Path in which animation file are saved. [default:animationR folder]
......
......@@ -23,6 +23,14 @@ typedef struct MbsModalOptions
int compute_JS; ///< no = 0, yes = 1, flag to compute the non diagonal terms [default = 0] TO IMPLEMENT !
int norm_evec; ///< no = 0, yes = 1, flag to normalize the eigenvectors [default = 1]
// lpk linearization options
int lpk_itermax; ///< options for the lpk gradient [default = 10]
double lpk_relincr; ///< options for the lpk gradient [default = 1e-2]
double lpk_absincr; ///< options for the lpk gradient [default = 1e-3]
double lpk_equitol; ///< options for the lpk gradient [default = 1e-6]
double lpk_lintol; ///< options for the lpk gradient [default = 1e-3]
/* The keyword used for determining the name of result files */
char* resfilename;
......
......@@ -193,7 +193,7 @@ int mbs_lut_res_interp(MbsLutRes *lut, double t, double *coord)
// checking time validity
if (lut->tq[max][0] < t || lut->tq[min][0] > t){ // Error management
mbs_msg("\t >LOOK UP TABLE> error in mbs_lut_res_interp. Requested time is outside the domain of the provided file.\n");
mbs_msg("\t > Time is %f, min. and max. time allowed are %f and %f.\n",t, lut->tq[min, 0], lut->tq[max, 0]);
mbs_msg("\t > Time is %f, min. and max. time allowed are %f and %f.\n",t, lut->tq[min][0], lut->tq[max][0]);
return -2;
}
// checking easy solution
......
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