Program listing for file kernel/src/simulationTools/MoreauJeanOSI.hpp

Program listing for file kernel/src/simulationTools/MoreauJeanOSI.hpp#

  1#ifndef MoreauJeanOSI_H
  2#define MoreauJeanOSI_H
  3
  4#include "OneStepIntegrator.hpp"
  5
  6#include <limits>
  7
  8const unsigned int MOREAUSTEPSINMEMORY = 1;
  9
 10
 11class MoreauJeanOSI : public OneStepIntegrator {
 12protected:
 13  ACCEPT_SERIALIZATION(MoreauJeanOSI);
 14
 15
 16  double _theta;
 17
 18
 19  double _gamma;
 20
 21
 22  bool _useGamma;
 23
 24
 25  double _constraintActivationThreshold;
 26
 27
 28  bool _useGammaForRelation;
 29
 30
 31  bool _explicitNewtonEulerDSOperators;
 32
 33
 34  bool _isWSymmetricDefinitePositive;
 35
 36
 37  bool _activateWithNegativeRelativeVelocity;
 38
 39
 40  double _constraintActivationThresholdVelocity;
 41
 42
 43  std::vector<std::size_t> _selected_coordinates;
 44
 45
 46
 47  struct _NSLEffectOnFreeOutput : public SiconosVisitor {
 48    using SiconosVisitor::visit;
 49
 50    OneStepNSProblem &_osnsp;
 51    Interaction &_inter;
 52    InteractionProperties &_interProp;
 53    double _theta;
 54
 55    _NSLEffectOnFreeOutput(OneStepNSProblem &p, Interaction &inter,
 56                           InteractionProperties &interProp,
 57                       double theta)
 58      : _osnsp(p), _inter(inter), _interProp(interProp), _theta(theta) {};
 59
 60    void visit(const NewtonImpactNSL &nslaw);
 61    void visit(const RelayNSL &nslaw);
 62    void visit(const NewtonImpactFrictionNSL &nslaw);
 63    void visit(const FremondImpactFrictionNSL &nslaw);
 64    void visit(const NewtonImpactRollingFrictionNSL &nslaw);
 65    void visit(const EqualityConditionNSL &nslaw);
 66    void visit(const MixedComplementarityConditionNSL &nslaw);
 67    void visit(const ComplementarityConditionNSL &nslaw);
 68  };
 69
 70  friend struct _NSLEffectOnFreeOutput;
 71
 72public:
 73  enum MoreauJeanOSI_ds_workVector_id {
 74    RESIDU_FREE,
 75    VFREE,
 76    BUFFER,
 77    QTMP,
 78    WORK_LENGTH
 79  };
 80
 81  enum MoreauJeanOSI_interaction_workVector_id {
 82    OSNSP_RHS,
 83    WORK_INTERACTION_LENGTH
 84  };
 85
 86  enum MoreauJeanOSI_interaction_workBlockVector_id {
 87    xfree,
 88    BLOCK_WORK_LENGTH
 89  };
 90
 91
 92  MoreauJeanOSI(double theta = 0.5,
 93                double gamma = std::numeric_limits<double>::quiet_NaN());
 94
 95
 96  virtual ~MoreauJeanOSI(){};
 97
 98
 99
100
101  const SimpleMatrix getW(SP::DynamicalSystem ds = SP::DynamicalSystem());
102
103
104  SP::SimpleMatrix W(SP::DynamicalSystem ds);
105
106  inline bool isWSymmetricDefinitePositive() const
107  {
108    return _isWSymmetricDefinitePositive;
109  };
110
111  inline void setIsWSymmetricDefinitePositive(bool b)
112  {
113    _isWSymmetricDefinitePositive = b;
114  };
115
116
117
118
119  const SimpleMatrix
120  getWBoundaryConditions(SP::DynamicalSystem ds = SP::DynamicalSystem());
121
122
123  SP::SiconosMatrix WBoundaryConditions(SP::DynamicalSystem ds);
124
125
126
127
128  inline double theta() { return _theta; };
129
130
131  inline void setTheta(double newTheta) { _theta = newTheta; };
132
133
134
135
136  inline double gamma() { return _gamma; };
137
138
139  inline void setGamma(double newGamma)
140  {
141    _gamma = newGamma;
142    _useGamma = true;
143  };
144
145
146
147
148  inline bool useGamma() { return _useGamma; };
149
150
151  inline void setUseGamma(bool newUseGamma) { _useGamma = newUseGamma; };
152
153
154  inline bool useGammaForRelation() { return _useGammaForRelation; };
155
156
157  inline void setUseGammaForRelation(bool newUseGammaForRelation)
158  {
159    _useGammaForRelation = newUseGammaForRelation;
160    if (_useGammaForRelation)
161      _useGamma = false;
162  };
163
164  inline void setConstraintActivationThreshold(double v)
165  {
166    _constraintActivationThreshold = v;
167  }
168
169
170  inline double constraintActivationThreshold()
171  {
172    return _constraintActivationThreshold;
173  }
174
175  inline void setConstraintActivationThresholdVelocity(double v)
176  {
177    _constraintActivationThresholdVelocity = v;
178  }
179
180
181  inline double constraintActivationThresholdVelocity()
182  {
183    return _constraintActivationThresholdVelocity;
184  }
185
186
187  inline bool explicitNewtonEulerDSOperators()
188  {
189    return _explicitNewtonEulerDSOperators;
190  };
191
192
193  inline void
194  setExplicitNewtonEulerDSOperators(bool newExplicitNewtonEulerDSOperators)
195  {
196    _explicitNewtonEulerDSOperators = newExplicitNewtonEulerDSOperators;
197  };
198
199
200  inline bool activateWithNegativeRelativeVelocity()
201  {
202    return _activateWithNegativeRelativeVelocity;
203  };
204
205
206  inline void
207  setActivateWithNegativeRelativeVelocity(bool newActivateWithNegativeRelativeVelocity)
208  {
209    _activateWithNegativeRelativeVelocity = newActivateWithNegativeRelativeVelocity;
210  };
211
212
213
214
215
216
217
218
219  void initialize_nonsmooth_problems() override;
220
221
222  void initializeWorkVectorsForDS(double t, SP::DynamicalSystem ds) override;
223
224
225  void initializeWorkVectorsForInteraction(Interaction &inter,
226                                           InteractionProperties &interProp,
227                                           DynamicalSystemsGraph &DSG) override;
228
229
230  unsigned int numberOfIndexSets() const override { return 2; };
231
232
233  void initializeIterationMatrixW(double time, SP::SecondOrderDS ds);
234
235
236  void computeW(double time, SecondOrderDS &ds, SiconosMatrix &W);
237
238
239  SP::SimpleMatrix Winverse(SP::SecondOrderDS ds, bool keepW = false);
240
241
242  void _computeWBoundaryConditions(SecondOrderDS &ds,
243                                   SiconosMatrix &WBoundaryConditions,
244                                   SiconosMatrix &iteration_matrix);
245
246
247  void _initializeIterationMatrixWBoundaryConditions(
248      SecondOrderDS &ds, const DynamicalSystemsGraph::VDescriptor &dsv);
249
250  void applyBoundaryConditions(SecondOrderDS &d, SiconosVector &residu,
251                               DynamicalSystemsGraph::VIterator dsi, double t,
252                               const SiconosVector &v);
253
254
255  void computeInitialNewtonState() override;
256
257
258  double computeResidu() override;
259
260
261  void computeFreeState() override;
262
263
264  void computeFreeOutput(InteractionsGraph::VDescriptor &vertex_inter,
265                         OneStepNSProblem *osnsp) override;
266
267
268  SiconosVector& osnsp_rhs(InteractionsGraph::VDescriptor& vertex_inter,   InteractionsGraph& indexSet) override
269  {
270    return *(*indexSet.properties(vertex_inter).workVectors)[MoreauJeanOSI::OSNSP_RHS];
271  };
272
273
274  bool addInteractionInIndexSet(SP::Interaction inter, unsigned int i) override;
275
276
277  bool removeInteractionFromIndexSet(SP::Interaction inter,
278                                     unsigned int i) override;
279
280
281  void prepareNewtonIteration(double time) override;
282
283
284  void integrate(double &tinit, double &tend, double &tout,
285                 int &notUsed) override;
286
287
288  virtual void updatePosition(DynamicalSystem &ds);
289
290
291  void updateState(const unsigned int level) override;
292
293
294  SP::SimpleMatrix computeWorkForces();
295
296
297  void display() override;
298
299  ACCEPT_STD_VISITORS();
300};
301
302#endif