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 MultipleImpact.hpp 19 * \brief Linear Complementarity Problem formulation and solving 20 */ 21 22 #ifndef _OSNSMULTIPLEIMPACT_ 23 #define _OSNSMULTIPLEIMPACT_ 24 25 #include "LinearOSNS.hpp" 26 #include "SiconosConst.hpp" // for MACHINE_PREC 27 #include <string> 28 29 #define DEFAULT__tolImpact MACHINE_PREC 30 #define DEFAULT_TOL_VEL MACHINE_PREC 31 #define DEFAULT_TOL_ENER MACHINE_PREC 32 33 /** Formalization and Resolution of a Multiple Impact Non-Smooth problem. 34 35 \todo write a short introduction about MultipleImpact ... 36 */ 37 class MultipleImpact : public LinearOSNS 38 { 39 private: 40 /** serialization hooks 41 */ 42 ACCEPT_SERIALIZATION(MultipleImpact); 43 44 //! Time-like variable (Impulse) 45 double _impulseVariable = 0.; 46 //! Time variable 47 double _timeVariable = 0.; 48 //! Number of contacts (only the active contacts) 49 unsigned int _nContact = 0; 50 //! Maximal number of steps for each computation 51 unsigned int _nStepMax = 100000; 52 //! Tolerance to define zero 53 double _tolImpact = DEFAULT__tolImpact; 54 //! Type of the compliance model 55 std::string _typeCompLaw = "BiStiffness"; 56 //Velocity of bodies during impact 57 //SP::SiconosVector VelAllBody; 58 // Relative velocity at all Interactions (with or without contact) 59 //SP::SiconosVector VelAllIteractions; 60 //! Relative velocity during impact (at the end of each calculation step) 61 SP::SiconosVector _velocityContact; 62 //! Relative velocity during impact (at the beginning of each calculation step) 63 SP::SiconosVector _oldVelocityContact; 64 //! Potential energy during impact (at the end of each calculation step) 65 SP::SiconosVector _energyContact; 66 //! Work done during the last compression phase at contact 67 SP::SiconosVector _WorkcContact; 68 //! Distribution vector to distribute the incremental impulse at contact 69 SP::SiconosVector _distributionVector; 70 /** State of contacts at the beginning of impact 71 if *_stateContact[i] = 0 => no impact at this contact (at contact with positive relative velocity and no potential energy, may be the impact has been terminated at this contact) 72 if *_stateContact[i] = 1 => impact takes place at this contact without potential energy (beginning of impact or repeating impact) 73 if *_stateContact[i] = 2 => impact takes place with not-zero potential energy */ 74 SP::IndexInt _stateContact; 75 //!Stiffness at contacts 76 SP::SiconosVector _Kcontact; 77 //! Restitution coefficient of contacts 78 SP::SiconosVector _restitutionContact; 79 //! Elasticity coefficient of contacts 80 SP::SiconosVector _elasticyCoefficientcontact; 81 //! Incremental impulse at contacts 82 SP::SiconosVector _deltaImpulseContact; 83 //! Total impulse at contacts 84 SP::SiconosVector _tolImpulseContact; 85 //! Impulse at contacts for each update time 86 SP::SiconosVector _impulseContactUpdate; 87 //! Force at contacts 88 SP::SiconosVector _forceContact; 89 //! ID of the primary contact 90 unsigned int _primaryContactId = 0; 91 /** Indicator about the selection of the primary contact 92 true if primary contact is selected according to the potential energy 93 false if primary contact is selected according to the relative velocity */ 94 bool _isPrimaryContactEnergy = false; 95 //! Relative velocity at primary contact 96 double _relativeVelocityPrimaryContact = 0.; 97 //! Potential energy at primary contact 98 double _energyPrimaryContact = 0.; 99 //! Step size for the iterative calculation 100 double _deltaP = 0.; 101 //! Name of file into which the datat is writen 102 std::string _namefile = "DataMultipleImpact.dat"; 103 /** YesWriteData = true ==>save the data during impact 104 YesWriteData = false ==> not save the data during impact */ 105 bool _saveData = false; 106 /** bool variable to set the step size for multiple impact computation 107 If IsNumberOfStepsEst = true ==> estimate the step size from the state of the dynamic system before impact and the number of step needed 108 Number of steps after which the data is saved */ 109 unsigned int _nStepSave = 100; //! If IsNumberOfStepsEst = false ==> user choose the step size 110 //! Matrix on which the data during impact is saved 111 SP::SiconosMatrix _DataMatrix; 112 //! Number of points to be save during impacts 113 unsigned int _sizeDataSave = 1000; 114 /** indicator on the termination of the multiple impact process 115 _IsImpactEnd = true: impact is terminated 116 _IsImpactEnd = false: otherwise */ 117 bool _IsImpactEnd = true; 118 //! Tolerance to define a negligeble value for a velocity grandeur 119 double _Tol_Vel = DEFAULT_TOL_VEL; 120 //! Tolerance to define a negligeable value for a potential energy grandeur 121 double _Tol_Ener = DEFAULT_TOL_ENER; 122 //! Epsilon to define a zero value for relative velocity in termination condition 123 double _ZeroVel_EndIm = DEFAULT_TOL_VEL; 124 //! Epsilon to define a zero value for potential energy in termination condition 125 double _ZeroEner_EndIm = DEFAULT_TOL_ENER; 126 //! we start to save data from _stepMinSave to _stepMaxSave 127 unsigned int _stepMinSave = 1, _stepMaxSave = _nStepMax; 128 public: 129 130 /** default constructor 131 */ MultipleImpact()132 MultipleImpact():LinearOSNS(){}; 133 134 /** Constructor from data (step size is required here) 135 * \param type the type of the compliance law 136 * \param step step size estimated 137 */ 138 MultipleImpact(std::string type, double step = 1.0e-5); 139 140 //!Destructor ~MultipleImpact()141 ~MultipleImpact(){}; 142 143 /* To get the type of the compliance law at contact 144 * \return std::string 145 */ get_typeCompLaw() const146 inline std::string get_typeCompLaw() const 147 { 148 return _typeCompLaw; 149 }; 150 151 /** To set the type of the compliance law 152 * \param newTypeLaw 153 */ 154 void set_typeCompLaw(std::string newTypeLaw); 155 156 /** To set the tolerance to define zero 157 * \param newTolZero 158 */ 159 void setTolImpact(double newTolZero); 160 161 /** To get the tolerance to define zero 162 * \return double 163 */ getTolImpact()164 inline double getTolImpact() 165 { 166 return _tolImpact; 167 }; 168 169 /** To set the flag to save the data during impact or not 170 * \param var 171 */ 172 void SetSaveData(bool var); 173 174 /** To set the name for the output file 175 * \param file_name 176 */ 177 void SetNameOutput(std::string file_name); 178 179 /** To get step size 180 * \return double 181 */ GetStepSize()182 inline double GetStepSize() 183 { 184 return _deltaP; 185 }; 186 187 /* To get the duration of multiple impacts process 188 * \return double 189 */ DurationImpact()190 inline double DurationImpact() 191 { 192 return _timeVariable; 193 }; 194 195 /** To set the variable _nStepSave 196 * \param var 197 */ 198 void SetNstepSave(unsigned int var); 199 200 /** To set the maximal number of steps allowed for each computation 201 * \param var 202 */ 203 void SetNstepMax(unsigned int var); 204 205 /** Set number of points to be saved during impact 206 * \param var 207 */ 208 void SetSizeDataSave(unsigned int var); 209 210 /** Set tolerence to define whether or not a velocity is zero 211 * \param var 212 */ 213 void SetTolVel(double var); 214 215 /** Set tolerence to define whether or not a potential energy is zero 216 * \param var 217 */ 218 void SetTolEner(double var); 219 220 /** Set epsilon _ZeroVel_EndIm 221 * \param var 222 */ 223 void SetZeroVelEndImp(double var); 224 225 /** Set epsilon _ZeroEner_EndIm 226 * \param var 227 */ 228 void SetZeroEnerEndImp(double var); 229 230 /** Set the step number to start the data save and step number to stop save 231 * \param min 232 * \param max 233 */ 234 void SetStepMinMaxSave(unsigned int min, unsigned int max); 235 236 /** To compare a double number with zero 237 * \param var 238 * \return bool 239 */ 240 bool isZero(double var); 241 242 /** To compare a velocity value with zero 243 * \param var 244 * \return bool 245 */ 246 bool isVelNegative(double var); 247 248 /** To compare an energy value with zero 249 * \param var 250 * \return bool 251 */ 252 bool isEnerZero(double var); 253 254 /** To select the pramary contact 255 */ 256 void SelectPrimaContact(); 257 258 /** Calculate the vector of distributing rule */ 259 void Compute_distributionVector(); 260 261 /** Compute the normal imulse at contacts */ 262 void ComputeImpulseContact(); 263 264 /** Compute the relative velocity at contacts */ 265 void Compute_velocityContact(); 266 267 /** Compute the potential energy at contacts during each computation step */ 268 void Compute_energyContact(); 269 270 /** Compute the velocity of the bodies during impact */ 271 void UpdateDuringImpact(); 272 273 /** Run the iterative procedure to solve the multiple impact problem */ 274 void ComputeImpact(); 275 276 /** Post-compute for multiple impacts */ 277 void PostComputeImpact(); 278 279 /** Check if the multiple impacts process is terminated or not 280 * \return bool 281 */ 282 bool IsMulImpactTerminate(); 283 284 /** To allocate the memory */ 285 void AllocateMemory(); 286 287 /** To build the vector of stiffnesses and restitution coefficient at contacts */ 288 void BuildParaContact(); 289 290 /** To get the velocity of bodies, relative velocity and potential energy at the beginning of impact */ 291 void InitializeInput(); 292 293 /** To check the state of contacts during impact */ 294 void Check_stateContact(); 295 296 /** Pre-compute for multiple impacs */ 297 void PreComputeImpact(); 298 299 /** To get the primary contact according to the relative velocity 300 In this case, the primary contact correspond to the contact at which the relative velocity 301 is minimum (the relative velocity for two approching bodies is negative so the magnitude of 302 the relative velocity at the primary contact is maximum) 303 */ 304 void PrimConVelocity(); 305 306 /** To get the primary contact according to the potential energy. In this case, the primary 307 contact corresponds to the one at which the potential energy is maximum 308 */ 309 void PrimConEnergy(); 310 311 /** To decide if the primary contact is selected according to the relative velocity or to the 312 * potential energy. The first case happens when there is no potential energy at any contact 313 * \return bool 314 */ 315 bool IsEnermaxZero(); 316 317 /** Verify if the minimum relative velocity at contacts is negative or not 318 * \return bool 319 */ 320 bool IsVcminNegative(); 321 322 /** compute the unknown post-impact relative velocity and post-impact impulse 323 * \param time 324 * \return int 325 */ 326 int compute(double time); 327 328 /**initialize 329 * \param sim 330 */ 331 void initialize(SP::Simulation sim); 332 333 334 bool checkCompatibleNSLaw(NonSmoothLaw& nslaw); 335 /** print the data to the screen */ 336 void display() const; 337 338 /** To write a SiconosVector into a matrix 339 * \param v 340 * \param row position starting to write 341 * \param col position starting to write 342 */ 343 344 void WriteVectorIntoMatrix(const SiconosVector& v, const unsigned int row, const unsigned int col); 345 346 /** Save data for each step 347 * \param i pointer to be save */ 348 void SaveDataOneStep(unsigned int i); 349 350 /** Estimate size of data matrix 351 * \return unsigned int 352 */ 353 unsigned int EstimateNdataCols(); 354 355 ACCEPT_STD_VISITORS(); 356 }; 357 358 #endif 359