1 /* Siconos is a program dedicated to modeling, simulation and control
2  * of non smooth dynamical systems.
3  *
4  * Copyright 2021 INRIA.
5  *
6  * Licensed under the Apache License, Version 2.0 (the "License");
7  * you may not use this file except in compliance with the License.
8  * You may obtain a copy of the License at
9  *
10  * http://www.apache.org/licenses/LICENSE-2.0
11  *
12  * Unless required by applicable law or agreed to in writing, software
13  * distributed under the License is distributed on an "AS IS" BASIS,
14  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15  * See the License for the specific language governing permissions and
16  * limitations under the License.
17 */
18 /*! \file OneStepNSProblem.hpp
19   \brief Interface to formalize and solve Non-Sooth problems
20  */
21 
22 #ifndef ONESTEPNSPROBLEM_H
23 #define ONESTEPNSPROBLEM_H
24 
25 #include "SiconosFwd.hpp"
26 #include "SiconosVisitor.hpp"
27 #include "SimulationTypeDef.hpp"
28 #include "SimulationGraphs.hpp"
29 
30 /** Non Smooth Problem Formalization and Simulation
31 
32   This is an abstract class, that provides an interface to define a
33   non smooth problem:
34   - a formulation (ie the way the problem is written)
35   - a solver (algorithm and solving formulation, that can be
36         different from problem formulation)
37   - routines to compute the problem solution.
38 
39   Two types of problem formulation are available :
40    - Quadratic Problem
41    - Linear Problem
42 
43    See derived classes (QP and LinearOSNS) for details.
44 
45    For Linear problems, the following formulations exists:
46    - Linear Complementarity (LCP)
47    - Mixed Linear Complementarity (MLCP)
48    - Affine Variational Inequalities (AVI)
49    - FrictionContact
50    - Relay
51    - Equality
52    - GenericMechanical
53    - MultipleImpact
54    - GlobalFrictionContact
55 
56    The usual way to build and initialize a one-step nonsmooth problem is :
57    - call constructor with the id of the required Numerics solver.
58    (see Solver class or Numerics documentation for details on algorithm name and parameters).
59    - initialize(simulation)
60    Initialize process is usually done through model->initialize(simulation).
61    See Examples for practical details.
62 
63    \section osns_options Options for Numerics and the driver for solvers
64 
65    When the Numerics driver is called a set of solver options (name, tolerance, max. number of iterations ...)
66    is required --> SolverOptions.
67 
68    Default values are always set in solver options the OneStepNSProblem is built
69    but if you need to set them yourself, please check Users'guide, Numerics solvers part.
70 
71  */
72 class OneStepNSProblem
73 {
74 
75 protected:
76   /* serialization hooks */
77   ACCEPT_SERIALIZATION(OneStepNSProblem);
78 
79   /** Numerics solver properties */
80   SP::SolverOptions _numerics_solver_options;
81 
82   /** size of the nonsmooth problem */
83   unsigned int _sizeOutput = 0;
84 
85   /** link to the simulation that owns the nonsmooth problem */
86   SP::Simulation _simulation;
87 
88   /** level of index sets that is considered by this osnsp */
89   unsigned int _indexSetLevel = 0;
90 
91   /** level of input and output variables of the nonsmooth problems.
92    *  We consider that the osnsp computes y[_inputOutputLevel] and lambda[_inputOutputLevel]
93    */
94   unsigned int _inputOutputLevel = 0;
95 
96   /** maximum value for sizeOutput. Set to the number of declared
97       constraints by default (topology->getNumberOfConstraints());
98       This value is used to allocate memory for M during initialize
99       call. The best choice is to set maxSize to the estimated maximum
100       dimension of the problem. It must not exceed ...
101   */
102   unsigned int _maxSize = 0;
103 
104   /* set of nslaw types */
105   std::set<float> _nslawtype;
106 
107   /*During Newton it, this flag allows to update the numerics matrices only once if necessary.*/
108   bool _hasBeenUpdated = false;
109 
110   // --- CONSTRUCTORS/DESTRUCTOR ---
111   /** default constructor */
112   OneStepNSProblem() = default;
113 
114 private:
115 
116   /* copy constructor, forbidden */
117   OneStepNSProblem(const OneStepNSProblem&) = delete;
118 
119   /* assignment, forbidden */
120   OneStepNSProblem& operator=(const OneStepNSProblem&) = delete;
121 
122 public:
123   /**  constructor from a pre-defined solver options set.
124        \param options, the options set,
125        \rst
126        see :ref:`problems_and_solvers` for details.
127        \endrst
128   */
OneStepNSProblem(SP::SolverOptions options)129   OneStepNSProblem(SP::SolverOptions options): _numerics_solver_options(options){};
130 
131   /** destructor
132    */
~OneStepNSProblem()133   virtual ~OneStepNSProblem(){};
134 
135   // --- GETTERS/SETTERS ---
136 
137 
138   /** To get the SolverOptions structure
139    *  \return , the numerics structure used to save solver parameters
140    */
numericsSolverOptions() const141   inline SP::SolverOptions numericsSolverOptions() const
142   {
143     return _numerics_solver_options;
144   };
145 
146   /** returns the dimension of the nonsmooth problem
147    */
getSizeOutput() const148   inline unsigned int getSizeOutput() const
149   {
150     return _sizeOutput;
151   }
152 
153   /** get the simulation which owns this nonsmooth problem
154    *  \return a pointer on Simulation
155    */
simulation() const156   inline SP::Simulation simulation() const
157   {
158     return _simulation;
159   }
160 
161   /** set the Simulation of the OneStepNSProblem
162    *  \param newS a pointer to Simulation
163    */
setSimulationPtr(SP::Simulation newS)164   inline void setSimulationPtr(SP::Simulation newS)
165   {
166     _simulation = newS;
167   }
168 
169   /** get indexSetLevel
170    *  \return an unsigned int
171    */
indexSetLevel() const172   inline unsigned int indexSetLevel() const
173   {
174     return _indexSetLevel;
175   }
176 
177   /** set the value of level min
178    *  \param newVal an unsigned int
179    */
setIndexSetLevel(unsigned int newVal)180   inline void setIndexSetLevel(unsigned int newVal)
181   {
182     _indexSetLevel = newVal;
183   }
184 
185   /** get the Input/Output level
186    *  \return an unsigned int
187    */
inputOutputLevel() const188   inline unsigned int inputOutputLevel() const
189   {
190     return _inputOutputLevel;
191   }
192 
193   /** set the value of Input/Output level
194    *  \param newVal an unsigned int
195    */
setInputOutputLevel(unsigned int newVal)196   inline void setInputOutputLevel(unsigned int newVal)
197   {
198     _inputOutputLevel = newVal;
199   }
200 
201   /** get maximum value allowed for the dimension of the problem
202    *  \return an unsigned int
203    */
maxSize() const204   inline unsigned int maxSize() const
205   {
206     return _maxSize;
207   }
208 
209   /** set the value of maxSize
210    *  \param newVal an unsigned int
211    */
setMaxSize(const unsigned int newVal)212   inline void setMaxSize(const unsigned int newVal)
213   {
214     _maxSize = newVal;
215   }
216 
217   /** Turn on/off verbose mode in numerics solver*/
218   void setNumericsVerboseMode(bool vMode);
219 
220 
221   /**  set the verbose level in numerics solver*/
222   void setNumericsVerboseLevel(int level);
223 
224   /** Check if the OSNSPb has interactions.
225       \return bool = true if the  osnsp has interactions, i.e. indexSet(_indexSetLevel)->size >0
226    */
227   bool hasInteractions() const;
228 
display() const229   virtual void display() const {}
230 
231   /** Display the set of blocks for  a given indexSet
232    * \param  indexSet  the concerned index set
233    */
234   virtual void displayBlocks(SP::InteractionsGraph indexSet);
235 
236   /** compute interactionBlocks if necessary (this depends on the type of
237    * OSNS, on the indexSets ...)
238    */
239   virtual void updateInteractionBlocks();
240 
241   /** compute extra-diagonal interactionBlock-matrix
242    *  \param ed an edge descriptor
243    */
244   virtual void computeInteractionBlock(const InteractionsGraph::EDescriptor& ed ) = 0;
245 
246   /** compute diagonal Interaction block
247    * \param vd a vertex descriptor
248    */
249   virtual void computeDiagonalInteractionBlock(const InteractionsGraph::VDescriptor& vd) = 0;
250 
251   /**
252    * \return bool _hasBeenUpdated
253    */
hasBeenUpdated()254   bool hasBeenUpdated()
255   {
256     return _hasBeenUpdated;
257   }
258 
259   /**
260    * \param v to set _hasBeenUpdated.
261    */
setHasBeenUpdated(bool v)262   void setHasBeenUpdated(bool v)
263   {
264     _hasBeenUpdated = v;
265   }
266 
267   /** initialize the problem(compute topology ...)
268       \param sim the simulation, owner of this OSNSPB
269     */
270   virtual void initialize(SP::Simulation sim);
271 
272   /** prepare data of the osns for solving
273    *  \param time the current time
274    *  \return true if the computation of the OSNS has to be carry on, false otherwise
275    */
276   virtual bool preCompute(double time) = 0;
277 
278   /** To run the solver for ns problem
279    *  \param time  current time
280    *  \return int information about the solver convergence.
281    */
282   virtual int compute(double time) = 0;
283 
284   /** post treatment for output of the solver
285    */
286   virtual void postCompute() = 0;
287 
288   /** change the solver type and its default parameters.
289 
290       - clear memory for the existing options set
291       - create and initialize a new one
292 
293       \rst
294       Check the available solvers id in :ref:`problems_and_solvers`.
295       \endrst
296 
297       \param solverId the new solver.
298    */
299   void setSolverId(int solverId);
300 
301   /** get the OSI-related matrices used to compute the current InteractionBlock
302       (Ex: for MoreauJeanOSI, W)
303       \param osi the OSI of the concerned dynamical system
304       \param ds the concerned dynamical system
305       \return the required matrix.
306   */
307   SP::SimpleMatrix getOSIMatrix(OneStepIntegrator& osi, SP::DynamicalSystem ds);
308 
309   /* visitors hook */
310   ACCEPT_STD_VISITORS();
311 
312 };
313 
314 #endif // ONESTEPNSPROBLEM_H
315