File kernel/src/simulationTools/MoreauJeanOSI.hpp

Go to the source code of this file

Variables

const unsigned int MOREAUSTEPSINMEMORY = 1

class MoreauJeanOSI : public OneStepIntegrator

#include <MoreauJeanOSI.hpp>

One Step time Integrator, Moreau-Jean algorithm.

This integrator is the work horse of the event&#8212;capturing time stepping schemes for mechanical systems. It is mainly based on the pioneering works of M. Jean and J.J. Moreau for the time integration of mechanical systems with unilateral contact, impact and Coulomb’s friction with \( \theta \) scheme

For the linear Lagrangian system, the scheme reads as

\[\begin{split} \begin{cases} M (v_{k+1}-v_k) + h K q_{k+\theta} + h C v_{k+\theta} - h F_{k+\theta} = p_{k+1} = G P_{k+1},\label{eq:MoreauTS-motion}\\[1mm] q_{k+1} = q_{k} + h v_{k+\theta}, \quad \\[1mm] U_{k+1} = G^\top\, v_{k+1}, \\[1mm] \begin{array}{lcl} 0 \leq U^\alpha_{k+1} + e U^\alpha_{k} \perp P^\alpha_{k+1} \geq 0,& \quad&\alpha \in \mathcal I_1, \\[1mm] P^\alpha_{k+1} =0,&\quad& \alpha \in \mathcal I \setminus \mathcal I_1, \end{array} \end{cases} \end{split}\]

with \( \theta \in [0,1] \). The index set \( \mathcal I_1 \) is the discrete equivalent to the rule that allows us to apply the Signorini condition at the velocity level. In the numerical practice, we choose to define this set by

\[ \mathcal I_1 = \{\alpha \in \mathcal I \mid G^\top (q_{k} + h v_{k}) + w \leq 0\text{ and } U_k \leq 0 \}. \]

For more details, we refer to

M. Jean and J.J. Moreau. Dynamics in the presence of unilateral contacts and dry friction: a numerical approach. In G. Del Pietro and F. Maceri, editors, Unilateral problems in structural analysis. II, pages 151–196. CISM 304, Spinger Verlag, 1987.

J.J. Moreau. Unilateral contact and dry friction in finite freedom dynamics. In J.J. Moreau and Panagiotopoulos P.D., editors, Nonsmooth Mechanics and Applications, number 302 in CISM, Courses and lectures, pages 1–82. CISM 302, Spinger Verlag, Wien- New York, 1988a.

J.J. Moreau. Numerical aspects of the sweeping process. Computer Methods in Applied Mechanics and Engineering, 177:329–349, 1999.

M. Jean. The non smooth contact dynamics method. Computer Methods in Applied Mechanics and Engineering, 177:235–257, 1999.

and for a review :

V. Acary and B. Brogliato. Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics, volume 35 of Lecture Notes in Applied and Computational Mechanics. Springer Verlag, 2008.

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

Each DynamicalSystem is associated to a SiconosMatrix, named “W”, the “teration” 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 mechanical systems, the implementation uses _p for storing the the input due to the nonsmooth law. This MoreauJeanOSI scheme assumes that the relative degree is two.

For Lagrangian systems, the implementation uses _p[1] for storing the discrete impulse.

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.

Subclassed by MoreauJeanCombinedProjectionOSI, MoreauJeanDirectProjectionOSI, MoreauJeanGOSI

Public Types

enum MoreauJeanOSI_ds_workVector_id

Values:

enumerator RESIDU_FREE

enumerator VFREE

enumerator BUFFER

enumerator QTMP

enumerator WORK_LENGTH

enum MoreauJeanOSI_interaction_workVector_id

Values:

enumerator OSNSP_RHS

enumerator WORK_INTERACTION_LENGTH

enum MoreauJeanOSI_interaction_workBlockVector_id

Values:

enumerator xfree

enumerator BLOCK_WORK_LENGTH

Public Functions

MoreauJeanOSI(double theta = 0.5, double gamma = std::numeric_limits<double>::quiet_NaN())

constructor from theta value only

Parameters:

• theta – value for all linked DS (default = 0.5).

• gamma – value for all linked DS (default = NaN and gamma is not used).

inline virtual ~MoreauJeanOSI()

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:

SimpleMatrix

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

get W corresponding to DynamicalSystem ds

Parameters:

ds – a pointer to DynamicalSystem

Returns:

pointer to a SiconosMatrix

inline bool isWSymmetricDefinitePositive() const

inline void setIsWSymmetricDefinitePositive(bool b)

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:

SimpleMatrix

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 newUseGamma)

set the Boolean to indicate that we use gamma

Parameters:

newUseGamma – a Boolean variable

inline bool useGammaForRelation()

get bool gammaForRelation for the relation

Returns:

a Boolean

inline void setUseGammaForRelation(bool newUseGammaForRelation)

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

Parameters:

newUseGammaForRelation – a Boolean

inline void setConstraintActivationThreshold(double v)

set the constraint activation threshold

inline double constraintActivationThreshold()

get the constraint activation threshold

inline void setConstraintActivationThresholdVelocity(double v)

set the constraint activation threshold

inline double constraintActivationThresholdVelocity()

get the constraint activation threshold

inline bool explicitNewtonEulerDSOperators()

get boolean _explicitNewtonEulerDSOperators for the relation

Returns:

a Boolean

inline void setExplicitNewtonEulerDSOperators(bool newExplicitNewtonEulerDSOperators)

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

Parameters:

newExplicitNewtonEulerDSOperators – a Boolean

inline bool activateWithNegativeRelativeVelocity()

get boolean _activateWithNegativeRelativeVelocity

Returns:

a Boolean

inline void setActivateWithNegativeRelativeVelocity(bool newActivateWithNegativeRelativeVelocity)

set the boolean to perform activation with negative relative velocity

Parameters:

newActivateWithNegativeRealtiveVelocity – a Boolean

virtual void initialize_nonsmooth_problems() override

initialization of the MoreauJeanOSI integrator; for linear time invariant systems, we compute time invariant operator (example : W)

Initialization process of the nonsmooth problems linked to this OSI

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::SecondOrderDS ds)

initialize iteration matrix W MoreauJeanOSI matrix at time t

Parameters:

• time –

• ds – a pointer to DynamicalSystem

void computeW(double time, SecondOrderDS &ds, SiconosMatrix &W)

compute W MoreauJeanOSI matrix at time t

Parameters:

• time – (double)

• ds – a DynamicalSystem

• W – the result in W

SP::SimpleMatrix Winverse(SP::SecondOrderDS ds, bool keepW = false)

get and compute if needed W MoreauJeanOSI matrix at time t

Parameters:

• time – (double)

• ds – a DynamicalSystem

• Winverse – the result in Winverse

• keepW –

void _computeWBoundaryConditions(SecondOrderDS &ds, SiconosMatrix &WBoundaryConditions, SiconosMatrix &iteration_matrix)

compute WBoundaryConditionsMap[ds] MoreauJeanOSI matrix at time t

Parameters:

• ds – a pointer to DynamicalSystem

• WBoundaryConditions – write the result in WBoundaryConditions

• iteration_matrix – the OSI iteration matrix (W)

void _initializeIterationMatrixWBoundaryConditions(SecondOrderDS &ds, const DynamicalSystemsGraph::VDescriptor &dsv)

initialize iteration matrix WBoundaryConditionsMap[ds] MoreauJeanOSI

Parameters:

• ds – a pointer to DynamicalSystem

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

void applyBoundaryConditions(SecondOrderDS &d, SiconosVector &residu, DynamicalSystemsGraph::VIterator dsi, double t, const SiconosVector &v)

virtual void computeInitialNewtonState() override

compute the initial state of the Newton loop.

virtual double computeResidu() override

return the maximum of all norms for the “MoreauJeanOSI-discretized” residus of DS

Returns:

a double

virtual void computeFreeState() override

Perform the integration of the dynamical systems linked to this integrator without taking into account the nonsmooth input (_r or _p)

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 – vertex of the interaction graph

• osnsp – pointer to OneStepNSProblem

inline virtual SiconosVector &osnsp_rhs(InteractionsGraph::VDescriptor &vertex_inter, InteractionsGraph &indexSet) override

return the workVector corresponding to the right hand side of the OneStepNonsmooth problem

virtual bool addInteractionInIndexSet(SP::Interaction inter, unsigned int i) override

Apply the rule to one Interaction to know if it should be included in the IndexSet of level i.

Parameters:

• inter – the Interaction to test

• i – level of the IndexSet

Returns:

Boolean

virtual bool removeInteractionFromIndexSet(SP::Interaction inter, unsigned int i) override

Apply the rule to one Interaction to know if it should be removed from the IndexSet of level i.

Parameters:

• inter – the Interaction to test

• i – level of the IndexSet

Returns:

Boolean

virtual void prepareNewtonIteration(double time) override

method to prepare the fist Newton iteration

Parameters:

time –

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

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

Parameters:

• tinit – the initial time

• tend – the end time

• tout – the real end time

• notUsed – useless flag (for MoreauJeanOSI, used in LsodarOSI)

virtual void updatePosition(DynamicalSystem &ds)

update the state of the dynamical systems

Parameters:

ds – the dynamical to update

virtual void updateState(const unsigned int level) override

update the state of the dynamical systems

Parameters:

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

SP::SimpleMatrix computeWorkForces()

Compute the matrix of work of forces by ds.

Returns:

SP::Siconosmatrix

virtual void display() override

Displays the data of the MoreauJeanOSI’s integrator.

ACCEPT_STD_VISITORS()

Protected Functions

ACCEPT_SERIALIZATION(MoreauJeanOSI)

Protected Attributes

double _theta

theta-scheme parameter

double _gamma

A gamma parameter for the forecast of activation of constraints leap-frog estimation of the constraints \( \tilde y_k = y_k + \gamma * h * ydot \).

bool _useGamma

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

double _constraintActivationThreshold

Constraint activation threshold.

bool _useGammaForRelation

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

bool _explicitNewtonEulerDSOperators

a boolean to force the evaluation of T in an explicit way

bool _isWSymmetricDefinitePositive

a boolean to know if the matrix W is symmetric definite positive

bool _activateWithNegativeRelativeVelocity

a boolean to perform activation with negative relative velocity

double _constraintActivationThresholdVelocity

Constraint activation threshold.

std::vector<std::size_t> _selected_coordinates

A set of work indices for the selected coordinates when we subprod in computeFreeOuput.

Friends

friend struct _NSLEffectOnFreeOutput

struct _NSLEffectOnFreeOutput : public SiconosVisitor

#include <MoreauJeanOSI.hpp>

nslaw effects

Public Functions

inline _NSLEffectOnFreeOutput(OneStepNSProblem &p, Interaction &inter, InteractionProperties &interProp, double theta)

void visit(const NewtonImpactNSL &nslaw)

void visit(const RelayNSL &nslaw)

void visit(const NewtonImpactFrictionNSL &nslaw)

void visit(const FremondImpactFrictionNSL &nslaw)

void visit(const NewtonImpactRollingFrictionNSL &nslaw)

void visit(const EqualityConditionNSL &nslaw)

void visit(const MixedComplementarityConditionNSL &nslaw)

void visit(const ComplementarityConditionNSL &nslaw)

Public Members

OneStepNSProblem &_osnsp

Interaction &_inter

InteractionProperties &_interProp

double _theta