ClumpParticle.h
Go to the documentation of this file.
1 //Copyright (c) 2013-2023, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #ifndef MultiParticle_H
27 #define MultiParticle_H
28 
29 #include "Mercury3D.h"
30 #include "Math/SmallVector.h"
31 #include "BaseParticle.h"
32 #include "BaseInteractable.h"
33 #include "NonSphericalParticle.h"
35 #include "DPMBase.h"
41 {
42 public:
46  ClumpParticle();
47 
51  ClumpParticle(const ClumpParticle& p);
52 
58  ~ClumpParticle() override;
59 
60  ClumpParticle* copy() const override;
61 
62 
63  void write(std::ostream& os) const override
64  {
66  os << " DzhanibekovParticle " << DzhanibekovParticle_;
67  os << " VerticallyOriented " << VerticallyOriented_;
68 
69  }
70 
71  void read(std::istream& is) override;
72 
73  std::string getName() const override;
74 
75  bool isSphericalParticle() const override
76  {
77  return true;
78  }
79 
80  void computeMass(const ParticleSpecies& s) override;
81 
82  void setInitInertia(MatrixSymmetric3D inertia);
83 
84  void rotateTensorOfInertia();
85 
86  void addPebble(Vec3D position, Mdouble radius);
87 
88  void setClump();
89 
90  void setPrincipalDirections(Matrix3D directions);
91 
92  void setInitPrincipalDirections(Matrix3D directions);
93 
95  {
96 
98  }
100  {
102  }
104  {
106  }
107 
108  // Methods to obtain initial principal directions
110  {
112  }
114  {
116  }
118  {
120  }
121 
122 
123  int NPebble() const; // Number of pebbles (for a clump particle)
124 
125  void actionsAfterAddObject() override; // The function that updates clump quantities after adding a pebble
126 
127  void updatePebblesVelPos();
128 
129  void integrateBeforeForceComputation(double time, double timeStep) override;
130 
131  void integrateAfterForceComputation(double time, double timeStep) override;
132 
133  void angularAccelerateClumpIterative(double timeStep); // Specific method of time integration for a clump particle
134 
135  void rotatePrincipalDirections(Vec3D rotation);
136 
137  // Principal direction
139  {
143  }
145  {
149  }
151  {
155  }
156 
157  std::vector<Mdouble> getPebbleRadius() const {
158  return pebbleRadius_;
159  }
160 
161  // add pointer to pebble pointers list
162  void setPebble(int kPebble, ClumpParticle* pPebble) {
163  pebbleParticles_[kPebble] = pPebble;
164  }
165 
166  // Sets the particle to be a pebble of a given clump
167  void setClump(ClumpParticle* master) {
168  isClump_ = false;
169  isPebble_ = true;
170  clumpParticle = master;
171  }
172 
174  {
175  clumpMass_ = mass;
176  invMass_ = 1 / clumpMass_;
177  }
178 
179  // Extra viscous damping on a clump
180  void setDamping(Mdouble damp)
181  {
182  viscousDamping_ = damp;
183  }
184 
185  Mdouble getKineticEnergy() const override{
186  Mdouble res = 0;
187  if (isClump_) {
188  Vec3D v = getVelocity();
189  res = 0.5 * clumpMass_ * (v.X * v.X + v.Y * v.Y + v.Z * v.Z );
190  }
191  return res;
192  }
193 
194  Mdouble getRotationalEnergy() const override{
195  Mdouble res = 0;
196  if (isClump_) {
197  Vec3D nn = getAngularVelocity();
198  Mdouble nl = nn.getLength();
199  Mdouble tol = 1e-10;
200  if (nl > tol) {
201  nn /= nl;
203  Mdouble II = Vec3D::dot(nn, (getInertia() * nn));
204  res = 0.5 * II * ww;
205  }
206  }
207  return res;
208 
209  }
210 
211  std::vector <Vec3D> getPebblePositions(){
212  std::vector <Vec3D> globalPos;
216  for (int i = 1; i <= NPebble(); i++){
217  globalPos.push_back(getPosition() + e1 * pebblePos_[i - 1].X + e2 * pebblePos_[i - 1].Y + e3 * pebblePos_[i - 1].Z);
218  }
219  return globalPos;
220  }
221 
222  std::vector <Mdouble> getPebbleRadii(){ return pebbleRadius_; }
223 
224  // Methods setting and getting some extra Boolean properties
225 
226  // check if particle is in "Dzhanibekov" state
228  {
229  return DzhanibekovParticle_;
230  }
231 
232  // check if particle is "vertically oriented"
234  {
235  return VerticallyOriented_;
236  }
237 
238  // set the "Dzhanibekov" state
240  {
242  }
243 
244  // set "vertically oriented" state
245  void setVerticallyOriented( bool d)
246  {
248  }
249 
250  unsigned getNumberOfFieldsVTK() const override
251  {
252  return 2;
253  }
254 
255  std::string getTypeVTK(unsigned i) const override
256  {
257  return "Int8";
258  }
259 
260 
261  std::string getNameVTK(unsigned i) const override
262  {
263  if (i==0)
264  return "DzhanibekovParticle";
265  else
266  return "VerticallyOriented";
267  }
268 
269 
270  std::vector<Mdouble> getFieldVTK(unsigned i) const override
271  {
272  if (i==0)
273  return std::vector<Mdouble>(1, DzhanibekovParticle_);
274  else
275  return std::vector<Mdouble>(1, VerticallyOriented_);
276  }
277 
278  void updateExtraQuantities();
279 
280 
281 private:
282 
283  int nPebble_; // Number of pebbles
284 
285  bool DzhanibekovParticle_; // This property is needed to quantify Dzhanibekov gas properties
286 
287  bool VerticallyOriented_; // This property is useful for mechnical stability simulations (Gomboc, Dominos)
288 
289  Vec3D angularAcceleration_; // Clump angular acceleration
290 
291  Mdouble clumpMass_; // Clump mass
292 
293  Mdouble viscousDamping_; // Viscous damping parameter - for extra clump damping
294 
295  MatrixSymmetric3D clumpInertia_; // Clump tensor of inertia
296 
297  MatrixSymmetric3D clumpInitInertia_; // Clump initial tensor of inertia
298 
299  std::vector<Vec3D> pebblePos_; // positions of pebbles in a local coordinate system
300 
301  std::vector<Mdouble> pebbleRadius_; // radii of pebbles
302 
303  Matrix3D principalDirections_; // clump's principal directions
304 
305  Matrix3D initPrincipalDirections_; // clump's initial principal directions
306 
307  std::vector<ClumpParticle*> pebbleParticles_; // pointers to pebbles
308 
309  //Helper functions
310 
311  // Converts a Matrix3D into a MatrixSymmetric3D
312  MatrixSymmetric3D MtoS( Matrix3D M){ return MatrixSymmetric3D(M.XX, M.XY, M.XZ, M.YY, M.YZ, M.ZZ);}
313 
314  // Converts a MatrixSymmetric3D into a Matrix3D
315  Matrix3D StoM( MatrixSymmetric3D M){ return Matrix3D(M.XX, M.XY, M.XZ, M.XY, M.YY, M.YZ, M.XZ, M.YZ, M.ZZ);}
316 
317  // Transposes the matrix
318  Matrix3D transpose(Matrix3D M){ return Matrix3D(M.XX, M.YX, M.ZX, M.XY, M.YY, M.ZY, M.XZ, M.YZ, M.ZZ);}
319 };
320 
321 #endif
double Mdouble
Definition: GeneralDefine.h:34
virtual const Vec3D & getAngularVelocity() const
Returns the angular velocity of this interactable.
Definition: BaseInteractable.cc:341
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
Definition: BaseInteractable.cc:329
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Definition: BaseInteractable.h:218
void write(std::ostream &os) const override
Particle print function, which accepts an std::ostream as input.
Definition: BaseParticle.cc:336
bool isClump_
The particle is pebble.
Definition: BaseParticle.h:748
MatrixSymmetric3D getInertia() const
Definition: BaseParticle.h:331
BaseParticle * clumpParticle
Function that updates necessary quantities of a clump particle after adding a pebble.
Definition: BaseParticle.h:744
bool isPebble_
pointer to a clump particle (for a pebble)
Definition: BaseParticle.h:746
Mdouble invMass_
Particle radius_.
Definition: BaseParticle.h:685
Definition: ClumpParticle.h:41
Mdouble getKineticEnergy() const override
Definition: ClumpParticle.h:185
void setPrincipalDirections_e1(Vec3D e)
Definition: ClumpParticle.h:138
int nPebble_
Definition: ClumpParticle.h:283
Vec3D getInitPrincipalDirections_e1() const
Definition: ClumpParticle.h:109
void angularAccelerateClumpIterative(double timeStep)
Definition: ClumpParticle.cc:358
void updatePebblesVelPos()
Definition: ClumpParticle.cc:251
ClumpParticle * copy() const override
Definition: ClumpParticle.cc:97
bool isSphericalParticle() const override
Definition: ClumpParticle.h:75
void setPebble(int kPebble, ClumpParticle *pPebble)
Definition: ClumpParticle.h:162
MatrixSymmetric3D clumpInitInertia_
Definition: ClumpParticle.h:297
bool getDzhanibekovParticle()
Definition: ClumpParticle.h:227
bool DzhanibekovParticle_
Definition: ClumpParticle.h:285
Mdouble clumpMass_
Definition: ClumpParticle.h:291
std::vector< Mdouble > getPebbleRadii()
Definition: ClumpParticle.h:222
ClumpParticle()
Basic Particle constructor, creates a particle at (0,0,0) with radius, mass and inertia equal to 1.
Definition: ClumpParticle.cc:37
Vec3D getPrincipalDirections_e2() const
Definition: ClumpParticle.h:99
MatrixSymmetric3D MtoS(Matrix3D M)
Definition: ClumpParticle.h:312
MatrixSymmetric3D clumpInertia_
Definition: ClumpParticle.h:295
void integrateAfterForceComputation(double time, double timeStep) override
Definition: ClumpParticle.cc:334
bool getVerticallyOriented()
Definition: ClumpParticle.h:233
std::vector< Mdouble > pebbleRadius_
Definition: ClumpParticle.h:301
Vec3D getInitPrincipalDirections_e2() const
Definition: ClumpParticle.h:113
bool VerticallyOriented_
Definition: ClumpParticle.h:287
std::vector< Vec3D > pebblePos_
Definition: ClumpParticle.h:299
void rotatePrincipalDirections(Vec3D rotation)
Definition: ClumpParticle.cc:151
std::string getName() const override
Definition: ClumpParticle.cc:103
void write(std::ostream &os) const override
Particle print function, which accepts an std::ostream as input.
Definition: ClumpParticle.h:63
Mdouble viscousDamping_
Definition: ClumpParticle.h:293
void setClumpMass(Mdouble mass)
Definition: ClumpParticle.h:173
std::vector< Mdouble > getFieldVTK(unsigned i) const override
Definition: ClumpParticle.h:270
void setInitPrincipalDirections(Matrix3D directions)
Definition: ClumpParticle.cc:145
void computeMass(const ParticleSpecies &s) override
Computes the particle's (inverse) mass and inertia.
Definition: ClumpParticle.cc:391
int NPebble() const
Definition: ClumpParticle.cc:115
void setClump()
Definition: ClumpParticle.cc:120
void actionsAfterAddObject() override
Definition: ClumpParticle.cc:192
void rotateTensorOfInertia()
Definition: ClumpParticle.cc:224
void setInitInertia(MatrixSymmetric3D inertia)
Definition: ClumpParticle.cc:216
Vec3D angularAcceleration_
Definition: ClumpParticle.h:289
std::string getNameVTK(unsigned i) const override
Definition: ClumpParticle.h:261
Vec3D getPrincipalDirections_e3() const
Definition: ClumpParticle.h:103
void setVerticallyOriented(bool d)
Definition: ClumpParticle.h:245
void updateExtraQuantities()
Definition: ClumpParticle.cc:405
Mdouble getRotationalEnergy() const override
Calculates the particle's rotational kinetic energy.
Definition: ClumpParticle.h:194
void integrateBeforeForceComputation(double time, double timeStep) override
Definition: ClumpParticle.cc:287
void setDzhanibekovParticle(bool d)
Definition: ClumpParticle.h:239
~ClumpParticle() override
Destructor, needs to be implemented and checked to see if it is the largest or smallest particle curr...
void setDamping(Mdouble damp)
Definition: ClumpParticle.h:180
Matrix3D StoM(MatrixSymmetric3D M)
Definition: ClumpParticle.h:315
void setPrincipalDirections(Matrix3D directions)
Definition: ClumpParticle.cc:139
void setPrincipalDirections_e3(Vec3D e)
Definition: ClumpParticle.h:150
std::vector< Vec3D > getPebblePositions()
Definition: ClumpParticle.h:211
Matrix3D initPrincipalDirections_
Definition: ClumpParticle.h:305
void setPrincipalDirections_e2(Vec3D e)
Definition: ClumpParticle.h:144
void read(std::istream &is) override
Particle read function, which accepts an std::istream as input.
Definition: ClumpParticle.cc:108
unsigned getNumberOfFieldsVTK() const override
Definition: ClumpParticle.h:250
std::vector< Mdouble > getPebbleRadius() const
Definition: ClumpParticle.h:157
Matrix3D transpose(Matrix3D M)
Definition: ClumpParticle.h:318
void addPebble(Vec3D position, Mdouble radius)
Definition: ClumpParticle.cc:129
std::vector< ClumpParticle * > pebbleParticles_
Definition: ClumpParticle.h:307
std::string getTypeVTK(unsigned i) const override
Definition: ClumpParticle.h:255
Vec3D getInitPrincipalDirections_e3() const
Definition: ClumpParticle.h:117
void setClump(ClumpParticle *master)
Definition: ClumpParticle.h:167
Vec3D getPrincipalDirections_e1() const
Definition: ClumpParticle.h:94
Matrix3D principalDirections_
Definition: ClumpParticle.h:303
Implementation of a 3D matrix.
Definition: Matrix.h:38
Mdouble YX
Definition: Matrix.h:43
Mdouble ZX
Definition: Matrix.h:43
Mdouble XY
Definition: Matrix.h:43
Mdouble YY
Definition: Matrix.h:43
Mdouble ZY
Definition: Matrix.h:43
Mdouble ZZ
Definition: Matrix.h:43
Mdouble YZ
Definition: Matrix.h:43
Mdouble XZ
Definition: Matrix.h:43
Mdouble XX
all nine matrix elements
Definition: Matrix.h:43
Implementation of a 3D symmetric matrix.
Definition: MatrixSymmetric.h:37
Base class for all non-spherical particle types.
Definition: NonSphericalParticle.h:37
Definition: ParticleSpecies.h:37
Definition: Vector.h:51
Mdouble Y
Definition: Vector.h:66
Mdouble Z
Definition: Vector.h:66
static Mdouble getLengthSquared(const Vec3D &a)
Calculates the squared length of a Vec3D: .
Definition: Vector.h:332
Mdouble X
the vector components
Definition: Vector.h:66
static Mdouble getLength(const Vec3D &a)
Calculates the length of a Vec3D: .
Definition: Vector.cc:331
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:76
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51