MoreauJeanOSI.hpp#

Return to documentation for this file
#ifndef MoreauJeanOSI_H
#define MoreauJeanOSI_H

#include "OneStepIntegrator.hpp"

#include <limits>

const unsigned int MOREAUSTEPSINMEMORY = 1;

class MoreauJeanOSI : public OneStepIntegrator {
protected:
  ACCEPT_SERIALIZATION(MoreauJeanOSI);

  double _theta;

  double _gamma;

  bool _useGamma;

  double _constraintActivationThreshold;

  bool _useGammaForRelation;

  bool _explicitNewtonEulerDSOperators;

  bool _isWSymmetricDefinitePositive;

  bool _activateWithNegativeRelativeVelocity;

  double _constraintActivationThresholdVelocity;

  std::vector<std::size_t> _selected_coordinates;

  struct _NSLEffectOnFreeOutput : public SiconosVisitor {
    using SiconosVisitor::visit;

    OneStepNSProblem &_osnsp;
    Interaction &_inter;
    InteractionProperties &_interProp;
    double _theta;

    _NSLEffectOnFreeOutput(OneStepNSProblem &p, Interaction &inter,
                           InteractionProperties &interProp,
                       double theta)
      : _osnsp(p), _inter(inter), _interProp(interProp), _theta(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);
  };

  friend struct _NSLEffectOnFreeOutput;

public:
  enum MoreauJeanOSI_ds_workVector_id {
    RESIDU_FREE,
    VFREE,
    BUFFER,
    QTMP,
    WORK_LENGTH
  };

  enum MoreauJeanOSI_interaction_workVector_id {
    OSNSP_RHS,
    WORK_INTERACTION_LENGTH
  };

  enum MoreauJeanOSI_interaction_workBlockVector_id {
    xfree,
    BLOCK_WORK_LENGTH
  };

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

  virtual ~MoreauJeanOSI(){};

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

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

  inline bool isWSymmetricDefinitePositive() const
  {
    return _isWSymmetricDefinitePositive;
  };

  inline void setIsWSymmetricDefinitePositive(bool b)
  {
    _isWSymmetricDefinitePositive = b;
  };

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

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

  inline double theta() { return _theta; };

  inline void setTheta(double newTheta) { _theta = newTheta; };

  inline double gamma() { return _gamma; };

  inline void setGamma(double newGamma)
  {
    _gamma = newGamma;
    _useGamma = true;
  };

  inline bool useGamma() { return _useGamma; };

  inline void setUseGamma(bool newUseGamma) { _useGamma = newUseGamma; };

  inline bool useGammaForRelation() { return _useGammaForRelation; };

  inline void setUseGammaForRelation(bool newUseGammaForRelation)
  {
    _useGammaForRelation = newUseGammaForRelation;
    if (_useGammaForRelation)
      _useGamma = false;
  };

  inline void setConstraintActivationThreshold(double v)
  {
    _constraintActivationThreshold = v;
  }

  inline double constraintActivationThreshold()
  {
    return _constraintActivationThreshold;
  }

  inline void setConstraintActivationThresholdVelocity(double v)
  {
    _constraintActivationThresholdVelocity = v;
  }

  inline double constraintActivationThresholdVelocity()
  {
    return _constraintActivationThresholdVelocity;
  }

  inline bool explicitNewtonEulerDSOperators()
  {
    return _explicitNewtonEulerDSOperators;
  };

  inline void
  setExplicitNewtonEulerDSOperators(bool newExplicitNewtonEulerDSOperators)
  {
    _explicitNewtonEulerDSOperators = newExplicitNewtonEulerDSOperators;
  };

  inline bool activateWithNegativeRelativeVelocity()
  {
    return _activateWithNegativeRelativeVelocity;
  };

  inline void
  setActivateWithNegativeRelativeVelocity(bool newActivateWithNegativeRelativeVelocity)
  {
    _activateWithNegativeRelativeVelocity = newActivateWithNegativeRelativeVelocity;
  };

  void initialize_nonsmooth_problems() override;

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

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

  unsigned int numberOfIndexSets() const override { return 2; };

  void initializeIterationMatrixW(double time, SP::SecondOrderDS ds);

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

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

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

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

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

  void computeInitialNewtonState() override;

  double computeResidu() override;

  void computeFreeState() override;

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

  SiconosVector& osnsp_rhs(InteractionsGraph::VDescriptor& vertex_inter,   InteractionsGraph& indexSet) override
  {
    return *(*indexSet.properties(vertex_inter).workVectors)[MoreauJeanOSI::OSNSP_RHS];
  };

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

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

  void prepareNewtonIteration(double time) override;

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

  virtual void updatePosition(DynamicalSystem &ds);

  void updateState(const unsigned int level) override;

  SP::SimpleMatrix computeWorkForces();

  void display() override;

  ACCEPT_STD_VISITORS();
};

#endif