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 LinearOSNS.hpp 19 \brief Linear Complementarity Problem formulation and solving 20 */ 21 22 #ifndef LinearOSNS_H 23 #define LinearOSNS_H 24 25 #include "OneStepNSProblem.hpp" 26 #include "SiconosVector.hpp" 27 #include "NumericsMatrix.h" // For NM_DENSE 28 29 /** stl vector of double */ 30 typedef std::vector<double> MuStorage; 31 TYPEDEF_SPTR(MuStorage) 32 33 /** Base (abstract) class for linear non-smooth problems 34 35 Base class for linear non-smooth problems, usually in the form: 36 37 \f$ w = q + M z \f$ 38 39 where 40 - \f$ w \in R^{n} \f$ and \f$z \in R^{n} \f$ are the unknowns, 41 - \f$ M \in R^{n \times n } \f$ and \f$q \in R^{n} \f$ 42 43 examples: LCP, FrictionContact ... 44 45 */ 46 class LinearOSNS : public OneStepNSProblem 47 { 48 protected: 49 /* serialization hooks */ 50 ACCEPT_SERIALIZATION(LinearOSNS); 51 52 /** vector w of a LinearOSNS system */ 53 SP::SiconosVector _w; 54 55 /** vector z of a LinearOSNS system */ 56 SP::SiconosVector _z; 57 58 /** matrix M of a LinearOSNS system */ 59 SP::OSNSMatrix _M; 60 61 /** vector q of a LinearOSNS system */ 62 SP::SiconosVector _q; 63 64 /** Storage type for M - NM_DENSE: SiconosMatrix (dense), NM_SPARSE_BLOCK: Sparse Storage 65 (embedded into OSNSMatrix) */ 66 NM_types _numericsMatrixStorageType = NM_DENSE; 67 68 /** a boolean to decide if _w and _z vectors are initialized with 69 previous values of Y and Lambda when a change occurs in problem 70 size */ 71 bool _keepLambdaAndYState = true; 72 73 /** nslaw effects : visitors experimentation 74 */ 75 struct _TimeSteppingNSLEffect; 76 struct _EventDrivenNSLEffect; 77 struct _NSLEffectOnSim; 78 friend struct _TimeSteppingNSLEffect; 79 friend struct _EventDrivenNSLEffect; 80 friend struct _NSLEffectOnSim; 81 82 /** default constructor (private) 83 */ 84 LinearOSNS() = default ; 85 86 public: 87 88 /** constructor from a pre-defined solver options set. 89 \param options, the options set, 90 \rst 91 see :ref:`problems_and_solvers` for details. 92 \endrst 93 */ LinearOSNS(SP::SolverOptions options)94 LinearOSNS(SP::SolverOptions options): OneStepNSProblem(options){}; 95 96 /** destructor 97 */ ~LinearOSNS()98 virtual ~LinearOSNS() {}; 99 100 // --- W --- 101 /** copy of the current value of vector w 102 \return a SiconosVector 103 */ getW() const104 inline const SiconosVector getW() const 105 { 106 return *_w; 107 } 108 109 /** current w vector (pointer link) 110 \return pointer on a SiconosVector 111 */ w() const112 inline SP::SiconosVector w() const 113 { 114 return _w; 115 } 116 117 /** set w vector (pointer link) 118 \param newPtr the new SP::SiconosVector 119 */ setWPtr(SP::SiconosVector newPtr)120 inline void setWPtr(SP::SiconosVector newPtr) 121 { 122 _w = newPtr; 123 } 124 125 // --- Z --- 126 /** copy of the current value of vector z 127 \return a SiconosVector 128 */ getz() const129 inline const SiconosVector getz() const 130 { 131 return *_z; 132 } 133 134 /** current z vector (pointer link) 135 \return pointer on a SiconosVector 136 */ z() const137 inline SP::SiconosVector z() const 138 { 139 return _z; 140 } 141 142 /** set z vector (pointer link) 143 \param newPtr the new SP::SiconosVector 144 */ setzPtr(SP::SiconosVector newPtr)145 inline void setzPtr(SP::SiconosVector newPtr) 146 { 147 _z = newPtr; 148 } 149 150 /** M matrix (pointer link) 151 \return pointer on a OSNSMatrix 152 */ M() const153 inline SP::OSNSMatrix M() const 154 { 155 return _M; 156 } 157 158 /** set M to pointer newPtr 159 \param newM the new M matrix 160 */ setMPtr(SP::OSNSMatrix newM)161 inline void setMPtr(SP::OSNSMatrix newM) 162 { 163 _M = newM; 164 } 165 166 /** get the value of q, the constant vector in the LinearOSNS 167 \return SiconosVector 168 */ getQ() const169 inline const SiconosVector getQ() const 170 { 171 return *_q; 172 } 173 174 /** get q, the the constant vector in the LinearOSNS 175 \return pointer on a SiconosVector 176 */ q() const177 inline SP::SiconosVector q() const 178 { 179 return _q; 180 } 181 182 /** set q to pointer newPtr 183 \param newQ the new q vector 184 */ setQPtr(SP::SiconosVector newQ)185 inline void setQPtr(SP::SiconosVector newQ) 186 { 187 _q = newQ; 188 } 189 190 /** get the type of storage used for M 191 \return NM_types (NM_DENSE, NM_SPARSE_BLOCK) 192 */ getMStorageType() const193 inline NM_types getMStorageType() const 194 { 195 return _numericsMatrixStorageType; 196 }; 197 198 /** set which type of storage will be used for M 199 * \warning this function does not allocate any memory for M, 200 * it just sets an indicator for future use 201 * \param i (NM_DENSE, NM_SPARSE_BLOCK) 202 */ setMStorageType(NM_types i)203 inline void setMStorageType(NM_types i) 204 { 205 _numericsMatrixStorageType = i; 206 }; 207 208 /** Memory allocation or resizing for z,w,q */ 209 void initVectorsMemory(); 210 211 /** initialize the _M matrix */ 212 virtual void initOSNSMatrix(); 213 214 /** To initialize the LinearOSNS problem(computes topology ...) 215 \param sim the simulation owning this OSNSPB 216 */ 217 virtual void initialize(SP::Simulation sim); 218 219 /** compute extra-diagonal interactionBlock-matrix 220 * \param ed an edge descriptor 221 */ 222 virtual void computeInteractionBlock(const InteractionsGraph::EDescriptor& ed); 223 224 /** compute diagonal Interaction block 225 * \param vd a vertex descriptor 226 */ 227 virtual void computeDiagonalInteractionBlock(const InteractionsGraph::VDescriptor& vd); 228 229 /** To compute a part of the "q" vector of the OSNS 230 \param vertex, vertex (interaction) which corresponds to the considered block 231 \param pos the position of the first element of yOut to be set 232 */ 233 virtual void computeqBlock(InteractionsGraph::VDescriptor& vertex, unsigned int pos); 234 235 /** compute vector q 236 * \param time the current time 237 */ 238 virtual void computeq(double time); 239 240 /** build problem coefficients (if required) 241 \param time the current time 242 \return true if the indexSet is not empty 243 */ 244 virtual bool preCompute(double time); 245 246 /** Compute the unknown z and w and update the Interaction (y and lambda ) 247 * \param time the current time 248 * \return information about the solver convergence. 249 */ 250 virtual int compute(double time) = 0; 251 252 /** update interactions variables (y and lambda) according to current problem found solutions. 253 */ 254 virtual void postCompute(); 255 256 /** print the data to the screen 257 */ 258 virtual void display() const; 259 260 /** choose initialisation behavior for w and z. 261 \param val true: init w and z with previous values 262 of y and lambda saved in interactions, false: init to 0. 263 */ setKeepLambdaAndYState(bool val)264 void setKeepLambdaAndYState(bool val) 265 { 266 _keepLambdaAndYState = val ; 267 } 268 269 virtual bool checkCompatibleNSLaw(NonSmoothLaw& nslaw) =0; 270 271 /* visitors hook */ 272 ACCEPT_STD_VISITORS(); 273 274 275 }; 276 277 #endif // LinearOSNS_H 278