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 OneStepIntegrator.hpp 19 20 Base class (i.e. common interface) for dynamical system integration over a time step. 21 */ 22 23 #ifndef ONESTEPINTEGRATOR_H 24 #define ONESTEPINTEGRATOR_H 25 26 #include "SiconosVisitor.hpp" // for VIRTUAL_ACCEPT_VISITORS 27 #include "SiconosException.hpp" 28 #include "SimulationTypeDef.hpp" 29 #include "OneStepIntegratorTypes.hpp" 30 #include "SimulationGraphs.hpp" 31 32 /** Generic class to manage DynamicalSystem(s) time-integration 33 * 34 * !!! This is a virtual class, interface for some specific integrators !!! 35 * 36 * At the time, available integrators are: 37 * <ul> 38 * <li> EulerMoreauOSI </li> 39 * <li> MoreauJeanOSI </li> 40 * <li> MoreauJeanCombinedProjectionOSI </li> 41 * <li> MoreauJeanDirectProjectionOSI </li> 42 * <li> MoreauJeanBilbaoOSI </li> 43 * <li> D1MinusLinearOSI </li> 44 * <li> SchatzmanPaoliOSI </li> 45 * <li> LsodarOSI </li> 46 * <li> Hem5OSI </li> 47 * <li> NewMarkAlphaOSI </li> 48 * <li> ZeroOrderHoldOSI </li> 49 * </ul> 50 * 51 */ 52 class OneStepIntegrator :public std::enable_shared_from_this<OneStepIntegrator> 53 { 54 55 protected: 56 /* serialization hooks */ 57 ACCEPT_SERIALIZATION(OneStepIntegrator); 58 59 /** type/name of the Integrator */ 60 OSI::TYPES _integratorType; 61 62 /** a graph of dynamical systems to integrate 63 * For the moment, we point to the graph of dynamical systems in 64 * the topology. We use the properties "osi" to check if the dynamical 65 * system is integrated by this osi. It has to be improved by using a subgraph 66 * to avoid the use of checkOSI 67 */ 68 SP::DynamicalSystemsGraph _dynamicalSystemsGraph; 69 70 /** size of the memory for the integrator */ 71 unsigned int _sizeMem; 72 73 /** steps of the integrator */ 74 unsigned int _steps; 75 76 /** _levelMinForOutput is the minimum level for the output 77 * needed by the OneStepIntegrator 78 */ 79 unsigned int _levelMinForOutput; 80 81 /** _levelMaxForOutput is the maximum level for the output 82 * needed by the OneStepIntegrator 83 */ 84 unsigned int _levelMaxForOutput; 85 86 /** _levelMinForInput is the minimum level for the input 87 * needed by the OneStepIntegrator 88 */ 89 unsigned int _levelMinForInput; 90 91 /** _levelMaxForInput is the maximum level for the input 92 * needed by the OneStepIntegrator 93 */ 94 unsigned int _levelMaxForInput; 95 96 bool _isInitialized; 97 98 /** boolean variable to force an explicit evaluation of the Jacobians 99 * mapping of relations only at the beginning of the time--step and 100 * not in the Newton iteration 101 */ 102 103 bool _explicitJacobiansOfRelation; 104 105 106 107 /** A link to the simulation that owns this OSI */ 108 SP::Simulation _simulation; 109 110 /** basic constructor with OSI Id 111 * \param type integrator type/name 112 */ OneStepIntegrator(const OSI::TYPES & type)113 OneStepIntegrator(const OSI::TYPES& type) 114 : _integratorType(type), _sizeMem(1), _steps(0), 115 _levelMinForOutput(0), _levelMaxForOutput(0), 116 _levelMinForInput(0), _levelMaxForInput(0), 117 _isInitialized(false), _explicitJacobiansOfRelation(false) {}; 118 119 /** struct to add terms in the integration. Useful for Control */ 120 SP::ExtraAdditionalTerms _extraAdditionalTerms; 121 122 /** Compare interaction and current OSI levels for input and output. 123 Reset interaction if they are not compliant. 124 \param inter a reference to an Interaction 125 */ 126 void _check_and_update_interaction_levels(Interaction& inter); 127 128 /** initialization of the work vectors and matrices (properties) related to 129 * one dynamical system on the graph and needed by the osi -- common code. 130 * \param ds the dynamical system 131 */ 132 SP::VectorOfVectors _initializeDSWorkVectors(SP::DynamicalSystem ds); 133 134 /** default constructor */ OneStepIntegrator()135 OneStepIntegrator() {}; 136 137 private: 138 139 /** copy constructor, private, no copy nor pass-by value allowed */ 140 OneStepIntegrator(const OneStepIntegrator&); 141 142 /** assignment (private => forbidden) 143 * \param OSI 144 * \return OneStepIntegrator& 145 */ 146 OneStepIntegrator& operator=(const OneStepIntegrator& OSI); 147 148 149 public: 150 151 /** destructor 152 */ ~OneStepIntegrator()153 virtual ~OneStepIntegrator() {}; 154 155 /*! @name Attributes access 156 @{ */ 157 158 /** id of the integrator (see list in OSI::TYPES enum) 159 \return int 160 */ getType() const161 inline OSI::TYPES getType() const 162 { 163 return _integratorType; 164 } 165 166 /** get the graph of dynamical systems associated with the Integrator 167 warning: returns the whole ds graph, not only ds integrated by the present osi. 168 \return a SP::DynamicalSystemsGraph 169 */ dynamicalSystemsGraph() const170 inline SP::DynamicalSystemsGraph dynamicalSystemsGraph() const 171 { 172 return _dynamicalSystemsGraph; 173 }; 174 175 /** set the graph of dynamical systems associated with the Integrator 176 */ setDynamicalSystemsGraph(SP::DynamicalSystemsGraph dsg)177 inline void setDynamicalSystemsGraph(SP::DynamicalSystemsGraph dsg) 178 { 179 _dynamicalSystemsGraph = dsg; 180 }; 181 182 /** get number of internal memory vectors needed in dynamical systems integrated with this osi. 183 * \return an unsigned int 184 */ getSizeMem() const185 inline unsigned int getSizeMem() const 186 { 187 return _sizeMem; 188 }; 189 190 /** get the Simulation that owns the OneStepIntegrator (pointer link) 191 * \return a pointer to Simulation 192 */ simulation() const193 inline SP::Simulation simulation() const 194 { 195 return _simulation; 196 } 197 198 /** set the Simulation of the OneStepIntegrator 199 * \param newS a pointer to Simulation 200 */ setSimulationPtr(SP::Simulation newS)201 inline void setSimulationPtr(SP::Simulation newS) 202 { 203 _simulation = newS; 204 } 205 206 /** minimal level required for output var used with this integration scheme. 207 var[level] is the derivative of order 'level' of var. 208 */ levelMinForOutput()209 virtual unsigned int levelMinForOutput() 210 { 211 return _levelMinForOutput; 212 } 213 214 /** maximal level required for output var used with this integration scheme. 215 var[level] is the derivative of order 'level' of var. 216 */ levelMaxForOutput()217 virtual unsigned int levelMaxForOutput() 218 { 219 return _levelMaxForOutput; 220 } 221 222 /** minimal level required for input var used with this integration scheme. 223 var[level] is the derivative of order 'level' of var. 224 */ levelMinForInput()225 virtual unsigned int levelMinForInput() 226 { 227 return _levelMinForInput; 228 } 229 230 /** maximal level required for input var used with this integration scheme. 231 var[level] is the derivative of order 'level' of var. 232 */ levelMaxForInput()233 virtual unsigned int levelMaxForInput() 234 { 235 return _levelMaxForInput; 236 } 237 238 /** get the number of index sets required for the simulation 239 * \return unsigned int 240 */ 241 virtual unsigned int numberOfIndexSets() const = 0; 242 243 /** @} end of members access group. */ 244 245 /*! @name internal memory (graph properties) management 246 @{ */ 247 248 isInitialized()249 inline bool isInitialized(){return _isInitialized;}; 250 setIsInitialized(bool value)251 inline void setIsInitialized( bool value) {_isInitialized = value;}; 252 explicitJacobiansOfRelation()253 bool explicitJacobiansOfRelation() 254 { 255 return _explicitJacobiansOfRelation; 256 } 257 setExplicitJacobiansOfRelation(bool newval)258 void setExplicitJacobiansOfRelation(bool newval) 259 { 260 _explicitJacobiansOfRelation = newval; 261 }; 262 263 /** initialise the integrator 264 */ 265 virtual void initialize(); 266 267 /** Initialization process of the nonsmooth problems 268 linked to this OSI*/ initialize_nonsmooth_problems()269 virtual void initialize_nonsmooth_problems(){}; 270 271 /** initialization of the work vectors and matrices (properties) related to 272 * one dynamical system on the graph and needed by the osi 273 * \param t time of initialization 274 * \param ds the dynamical system 275 */ 276 virtual void initializeWorkVectorsForDS(double t, SP::DynamicalSystem ds) = 0 ; 277 278 /** initialization of the work vectors and matrices (properties) related to 279 * one interaction on the graph and needed by the osi 280 * \param inter the interaction 281 * \param interProp the properties on the graph 282 * \param DSG the dynamical systems graph 283 */ 284 virtual void initializeWorkVectorsForInteraction(Interaction &inter, 285 InteractionProperties& interProp, 286 DynamicalSystemsGraph & DSG) = 0 ; 287 288 289 /** @} end of internal memory management functions */ 290 291 /*! @name Computation functions 292 @{ */ 293 294 /** compute interaction output (y) for all levels and swaps in memory 295 \param inter the interaction to update 296 \param time value for output computation 297 \param interaction_properties properties of the interaction, in the Interaction Graph I0 298 */ 299 void update_interaction_output(Interaction& inter, double time, InteractionProperties& interaction_properties); 300 301 /** compute the initial state (for dynamical system variables) of the Newton loop. */ computeInitialNewtonState()302 virtual void computeInitialNewtonState(){ 303 // Default behavior : does nothing and used the current state as starting state of the Newton iteration 304 } 305 306 /** return the maximum of all norms for the discretized residus of DS 307 * \return a double 308 */ computeResidu()309 virtual double computeResidu(){ 310 // default : error 311 THROW_EXCEPTION("OneStepIntegrator::computeResidu not implemented for integrator of type " + std::to_string(_integratorType)); 312 return 0.0; 313 } 314 315 /** integrates the Dynamical System linked to this integrator, without taking constraints 316 into account. 317 */ computeFreeState()318 virtual void computeFreeState() 319 { 320 // default : error 321 THROW_EXCEPTION("OneStepIntegrator::computeFreeState not implemented for integrator of type " + std::to_string(_integratorType)); 322 } 323 324 /** integrates the Interaction linked to this integrator, without taking non-smooth effects into account 325 * \param vertex_inter of the interaction graph 326 * \param osnsp pointer to OneStepNSProblem 327 */ computeFreeOutput(InteractionsGraph::VDescriptor & vertex_inter,OneStepNSProblem * osnsp)328 virtual void computeFreeOutput(InteractionsGraph::VDescriptor& vertex_inter, OneStepNSProblem* osnsp) 329 { 330 // default : error 331 THROW_EXCEPTION("OneStepIntegrator::computeFreeOutput not implemented for integrator of type " + std::to_string(_integratorType)); 332 } 333 334 /** compute the residu of the output of the relation (y) 335 * This computation depends on the type of OSI 336 * \param time time of computation 337 * \param indexSet the index set of the interaction that are concerned 338 */ 339 virtual double computeResiduOutput(double time, SP::InteractionsGraph indexSet); 340 /** compute the residu of the input of the relation (R or p) 341 * This computation depends on the type of OSI 342 * \param time time of computation 343 * \param indexSet the index set of the interaction that are concerned 344 */ 345 virtual double computeResiduInput(double time, SP::InteractionsGraph indexSet); 346 347 /** integrate the system, between tinit and tend, with possible stop at tout 348 * \param tinit start time 349 * \param tend expected end time 350 * \param tout real end time 351 * \param idid extra flag, meaningful only for OSI used in EventDriven schemes 352 */ 353 virtual void integrate(double& tinit, double& tend, double& tout, int& idid) = 0; 354 355 /** set to zero all the r vectors of the DynamicalSystems integrated by this OSI 356 */ 357 void resetAllNonSmoothParts(); 358 359 /** set to zero all the r vectors of the DynamicalSystems of the present OSI for a given level 360 * \param level 361 */ 362 void resetNonSmoothPart(unsigned int level); 363 364 /** update the state of the DynamicalSystem attached to this Integrator 365 * \param level level of interest for the dynamics 366 * level is set to 0 by default since in all time-stepping schemes we update all the state 367 * whatever the value of level 368 */ 369 virtual void updateState(const unsigned int level) = 0; 370 371 /** update the state of the DynamicalSystem attached to this Integrator 372 * level is set to 0 by default since in all time-stepping schemes we update all the state 373 * whatever the value of level 374 */ updateState()375 void updateState() { updateState(0); } 376 377 /** update the output of the Interaction attached to this Integrator 378 */ 379 virtual void updateOutput(double time); 380 381 /** update the input of the Interaction attached to this Integrator 382 */ 383 virtual void updateInput(double time); 384 385 /** update the output of the Interaction attached to this Integrator 386 * \param time current time 387 * \param level level of interest for the dynamics 388 */ 389 virtual void updateOutput(double time, unsigned int level); 390 391 /** update the input of the Interaction attached to this Integrator 392 * \param time current time 393 * \param level level of interest for the dynamics 394 */ 395 virtual void updateInput(double time, unsigned int level); 396 397 /** */ 398 virtual void prepareNewtonIteration(double time) = 0; 399 /** @} end of computation functions */ 400 401 /*! @name Misc 402 @{ */ 403 404 /** print the data to the screen 405 */ 406 virtual void display() = 0; 407 408 /** Apply the rule to one Interaction to known if is it should be included 409 * in the IndexSet of level i 410 * \param inter 411 * \param i 412 * \return bool 413 */ addInteractionInIndexSet(SP::Interaction inter,unsigned int i)414 virtual bool addInteractionInIndexSet(SP::Interaction inter, unsigned int i) 415 { 416 THROW_EXCEPTION("OneStepIntegrator::addInteractionInIndexSet - Should be called at this level"); 417 return 0; 418 } 419 ; 420 421 /** Apply the rule to one Interaction to know if is it should be removed 422 * from the IndexSet of level i 423 * \param inter 424 * \param i 425 * \return bool 426 */ removeInteractionFromIndexSet(SP::Interaction inter,unsigned int i)427 virtual bool removeInteractionFromIndexSet(SP::Interaction inter, unsigned int i) 428 { 429 THROW_EXCEPTION("OneStepIntegrator::removeInteractionFromIndexSet - Should not be called at this level"); 430 return 0; 431 }; 432 433 /** get the ExtraAdditionalTerms. 434 * \return the ExtraAdditionalTerms 435 */ extraAdditionalTerms()436 inline SP::ExtraAdditionalTerms extraAdditionalTerms() 437 { 438 return _extraAdditionalTerms; 439 } 440 441 /** set the ExtraAdditionalTerms to add smooth terms for the integration process. 442 * Useful when a control loop is added to a DynamicalSystem. 443 * \param eat the ExtraAdditionalTerms to use 444 */ setExtraAdditionalTerms(SP::ExtraAdditionalTerms eat)445 inline void setExtraAdditionalTerms(SP::ExtraAdditionalTerms eat) 446 { 447 _extraAdditionalTerms = eat; 448 } 449 /** True if the dynamical system (a vertex in the ds graph) is integrated by this osi. 450 \param dsi the iterator on the node of the graph corresponding to the dynamical system of interest. 451 */ checkOSI(DynamicalSystemsGraph::VIterator dsi)452 inline bool checkOSI(DynamicalSystemsGraph::VIterator dsi) 453 { 454 return (_dynamicalSystemsGraph->properties(*dsi).osi.get()) == this; 455 }; 456 457 /** True if the dynamical system (a vertex in the ds graph) is integrated by this osi. 458 \param dsgv the descriptor of the node in the graph corresponding to the dynamical system of interest. 459 */ checkOSI(DynamicalSystemsGraph::VDescriptor dsgv)460 inline bool checkOSI(DynamicalSystemsGraph::VDescriptor dsgv) 461 { 462 return (_dynamicalSystemsGraph->properties(dsgv).osi.get()) == this; 463 }; 464 465 /** @} end of misc */ 466 467 /* visitors hook */ 468 VIRTUAL_ACCEPT_VISITORS(OneStepIntegrator); 469 470 }; 471 472 #endif // ONESTEPINTEGRATOR_H 473