# File kernel/src/simulationTools/EulerMoreauOSI.hpp¶

Go to the source code of this file

class EulerMoreauOSI : public OneStepIntegrator
#include <>

One Step time Integrator for First Order Dynamical Systems.

This integrator is the work horse of the event&#8212;capturing time stepping schemes for first order systems. It is mainly based on some extensions of the Backward Euler and $$\theta-\gamma$$ schemes proposed in the pionnering work of J.J. Moreau for the sweeping process

J.J. Moreau. Evolution problem associated with a moving convex set in a Hilbert space. Journal of Differential Equations, 26, pp 347&#8212;374, 1977.

Variants are now used to integrate LCS, Relay systems, Higher order sweeping process see for instance

Consistency of a time-stepping method for a class of piecewise linear networks

M.K. Camlibel, W.P.M.H. Heemels, and J.M. Schumacher IEEE Transactions on Circuits and Systems I, 2002, 49(3):349&#8212;357

Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics

V Acary, B Brogliato Springer Verlag 2008

Convergence of time-stepping schemes for passive and extended linear complementarity systems L. Han, A. Tiwari, M.K. Camlibel, and J.-S. Pang SIAM Journal on Numerical Analysis 2009, 47(5):3768-3796

On preserving dissipativity properties of linear complementarity dynamical systems with the &theta-method

Greenhalgh Scott, Acary Vincent, Brogliato Bernard Numer. Math., , 2013.

Main time&#8212;integration schemes are based on the following $$\theta-\gamma$$ scheme

$\begin{split} \begin{cases} \label{eq:toto1} M x_{k+1} = M x_{k} +h\theta f(x_{k+1},t_{k+1})+h(1-\theta) f(x_k,t_k) + h \gamma r(t_{k+1}) + h(1-\gamma)r(t_k) \\[2mm] y_{k+1} = h(t_{k+1},x_{k+1},\lambda _{k+1}) \\[2mm] r_{k+1} = g(x_{k+1},\lambda_{k+1},t_{k+1})\\[2mm] \mbox{nslaw} ( y_{k+1} , \lambda_{k+1}) \end{cases} \end{split}$

where $$\theta = [0,1]$$ and $$\gamma \in [0,1]$$. As in Acary & Brogliato 2008, we call the previous problem the `one–step nonsmooth problem’.

Another variant can also be used (FullThetaGamma scheme)

$\begin{split} \begin{cases} M x_{k+1} = M x_{k} +h f(x_{k+\theta},t_{k+1}) + h r(t_{k+\gamma})\\[2mm] y_{k+\gamma} = h(t_{k+\gamma},x_{k+\gamma},\lambda _{k+\gamma})\\[2mm] r_{k+\gamma} = g(x_{k+\gamma},\lambda_{k+\gamma},t_{k+\gamma})\\[2mm] \mbox{nslaw} ( y_{k+\gamma} , \lambda_{k+\gamma}) \end{cases} \end{split}$

EulerMoreauOSI class is used to define some time-integrators methods for a list of first order dynamical systems. A EulerMoreauOSI instance is defined by the value of theta and possibly gamma and the list of concerned dynamical systems.

Each DynamicalSystem is associated to a SiconosMatrix, named “W”, which is the “iteration” matrix. W matrices are initialized and computed in initializeIterationMatrixW and computeW. Depending on the DS type, they may depend on time t and DS state x.

For first order systems, the implementation uses _r for storing the the input due to the nonsmooth law. This EulerMoreauOSI scheme assumes that the relative degree is zero or one and one level for _r is sufficient

Main functions:

• computeFreeState(): computes xfree (or vfree), dynamical systems state without taking non-smooth part into account

• updateState(): computes x (q,v), the complete dynamical systems states.

See User’s guide, for details.

Public Types

enum EulerMoreauOSI_ds_workVector_id

Values:

enumerator RESIDU
enumerator RESIDU_FREE
enumerator FREE
enumerator X_PARTIAL_NS_FOR_RELATION
enumerator DELTA_X_FOR_RELATION
enumerator LOCAL_BUFFER
enumerator WORK_LENGTH
enum EulerMoreauOSI_interaction_workVector_id

Values:

enumerator OSNSP_RHS
enumerator VEC_X
enumerator H_ALPHA
enumerator VEC_RESIDU_Y
enumerator VEC_RESIDU_R
enumerator YOLD
enumerator LAMBDAOLD
enumerator WORK_INTERACTION_LENGTH
enum EulerMoreauOSI_interaction_workBlockVector_id

Values:

enumerator XFREE
enumerator X_PARTIAL_NS
enumerator DELTA_X
enumerator G_ALPHA
enumerator BLOCK_WORK_LENGTH
enum EulerMoreauOSI_interaction_workMat_id

Values:

enumerator MAT_KHAT
enumerator MAT_KTILDE
enumerator MAT_WORK_LENGTH

Public Functions

EulerMoreauOSI(double theta)

constructor from theta value only

Parameters

theta – value for all DS.

EulerMoreauOSI(double theta, double gamma)

constructor from theta value only

Parameters
• theta – value for all linked DS.

• gamma – value for all linked DS.

inline virtual ~EulerMoreauOSI()

destructor

const SimpleMatrix getW(SP::DynamicalSystem ds = SP::DynamicalSystem())

get the value of W corresponding to DynamicalSystem ds

Parameters

ds – a pointer to DynamicalSystem, optional, default = nullptr. get W[0] in that case

Returns
SP::SimpleMatrix W(SP::DynamicalSystem ds)

get W corresponding to DynamicalSystem ds

Parameters

ds – a pointer to DynamicalSystem

Returns

pointer to a SiconosMatrix

const SimpleMatrix getWBoundaryConditions(SP::DynamicalSystem ds = SP::DynamicalSystem())

get the value of WBoundaryConditions corresponding to DynamicalSystem ds

Parameters

ds – a pointer to DynamicalSystem, optional, default = nullptr. get WBoundaryConditions[0] in that case

Returns
SP::SiconosMatrix WBoundaryConditions(SP::DynamicalSystem ds)

get WBoundaryConditions corresponding to DynamicalSystem ds

Parameters

ds – a pointer to DynamicalSystem, optional, default = nullptr. get WBoundaryConditions[0] in that case

Returns

pointer to a SiconosMatrix

inline double theta()

get theta

Returns

a double

inline void setTheta(double newTheta)

set the value of theta

Parameters

newTheta – a double

inline double gamma()

get gamma

Returns

a double

inline void setGamma(double newGamma)

set the value of gamma

Parameters

newGamma – a double

inline bool useGamma()

get bool useGamma

Returns

a bool

inline void setUseGamma(bool b)

set the boolean to indicate that we use gamma

Parameters

b – true if gamma has to be used, false otherwise

inline bool useGammaForRelation()

get bool gammaForRelation for the relation

Returns

a

inline void setUseGammaForRelation(bool newUseGammaForRelation)

set the boolean to indicate that we use gamma for the relation

Parameters

newUseGammaForRelation – a bool

virtual void initializeWorkVectorsForDS(double t, SP::DynamicalSystem ds) override

initialization of the work vectors and matrices (properties) related to one dynamical system on the graph and needed by the osi

Parameters
• t – time of initialization

• ds – the dynamical system

virtual void initializeWorkVectorsForInteraction(Interaction &inter, InteractionProperties &interProp, DynamicalSystemsGraph &DSG) override

initialization of the work vectors and matrices (properties) related to one interaction on the graph and needed by the osi

Parameters
• inter – the interaction

• interProp – the properties on the graph

• DSG – the dynamical systems graph

inline virtual unsigned int numberOfIndexSets() const override

get the number of index sets required for the simulation

Returns

unsigned int

void initializeIterationMatrixW(double time, SP::DynamicalSystem ds)

initialize iteration matrix W EulerMoreauOSI matrix at time t

Parameters
void computeW(double time, DynamicalSystem &ds, DynamicalSystemsGraph::VDescriptor &dsv, SiconosMatrix &W)

compute W EulerMoreauOSI matrix at time t

Parameters
• time – the current time

• ds – the DynamicalSystem

• dsv – a descriptor of the ds on the graph (redundant to avoid invocation)

• W – the matrix to compute

void computeKhat(Interaction &inter, SiconosMatrix &m, VectorOfSMatrices &workM, double h) const
void computeWBoundaryConditions(SP::DynamicalSystem ds)

compute WBoundaryConditionsMap[ds] EulerMoreauOSI matrix at time t

Parameters

ds – a pointer to DynamicalSystem

void initializeIterationMatrixWBoundaryConditions(SP::DynamicalSystem ds)

initialize iteration matrix WBoundaryConditionsMap[ds] EulerMoreauOSI

Parameters

ds – a pointer to DynamicalSystem

virtual double computeResidu() override

Computes the residuFree and residu of all the DynamicalSystems.

Returns

the maximum of the 2-norm over all the residu

virtual void computeFreeState() override

Perform the integration of the dynamical systems linked to this integrator without taking into account the nonsmooth input r.

virtual void updateOutput(double time) override

update the output of the Interaction attached to this Integrator

virtual void updateInput(double time) override

update the input of the Interaction attached to this Integrator

virtual void updateOutput(double time, unsigned int level) override

update the output of the Interaction attached to this Integrator

Parameters
• time – current time

• level – level of interest for the dynamics

virtual void updateInput(double time, unsigned int level) override

update the input of the Interaction attached to this Integrator

Parameters
• time – current time

• level – level of interest for the dynamics

virtual double computeResiduOutput(double time, SP::InteractionsGraph indexSet) override

compute the residu of the output of the relation (y) This computation depends on the type of OSI

Parameters
• time – time of computation

• indexSet – the index set of the interaction that are concerned

virtual double computeResiduInput(double time, SP::InteractionsGraph indexSet) override

compute the residu of the input of the relation (R or p) This computation depends on the type of OSI

Parameters
• time – time of computation

• indexSet – the index set of the interaction that are concerned

virtual void computeFreeOutput(InteractionsGraph::VDescriptor &vertex_inter, OneStepNSProblem *osnsp) override

integrates the Interaction linked to this integrator, without taking non-smooth effects into account

Parameters
• vertex_inter – of the interaction graph

• osnsp – pointer to OneStepNSProblem

virtual void prepareNewtonIteration(double time) override

computes all the W matrices

Parameters

time – current time

virtual void integrate(double &tinit, double &tend, double &tout, int &useless) override

integrate the system, between tinit and tend (->iout=true), with possible stop at tout (->iout=false)

Parameters
• tinit – initial time

• tend – end time

• tout – real end time

• useless – flag (for EulerMoreauOSI, used in LsodarOSI)

virtual void updateState(const unsigned int level) override

updates the state of the Dynamical Systems

Parameters

level – the level of interest for the dynamics: not used at the time

virtual void display() override

Displays the data of the EulerMoreauOSI’s integrator.

ACCEPT_STD_VISITORS()

Protected Functions

ACCEPT_SERIALIZATION(EulerMoreauOSI)
inline EulerMoreauOSI()

Default constructor.

Protected Attributes

double _theta

Stl map that associates a theta parameter for the integration scheme to each DynamicalSystem of the OSI.

double _gamma

A gamma parameter for the integration scheme to each DynamicalSystem of the OSI This parameter is used to apply a theta-method to the input $r$.

bool _useGamma

a boolean to know if the parameter must be used or not

bool _useGammaForRelation

a boolean to know if the parameter must be used or not

Friends

friend struct _NSLEffectOnFreeOutput