MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BaseInteractable.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, 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 
27 #include "Particles/BaseParticle.h"
28 #include "SpeciesHandler.h"
29 
30 #ifdef MERCURY_USE_OMP
31  #include "omp.h"
32 #endif
33 
43  BaseObject()
44 {
45  //setting all vectors to zero
50  force_.setZero();
51  torque_.setZero();
52  //setting the default species to species 0
53  indSpecies_ = 0;
54  //setting all relevant pointers to nullptr
55  species_ = nullptr;
56  prescribedPosition_ = nullptr;
57  prescribedVelocity_ = nullptr;
58  prescribedOrientation_ = nullptr;
60  logger(DEBUG, "BaseInteractable::BaseInteractable() finished");
61 }
62 
76  : BaseObject(p)
77 {
78  interactions_.clear();
79  position_ = p.position_;
81  velocity_ = p.velocity_;
83  force_ = p.force_;
84  torque_ = p.torque_;
85  species_ = p.species_;
91  logger(DEBUG, "BaseInteractable::BaseInteractable(const BaseInteractable &p finished");
92 }
93 
98 {
99  logger(VERBOSE, "Deleting BaseInteractable with index= % and id = %, size = %", getIndex(), getId(),
100  interactions_.size());
101  while (!interactions_.empty())
102  {
103  interactions_.front()->removeFromHandler();
104  }
105  logger(DEBUG, "BaseInteractable::~BaseInteractable() finished");
106 }
107 
108 
116 void BaseInteractable::addForce(const Vec3D& addForce)
117 {
118  if (OMP_THREAD_NUM==0) {
119  force_ += addForce;
120  } else {
122  }
123 }
124 
132 void BaseInteractable::addTorque(const Vec3D& addTorque)
133 {
134  if (OMP_THREAD_NUM==0) {
135  torque_ += addTorque;
136  } else {
138  }
139 }
140 
141 void BaseInteractable::resetForceTorque(int numberOfOMPthreads)
142 {
143  #ifdef MERCURY_USE_OMP
144  forceOMP_.resize(numberOfOMPthreads-1);
145  torqueOMP_.resize(numberOfOMPthreads-1);
146  for (auto& force : forceOMP_)
147  {
148  force.setZero();
149  }
150 
151  for (auto& torque : torqueOMP_)
152  {
153  torque.setZero();
154  }
155  #endif
156 
157  force_.setZero();
158  torque_.setZero();
159 }
160 
161 
163 {
164  #ifdef MERCURY_USE_OMP
165  for (const auto force : forceOMP_)
166  {
167  force_ += force;
168  }
169 
170  for (const auto torque : torqueOMP_)
171  {
172  torque_ += torque;
173  }
174  #endif
175 }
176 
186 {
187  logger.assert(species->getHandler() != nullptr && species->getHandler()->getObject(species->getIndex())==species, "Error: Species is not part of any handler yet");
188  species_ = species;
189  indSpecies_ = species->getIndex();
190 }
191 
200 {
202 }
203 
205 {
206  orientation_.setEuler(eulerAngle);
207 }
208 
215 void BaseInteractable::move(const Vec3D& move)
216 {
217  position_ += move;
218 }
219 
230 void BaseInteractable::rotate(const Vec3D& angularVelocityDt)
231 {
232  //if (!angularVelocityDt.isZero()) {
233  orientation_.updateAngularDisplacement(angularVelocityDt);
234  //}
235 }
236 
244 void BaseInteractable::read(std::istream& is)
245 {
246  BaseObject::read(is);
247  std::string dummy;
248  is >> dummy >> indSpecies_;
249  is >> dummy >> position_;
250  is >> dummy;
251  //this if-statement is added to read Kasper van der Vaart's data files, which contain an additional variable named positionAbsolute
252  if (dummy == "positionAbsolute")
253  {
254  is >> dummy >> dummy >> dummy >> dummy;
255  }
256  is >> orientation_;
257  is >> dummy >> velocity_;
258  is >> dummy >> angularVelocity_;
259  is >> dummy;
260  if (dummy == "0")
261  is >> dummy;
262  is >> force_;
263  is >> dummy >> torque_;
264 }
265 
274 void BaseInteractable::write(std::ostream& os) const
275 {
276  BaseObject::write(os);
277  os << " indSpecies " << indSpecies_
278  << " position " << position_
279  << " orientation " << orientation_
280  //<< " axis " << orientation_.getAxis()
281  << " velocity " << velocity_
282  << " angularVelocity " << angularVelocity_
283  << " force " << force_
284  << " torque " << torque_;
285 }
286 
293 {
294  interactions_.push_back(I);
295 }
296 
309 {
310  for (std::vector<BaseInteraction*>::iterator it = interactions_.begin(); it != interactions_.end(); ++it)
311  {
312  if (I == (*it))
313  {
314  interactions_.erase(it);
315  return true;
316  }
317  }
318  logger(WARN, "Error in BaseInteractable::removeInteraction: Interaction could not be removed");
319  return false;
320 }
321 
330 {
331  return velocity_;
332 }
333 
342 {
343  return angularVelocity_;
344 }
345 
351 {
352  velocity_ = velocity;
353 }
354 
360 void BaseInteractable::setAngularVelocity(const Vec3D& angularVelocity)
361 {
362  angularVelocity_ = angularVelocity;
363 }
364 
370 void BaseInteractable::addAngularVelocity(const Vec3D& angularVelocity)
371 {
372  angularVelocity_ += angularVelocity;
373 }
374 
376 {
377  return getVelocity() - Vec3D::cross(contact - getPosition(), getAngularVelocity());
378 }
379 
387 {
388  for (BaseInteraction* interaction : pOriginal.interactions_)
389  {
390  //So here this is the ghost and it is the interaction of the ghost/
391  interaction->copySwitchPointer(&pOriginal, this);
392  }
393 }
394 
413 void BaseInteractable::setPrescribedPosition(const std::function<Vec3D(double)>& prescribedPosition)
414 {
415  prescribedPosition_ = prescribedPosition;
416 }
417 
424 {
426  {
428  }
429 }
430 
444 void BaseInteractable::setPrescribedVelocity(const std::function<Vec3D(double)>& prescribedVelocity)
445 {
446  prescribedVelocity_ = prescribedVelocity;
447 }
448 
455 {
457  {
459  }
460 }
461 
474 void BaseInteractable::setPrescribedOrientation(const std::function<Quaternion(double)>& prescribedOrientation)
475 {
476  prescribedOrientation_ = prescribedOrientation;
477 }
478 
486 {
488  {
490  }
491 }
492 
493 /*
494  * \details This functions sets a std::function which specifies the
495  * AngularVeclocity of an interactable. This function takes a double
496  * the time and returns a Vec3D which is the new angular velocity.
497  * See also BaseInteractable::setPrescribedOrientation
498  * Note this functions works the same way as
499  * BaseInteractable::setPosition and BaseInteractable::setVeclocity.
500  */
501 void BaseInteractable::setPrescribedAngularVelocity(const std::function<Vec3D(double)>& prescribedAngularVelocity)
502 {
503  prescribedAngularVelocity_ = prescribedAngularVelocity;
504 }
505 
513 {
515  {
517  }
518 }
519 
538 void BaseInteractable::integrateBeforeForceComputation(double time, double timeStep)
539 {
541  {
543  {
544  //Both the velocity and position are defined; as we are using leap
545  //frog method so the velocity is evaluated half a time later.
546  applyPrescribedPosition(time + timeStep);
547  applyPrescribedVelocity(time + 0.5 * timeStep);
548  }
549  else
550  {
551  //Only the position is defined.
552  //Velocity is evaluated from a finite different of the Position
553  //Note, we use 0.5 +- 0.1 timeStep for the velocity eval.
554  applyPrescribedPosition(time + timeStep);
555  setVelocity((prescribedPosition_(time + 0.6 * timeStep) - prescribedPosition_(time + 0.4 * timeStep)) /
556  (0.2 * timeStep));
557  }
558  }
559  else
560  {
562  {
563  //Only the velocity is set. The position is calculated from the
564  //the integral of velocity.
565  applyPrescribedVelocity(time + 0.5 * timeStep);
566  move(getVelocity() * timeStep);
567  }
568  else
569  {
570  //Neither is set move based on the computed velocity of the object.
571  move(getVelocity() * timeStep);
572  }
573  }
575  {
577  {
578  applyPrescribedOrientation(time + timeStep);
579  applyPrescribedAngularVelocity(time + 0.5 * timeStep);
580  }
581  else
582  {
583  applyPrescribedOrientation(time + timeStep);
584  setAngularVelocity(getOrientation().applyCInverse(
585  (prescribedOrientation_(time + 0.6 * timeStep) - prescribedOrientation_(time + 0.4 * timeStep)) /
586  (0.2 * timeStep)
587  ));
588  }
589  }
590  else
591  {
593  {
594  applyPrescribedAngularVelocity(time + 0.5 * timeStep);
595  rotate(getAngularVelocity() * timeStep);
596  }
597  else
598  {
599  rotate(getAngularVelocity() * timeStep);
600  }
601  }
602 }
603 
611 void BaseInteractable::integrateAfterForceComputation(double time, double timeStep)
612 {
614  {
616  {
617  applyPrescribedVelocity(time + timeStep);
618  }
619  else
620  {
621  setVelocity((prescribedPosition_(time + 1.1 * timeStep) - prescribedPosition_(time + 0.9 * timeStep)) /
622  (0.2 * timeStep));
623  }
624  }
625  else
626  {
628  {
629  applyPrescribedVelocity(time + 0.5 * timeStep);
630  }
631  }
633  {
635  {
636  applyPrescribedAngularVelocity(time + timeStep);
637  }
638  else
639  {
640  setAngularVelocity(getOrientation().applyCInverse(
641  (prescribedOrientation_(time + 1.1 * timeStep) - prescribedOrientation_(time + 0.9 * timeStep)) /
642  (0.2 * timeStep)
643  ));
644  }
645  }
646  else
647  {
649  {
650  applyPrescribedAngularVelocity(time + 0.5 * timeStep);
651  }
652  }
653 }
654 
void setPrescribedVelocity(const std::function< Vec3D(double)> &prescribedVelocity)
Allows the velocity of an infinite mass interactable to be prescribed.
Implementation of a 3D quaternion (by Vitaliy).
Definition: Quaternion.h:62
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:125
unsigned int getIndex() const
Returns the index of the object in the handler.
Definition: BaseObject.h:118
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
void copySwitchPointer(const BaseInteractable *original, BaseInteractable *ghost) const
This copies the interactions of the original particle and replaces the original with the ghost copy...
std::function< Vec3D(double)> prescribedAngularVelocity_
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
unsigned int indSpecies_
void applyPrescribedPosition(double time)
Computes the position from the user defined prescribed position function.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
void addForce(const Vec3D &addForce)
Adds an amount to the force on this BaseInteractable.
void resetForceTorque(int numberOfOMPthreads)
bool removeInteraction(BaseInteraction *I)
Removes an interaction from this BaseInteractable.
void addInteraction(BaseInteraction *I)
Adds an interaction to this BaseInteractable.
void applyPrescribedAngularVelocity(double time)
Computes the angular velocity from the user defined prescribed angular velocity.
It is an abstract base class due to the purely virtual functions declared below. Even if the function...
Definition: BaseObject.h:50
void setPrescribedAngularVelocity(const std::function< Vec3D(double)> &prescribedAngularVelocity)
Allows the angular velocity of the infinite mass interactable to be prescribed.
virtual const Vec3D getVelocityAtContact(const Vec3D &contact) const
Returns the velocity at the contact point, use by many force laws.
void applyPrescribedOrientation(double time)
Computes the orientation from the user defined prescribed orientation function.
void integrateBeforeForceComputation(double time, double timeStep)
This is part of integrate routine for objects with infinite mass.
void setZero()
Sets all elements to zero.
Definition: Vector.cc:43
std::function< Quaternion(double)> prescribedOrientation_
void setOrientationViaEuler(Vec3D eulerAngle)
Sets the orientation of this BaseInteractable by defining the euler angles.
void setPrescribedOrientation(const std::function< Quaternion(double)> &prescribedOrientation)
Allows the orientation of the infinite mass interactbale to be prescribed.
virtual const Vec3D & getAngularVelocity() const
Returns the angular velocity of this interactable.
std::function< Vec3D(double)> prescribedPosition_
void setPrescribedPosition(const std::function< Vec3D(double)> &prescribedPosition)
Allows the position of an infinite mass interactable to be prescribed.
Stores information about interactions between two interactable objects; often particles but could be ...
void addAngularVelocity(const Vec3D &angularVelocity)
add an increment to the angular velocity.
SpeciesHandler * getHandler() const
Returns the pointer to the handler to which this species belongs.
Definition: BaseSpecies.cc:99
void setSpecies(const ParticleSpecies *species)
Sets the species of this BaseInteractable.
void write(std::ostream &os) const override
Write a BaseInteractable to an output stream.
const ParticleSpecies * species_
std::vector< Vec3D > torqueOMP_
#define OMP_THREAD_NUM
Definition: GeneralDefine.h:69
std::function< Vec3D(double)> prescribedVelocity_
void setEuler(const Vec3D &e)
Convert Euler angles to a quaternion. See Wikipedia for details.
Definition: Quaternion.cc:473
Quaternion orientation_
std::vector< BaseInteraction * > interactions_
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
void setOrientationViaNormal(Vec3D normal)
Used to the the normal of an InfiniteWall that has a normal into the x-direction by default...
Definition: Quaternion.cc:536
static Vec3D cross(const Vec3D &a, const Vec3D &b)
Calculates the cross product of two Vec3D: .
Definition: Vector.cc:163
virtual void read(std::istream &is)=0
Definition: BaseObject.cc:81
void copyInteractionsForPeriodicParticles(const BaseInteractable &p)
Copies interactions to this BaseInteractable whenever a periodic copy made.
std::vector< Vec3D > forceOMP_
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
virtual void rotate(const Vec3D &angularVelocityDt)
Rotates this BaseInteractable.
void applyPrescribedVelocity(double time)
Computes the velocity from the user defined prescribed velocity function.
void setOrientationViaNormal(Vec3D normal)
Sets the orientation of this BaseInteractable by defining the vector that results from the rotation o...
Defines the basic properties that a interactable object can have.
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
void addTorque(const Vec3D &addTorque)
Adds an amount to the torque on this BaseInteractable.
virtual const Vec3D & getVelocity() const
Returns the velocity of this interactable.
void updateAngularDisplacement(Vec3D angularVelocityDt)
Definition: Quaternion.cc:431
Definition: Vector.h:49
virtual void move(const Vec3D &move)
Moves this BaseInteractable by adding an amount to the position.
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
void setAngularVelocity(const Vec3D &angularVelocity)
set the angular velocity of the BaseInteractble.
void read(std::istream &is) override
Reads a BaseInteractable from an input stream.
virtual void write(std::ostream &os) const =0
A purely virtual function which has an implementation which writes the name and the object id_ to the...
Definition: BaseObject.cc:91
BaseInteractable()
Default BaseInteractable constructor.
void integrateAfterForceComputation(double time, double timeStep)
This is part of the integration routine for objects with infinite mass.
void setUnity()
Sets quaternion value to (1,0,0,0)
Definition: Quaternion.cc:53
~BaseInteractable() override
Destructor, it simply destructs the BaseInteractable and all the objects it contains.