revision: v0.14
SuperQuadricParticle Class Referencefinal

#include <SuperQuadricParticle.h>

+ Inheritance diagram for SuperQuadricParticle:

Public Member Functions

 SuperQuadricParticle ()
 Basic Particle constructor, creates a superquadric with axes (1,1,1) and exponents (2,2), so it creates a sphere with radius 1. More...
 
 SuperQuadricParticle (const SuperQuadricParticle &p)
 Copy constructor, which accepts as input a reference to a Superquadric. It creates a copy of this Particle and all it's information. Usually it is better to use the copy() function for polymorphism. More...
 
 SuperQuadricParticle (const BaseParticle &p)
 
 ~SuperQuadricParticle () override
 Destructor, needs to be implemented and checked to see if it is the largest or smallest particle currently in its particleHandler. More...
 
SuperQuadricParticlecopy () const override
 Copy method. It calls to copy constructor of this superquadric, useful for polymorphism. More...
 
void write (std::ostream &os) const override
 Write function: write this superquadric to the given output-stream, for example a restart-file. More...
 
void read (std::istream &is) override
 Read function: read in the information for this superquadric from the given input-stream, for example a restart file. More...
 
std::string getName () const override
 Returns the name of the class, here "SuperQuadricParticle". More...
 
void setAxesAndExponents (const Mdouble &a1, const Mdouble &a2, const Mdouble &a3, const Mdouble &eps1, const Mdouble &eps2)
 Set the geometrical properties of the superquadrics, namely the axes-lengths a1, a2 and a3, and the exponents epsilon1 and epsilon2. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al. More...
 
void setAxesAndExponents (const Vec3D &axes, const Mdouble &eps1, const Mdouble &eps2)
 Set the geometrical properties of the superquadrics, namely the axes-lengths axes, and the exponents epsilon1 and epsilon2. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al. More...
 
void setAxes (const Mdouble &a1, const Mdouble &a2, const Mdouble &a3)
 Set the axes-lengths to a1, a2 and a3 for this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al. More...
 
void setAxes (const Vec3D &axes) override
 Set the axes-lengths to axes for this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al. More...
 
void setExponents (const Mdouble &eps1, const Mdouble &eps2) override
 Set the exponents to eps1 and eps2 for this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al. More...
 
Vec3D getAxes () const override
 Get the axes-lengths of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al. More...
 
Mdouble getExponentEps1 () const override
 Get the first exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al. More...
 
Mdouble getExponentEps2 () const override
 Get the second exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al. More...
 
Mdouble getVolume () const override
 
void setInertia () override
 Compute and set the inertia-tensor for this superquadric. For internal use only. More...
 
void setRadius (const Mdouble radius) override
 Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species) More...
 
BaseInteractiongetInteractionWith (BaseParticle *P, unsigned timeStamp, InteractionHandler *interactionHandler) override
 Checks if this superquadric is in interaction with the given particle, and if so, returns vector of pointer to the associated BaseInteraction object (else returns empty vector). More...
 
Mdouble getCurvature (const LabFixedCoordinates &labFixedCoordinates) const override
 Get the mean curvature of this superquadric at the given (lab-fixed) position, see Podlozhyuk et al. (2017) eq (39) More...
 
bool isInContactWith (const BaseParticle *p) const override
 Get whether or not this superquadric is in contact with the given particle. More...
 
SmallVector< 3 > computeShapeGradientLabFixed (const LabFixedCoordinates &labFixedCoordinates) const
 Compute and get the gradient of the shape-function at the given (lab-fixed) position. More...
 
BaseInteractiongetInteractionWithSuperQuad (SuperQuadricParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler)
 Checks if this superquadric is in interaction with the given superquadric, and if so, returns vector of pointer to the associated BaseInteraction object (else returns empty vector). More...
 
SmallMatrix< 3, 3 > computeHessianLabFixed (const LabFixedCoordinates &labFixedCoordinates) const
 Compute and get the hessian ("second derivative") of the shape-function at the given (lab-fixed) position. More...
 
Mdouble computeShape (const LabFixedCoordinates &labFixedCoordinates) const
 Compute and get the shape-functiion at the given (lab-fixed) position. More...
 
SmallVector< 4 > computeResidualContactDetection (const SmallVector< 4 > &position, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
 Objective function for contact detection between the two given superquadrics. See Podlozhyuk et al. (2017) eq (22). More...
 
SmallMatrix< 4, 4 > getJacobianOfContactDetectionObjective (const SmallVector< 4 > &contactPoint, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
 Compute and return the derivative of functionThatShouldBecomeZeroForContactDetection, both to the position and the Lagrange multiplier, and evaluated at the contact point. More...
 
SmallVector< 4 > getInitialGuessForContact (const SuperQuadricParticle *pQuad, BaseInteraction *C) const
 Get an initial guess for the contact-point between this particle and the given particle. More...
 
Mdouble overlapFromContactPoint (const LabFixedCoordinates &contactPoint, const LabFixedCoordinates &normal) const
 Compute the distance between the contact-point and surface of this superquadric particle. More...
 
SmallVector< 4 > getContactPoint (const SuperQuadricParticle *p, BaseInteraction *C) const
 Compute the contact point between this and the given superquadric particle. More...
 
SmallVector< 4 > getContactPointPlanB (const SuperQuadricParticle *pOther, unsigned numberOfSteps) const
 If the "normal" procedure fails to find a contact point, use an alternative approach that involves starting with two spheres to compute the interaction, and becoming less and less spherical. More...
 
bool computeContactPoint (SmallVector< 4 > &contactPoint, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
 Perform the actual Newton-iterations to find the contact point. Note, that it is given back as a parameter. More...
 
void writeDebugMessageStep1 (const SuperQuadricParticle *pQuad, const SmallVector< 4 > &contactPointPlanB) const
 
void writeDebugMessageStep2 (const SuperQuadricParticle *pQuad, const Vec3D &dAxesThis, const Mdouble &dn11, const Mdouble &dn12, const Vec3D &dAxesOther, const Mdouble &dn21, const Mdouble &dn22) const
 
void writeDebugMessageStep3 (const Vec3D &axesThis, const Mdouble &n11, const Mdouble &n12, const Vec3D &axesOther, const Mdouble &n21, const Mdouble &n22) const
 
void writeDebugMessageMiddleOfLoop (const SuperQuadricParticle &p1, const SuperQuadricParticle &p2, SmallVector< 4 > &contactPointPlanB, const unsigned int &counter) const
 
Mdouble getInteractionRadius (const BaseParticle *particle) const
 returns the radius plus half the interactionDistance of the mixed species More...
 
void computeMass (const ParticleSpecies &s) override
 Computes the particle's (inverse) mass and inertia. More...
 
- Public Member Functions inherited from BaseParticle
 BaseParticle ()
 Basic Particle constructor, creates an Particle at (0,0,0) with radius, mass and inertia equal to 1. More...
 
 BaseParticle (const BaseParticle &p)
 Particle copy constructor, which accepts as input a reference to a Particle. It creates a copy of this Particle and all it's information. Usually it is better to use the copy() function for polymorphism. More...
 
 BaseParticle (const ParticleSpecies *s)
 
 ~BaseParticle () override
 Particle destructor, needs to be implemented and checked if it removes tangential spring information. More...
 
void fixParticle ()
 Fix Particle function. It fixes a Particle by setting its inverse mass and inertia and velocities to zero. More...
 
bool isFixed () const override
 Is fixed Particle function. It returns whether a Particle is fixed or not, by checking its inverse Mass. More...
 
bool isMPIParticle () const
 Indicates if this particle is a ghost in the MPI domain. More...
 
void setMPIParticle (bool flag)
 Flags the mpi particle status. More...
 
bool isInMPIDomain ()
 Indicates if the particle is in the communication zone of the mpi domain. More...
 
void setInMPIDomain (bool flag)
 Flags the status of the particle if wether it is in the communication zone or not. More...
 
bool isInPeriodicDomain () const
 Indicates if the particle is in the periodic boundary communication zone. More...
 
void setInPeriodicDomain (bool flag)
 Flags the status of the particle wether it is in the periodic communication zone or not. More...
 
bool isPeriodicGhostParticle () const
 Indicates if this particle is a ghost in the periodic boundary. More...
 
void setPeriodicGhostParticle (bool flag)
 Flags the status of the particle to be a ghost in periodic boundary or not. More...
 
bool isMaserParticle () const
 Indicates if this particle belongs to the maser boundary. More...
 
void setMaserParticle (bool flag)
 Flags the status of the particle if it belongs to the maser boundary or not. More...
 
void setCommunicationComplexity (unsigned complexity)
 Set the communication complexity of the particle. More...
 
unsigned getCommunicationComplexity ()
 Obtains the communication complexity of the particle. More...
 
void setPeriodicComplexity (std::vector< int > complexity)
 Set the periodic communication complexity of the particle. More...
 
void setPeriodicComplexity (int index, int value)
 Set the periodic communication complexity of the particle. More...
 
const std::vector< int > & getPeriodicComplexity ()
 Obtains the periodic communication complexity of the particle. More...
 
void setPreviousPeriodicComplexity (std::vector< int > complexity)
 Set the previous periodic communication complexity of the paritcle. More...
 
const std::vector< int > & getPreviousPeriodicComplexity () const
 Sets the previous periodic communication complexity of the particle. More...
 
int getPeriodicComplexity (int index)
 Gets the periodic communication complexity of a certain boundary. More...
 
void unfix ()
 Unfix Particle function, which required a reference to the Species vector. It unfixes a Particle by computing the Particles mass and inertia. More...
 
virtual void oldRead (std::istream &is)
 
virtual void setInfo (Mdouble info)
 Sets some user-defined information about this object (by default, species ID). More...
 
virtual Mdouble getInfo () const
 Returns some user-defined information about this object (by default, species ID). More...
 
void printHGrid (std::ostream &os) const
 Adds particle's HGrid level and cell coordinates to an ostream. More...
 
unsigned int getHGridLevel () const
 Returns particle's HGrid level. More...
 
BaseParticlegetHGridNextObject () const
 Returns pointer to next object in particle's HGrid level & cell. More...
 
BaseParticlegetHGridPrevObject () const
 Returns pointer to previous object in particle's HGrid level & cell. More...
 
int getHGridX () const
 Returns particle's HGrid cell X-coordinate. More...
 
int getHGridY () const
 Returns particle's HGrid cell Y-coordinate. More...
 
int getHGridZ () const
 Returns particle's HGrid cell Z-coordinate. More...
 
MatrixSymmetric3D getInvInertia () const
 Returns the inverse of the particle's inertia tensor. More...
 
Mdouble getInvMass () const override
 Returns the inverse of the particle's mass. More...
 
Mdouble getKineticEnergy () const
 Calculates the particle's translational kinetic energy. More...
 
Mdouble getRotationalEnergy () const
 Calculates the particle's rotational kinetic energy. More...
 
Mdouble getGravitationalEnergy () const
 Calculates the particle's gravitational energy. More...
 
Mdouble getMass () const
 Returns the particle's mass. More...
 
Mdouble getSurfaceArea () const
 
Vec3D getMomentum () const
 
MatrixSymmetric3D getInertia () const
 
Vec3D getAngularMomentum () const
 
BaseParticlegetPeriodicFromParticle () const
 Returns the 'original' particle this one's a periodic copy of. More...
 
Mdouble getRadius () const
 Returns the particle's radius. More...
 
Mdouble getMaxInteractionRadius () const
 Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles) More...
 
Mdouble getInteractionDistance (const BaseInteractable *i) const
 Returns the interactionDistance_ of the mixed species of this particle and the particle or wall i. More...
 
Mdouble getSumOfInteractionRadii (const BaseParticle *particle) const
 returns the sum of the radii plus the interactionDistance More...
 
Mdouble getWallInteractionRadius (const BaseWall *wall) const
 returns the radius plus the interactionDistance More...
 
const Vec3DgetDisplacement () const
 Returns the particle's displacement relative to the previous time step. More...
 
const Vec3DgetPreviousPosition () const
 Returns the particle's position in the previous time step. More...
 
const Vec3D getDisplacement2 (Mdouble xmin, Mdouble xmax, Mdouble ymin, Mdouble ymax, Mdouble zmin, Mdouble zmax, Mdouble t) const
 
void setInertia (MatrixSymmetric3D inertia)
 Sets the particle's inertia_ (and adjusts invInertia_ accordingly) More...
 
void setInverseInertia (MatrixSymmetric3D inverseInertia)
 Sets the particle's inertia_ (and adjusts invInertia_ accordingly) More...
 
void setInfiniteInertia ()
 Sets the particle's inertia_ to 'infinite' (1e20) and its invInertia_ to 0. More...
 
void setPeriodicFromParticle (BaseParticle *p)
 Assigns the pointer to the 'original' particle this one's a periodic copy of (used in periodic boundary condition implementations). More...
 
void setHGridX (const int x)
 Sets the particle's HGrid cell X-coordinate. More...
 
void setHGridY (const int y)
 Sets the particle's HGrid cell Y-coordinate. More...
 
void setHGridZ (const int z)
 Sets the particle's HGrid cell Z-coordinate. More...
 
void setHGridLevel (const unsigned int level)
 Sets the particle's HGrid level. More...
 
void setHGridNextObject (BaseParticle *p)
 Sets the pointer to the next object in the particle's HGrid cell & level. More...
 
void setHGridPrevObject (BaseParticle *p)
 Sets the pointer to the previous object in the particle's HGrid cell & level. More...
 
MERCURY_DEPRECATED void setMass (Mdouble mass)
 Sets the particle's mass. More...
 
void setMassForP3Statistics (Mdouble mass)
 Sets the particle's mass This function should not be used, but is necessary to extend the CG toolbox to non-spherical particles. More...
 
void setDisplacement (const Vec3D &disp)
 Sets the particle's displacement (= difference between current position and that of the previous time step) More...
 
void setPreviousPosition (const Vec3D &pos)
 Sets the particle's position in the previous time step. More...
 
void movePrevious (const Vec3D &posMove)
 Adds a vector to the particle's previousPosition_. More...
 
void accelerate (const Vec3D &vel)
 Increases the particle's velocity_ by the given vector. More...
 
void angularAccelerate (const Vec3D &angVel)
 Increases the particle's angularVelocity_ by the given vector. More...
 
void addDisplacement (const Vec3D &addDisp)
 Adds a vector to the particle's displacement_. More...
 
void setHandler (ParticleHandler *handler)
 Sets the pointer to the particle's ParticleHandler. More...
 
ParticleHandlergetHandler () const
 Returns pointer to the particle's ParticleHandler. More...
 
void integrateBeforeForceComputation (double time, double timeStep)
 First step of Velocity Verlet integration. More...
 
void integrateAfterForceComputation (double time, double timeStep)
 Second step of Velocity Verlet integration. More...
 
unsigned int getParticleDimensions () const
 Returns the particle's dimensions (either 2 or 3). More...
 
MERCURY_DEPRECATED void setIndSpecies (unsigned int indSpecies) override
 
void setSpecies (const ParticleSpecies *species)
 
virtual unsigned getNumberOfFieldsVTK () const
 
virtual std::string getTypeVTK (unsigned i) const
 
virtual std::string getNameVTK (unsigned i) const
 
virtual std::vector< MdoublegetFieldVTK (unsigned i) const
 
virtual void actionsAfterTimeStep ()
 
virtual bool isSphericalParticle () const
 
const HGridCellgetHGridCell () const
 
- Public Member Functions inherited from BaseInteractable
 BaseInteractable ()
 Default BaseInteractable constructor. More...
 
 BaseInteractable (const BaseInteractable &p)
 Copy constructor. More...
 
 ~BaseInteractable () override
 Destructor, it simply destructs the BaseInteractable and all the objects it contains. More...
 
unsigned int getIndSpecies () const
 Returns the index of the species associated with the interactable object. More...
 
const ParticleSpeciesgetSpecies () const
 Returns a pointer to the species of this BaseInteractable. More...
 
void setSpecies (const ParticleSpecies *species)
 Sets the species of this BaseInteractable. More...
 
const Vec3DgetForce () const
 Returns the force on this BaseInteractable. More...
 
const Vec3DgetTorque () const
 Returns the torque on this BaseInteractable. More...
 
void setForce (const Vec3D &force)
 Sets the force on this BaseInteractable. More...
 
void setTorque (const Vec3D &torque)
 Sets the torque on this BaseInteractable. More...
 
void addForce (const Vec3D &addForce)
 Adds an amount to the force on this BaseInteractable. More...
 
void addTorque (const Vec3D &addTorque)
 Adds an amount to the torque on this BaseInteractable. More...
 
void resetForceTorque (int numberOfOMPthreads)
 
void sumForceTorqueOMP ()
 
const Vec3DgetPosition () const
 Returns the position of this BaseInteractable. More...
 
const QuaterniongetOrientation () const
 Returns the orientation of this BaseInteractable. More...
 
void setPosition (const Vec3D &position)
 Sets the position of this BaseInteractable. More...
 
void setOrientationViaNormal (Vec3D normal)
 Sets the orientation of this BaseInteractable by defining the vector that results from the rotation of the (1,0,0) vector. More...
 
void setOrientationViaEuler (Vec3D eulerAngle)
 Sets the orientation of this BaseInteractable by defining the euler angles. More...
 
void setOrientation (const Quaternion &orientation)
 Sets the orientation of this BaseInteractable. More...
 
virtual void move (const Vec3D &move)
 Moves this BaseInteractable by adding an amount to the position. More...
 
virtual void rotate (const Vec3D &angularVelocityDt)
 Rotates this BaseInteractable. More...
 
const std::vector< BaseInteraction * > & getInteractions () const
 Returns a list of interactions which belong to this interactable. More...
 
void addInteraction (BaseInteraction *I)
 Adds an interaction to this BaseInteractable. More...
 
bool removeInteraction (BaseInteraction *I)
 Removes an interaction from this BaseInteractable. More...
 
void copyInteractionsForPeriodicParticles (const BaseInteractable &p)
 Copies interactions to this BaseInteractable whenever a periodic copy made. More...
 
void setVelocity (const Vec3D &velocity)
 set the velocity of the BaseInteractable. More...
 
void setAngularVelocity (const Vec3D &angularVelocity)
 set the angular velocity of the BaseInteractble. More...
 
void addVelocity (const Vec3D &velocity)
 adds an increment to the velocity. More...
 
void addAngularVelocity (const Vec3D &angularVelocity)
 add an increment to the angular velocity. More...
 
virtual const Vec3DgetVelocity () const
 Returns the velocity of this interactable. More...
 
virtual const Vec3DgetAngularVelocity () const
 Returns the angular velocity of this interactable. More...
 
void setPrescribedPosition (const std::function< Vec3D(double)> &prescribedPosition)
 Allows the position of an infinite mass interactable to be prescribed. More...
 
void applyPrescribedPosition (double time)
 Computes the position from the user defined prescribed position function. More...
 
void setPrescribedVelocity (const std::function< Vec3D(double)> &prescribedVelocity)
 Allows the velocity of an infinite mass interactable to be prescribed. More...
 
void applyPrescribedVelocity (double time)
 Computes the velocity from the user defined prescribed velocity function. More...
 
void setPrescribedOrientation (const std::function< Quaternion(double)> &prescribedOrientation)
 Allows the orientation of the infinite mass interactbale to be prescribed. More...
 
void applyPrescribedOrientation (double time)
 Computes the orientation from the user defined prescribed orientation function. More...
 
void setPrescribedAngularVelocity (const std::function< Vec3D(double)> &prescribedAngularVelocity)
 Allows the angular velocity of the infinite mass interactable to be prescribed. More...
 
void applyPrescribedAngularVelocity (double time)
 Computes the angular velocity from the user defined prescribed angular velocity. More...
 
virtual const Vec3D getVelocityAtContact (const Vec3D &contact) const
 Returns the velocity at the contact point, use by many force laws. More...
 
void integrateBeforeForceComputation (double time, double timeStep)
 This is part of integrate routine for objects with infinite mass. More...
 
void integrateAfterForceComputation (double time, double timeStep)
 This is part of the integration routine for objects with infinite mass. More...
 
- Public Member Functions inherited from BaseObject
 BaseObject ()=default
 Default constructor. More...
 
 BaseObject (const BaseObject &p)=default
 Copy constructor, copies all the objects BaseObject contains. More...
 
virtual ~BaseObject ()=default
 virtual destructor More...
 
virtual void moveInHandler (unsigned int index)
 Except that it is virtual, it does the same thing as setIndex() does. More...
 
void setIndex (unsigned int index)
 Allows one to assign an index to an object in the handler/container. More...
 
void setId (unsigned long id)
 Assigns a unique identifier to each object in the handler (container) which remains constant even after the object is deleted from the container/handler. More...
 
unsigned int getIndex () const
 Returns the index of the object in the handler. More...
 
unsigned int getId () const
 Returns the unique identifier of any particular object. More...
 
void setGroupId (unsigned groupId)
 
unsigned getGroupId () const
 

Private Member Functions

void setBoundingRadius ()
 Get the radius of the sphere that fits precisely around the particle. More...
 

Private Attributes

Mdouble eps1_
 Blockiness parameters. More...
 
Mdouble eps2_
 
Vec3D axes_
 Lengths of principal axes (a1, a2, a3). More...
 

Additional Inherited Members

- Protected Attributes inherited from BaseParticle
Mdouble radius_
 
Mdouble invMass_
 Particle radius_. More...
 
MatrixSymmetric3D invInertia_
 Inverse Particle mass (for computation optimization) More...
 

Constructor & Destructor Documentation

◆ SuperQuadricParticle() [1/3]

SuperQuadricParticle::SuperQuadricParticle ( )

Basic Particle constructor, creates a superquadric with axes (1,1,1) and exponents (2,2), so it creates a sphere with radius 1.

calls the default constructor of BaseParticle, and creates an SuperEllipsoid with axes (1,1,1) and exponents (2,2), so it creates a sphere with radius 1.

39  : BaseParticle()
40 {
41  axes_ = Vec3D(1.0, 1.0, 1.0);
42  eps1_ = 1.0;
43  eps2_ = 1.0;
44  logger(DEBUG, "SuperQuadricParticle::SuperQuadricParticle() finished");
45 }

References axes_, eps1_, eps2_, FATAL, and logger.

Referenced by copy(), getInteractionWith(), and isInContactWith().

◆ SuperQuadricParticle() [2/3]

SuperQuadricParticle::SuperQuadricParticle ( const SuperQuadricParticle p)

Copy constructor, which accepts as input a reference to a Superquadric. It creates a copy of this Particle and all it's information. Usually it is better to use the copy() function for polymorphism.

Constructor that copies most of the properties of the given particle. Please note that not everything is copied, for example the position in the HGrid is not determined yet by the end of this constructor. It also does not copy the interactions and the pointer to the handler that handles this particle. Use with care.

Parameters
[in,out]pReference to the SuperQuad this one should become a copy of.
56  : BaseParticle(p)
57 {
58  axes_ = p.axes_;
59  eps1_ = p.eps1_;
60  eps2_ = p.eps2_;
61 }

References axes_, eps1_, and eps2_.

◆ SuperQuadricParticle() [3/3]

SuperQuadricParticle::SuperQuadricParticle ( const BaseParticle p)
64  : BaseParticle(p)
65 {
66  Mdouble radius = p.getRadius();
67  axes_ = Vec3D(radius, radius, radius);
68  eps1_ = 1.0;
69  eps2_ = 1.0;
70 }

References axes_, eps1_, eps2_, and BaseParticle::getRadius().

◆ ~SuperQuadricParticle()

SuperQuadricParticle::~SuperQuadricParticle ( )
override

Destructor, needs to be implemented and checked to see if it is the largest or smallest particle currently in its particleHandler.

Destructor. It asks the ParticleHandler to check if this was the smallest or largest particle and adjust itself accordingly.

77 {
78  if (getHandler() != nullptr)
79  {
81  }
82  logger(DEBUG, "SuperQuadricParticle::SuperQuadricParticleParticle() of particle % finished.", getId());
83 
84 }

References ParticleHandler::checkExtremaOnDelete(), FATAL, BaseParticle::getHandler(), BaseObject::getId(), and logger.

Member Function Documentation

◆ computeContactPoint()

bool SuperQuadricParticle::computeContactPoint ( SmallVector< 4 > &  contactPoint,
const SuperQuadricParticle p1,
const SuperQuadricParticle p2 
) const

Perform the actual Newton-iterations to find the contact point. Note, that it is given back as a parameter.

Function that actually performs the Newton iterations for contact-detection. We use a Newton-method with adaptive damping, i.e. we start with an undamped Newton iteration, and if we "overshoot" we use a damping-factor that gets halved until there is no overshoot anymore. If the damping-factor gets too small, either plan B is initiated, or if this is already plan B, we give an error. After each iteration, the damping factor is set to 1 again.

Parameters
[in|out]contactPoint The contact point we are looking for. The input is an approximate contact point that is obtained by another method, and we update this point until we reach the defined tolerance.
[in]p1The first particle of the contact we are looking for
[in]p2The second particle of the contact we are looking for
Returns
Boolean for whether or not the contact-detection was successful.
819 {
820  // Damped Newton's method: (dampingFactor 1 is undamped)
821  Mdouble dampingFactor = 1;
822  int iteration = 1;
823  const Mdouble tolerance = 1e-5;
824  const unsigned maxIterations = 100;
825 
826  logger(VERBOSE, "func to be zero value: % \n",
827  computeResidualContactDetection(contactPoint, p1, p2));
828 
829  while (computeResidualContactDetection(contactPoint, p1, p2).length() > tolerance && iteration < maxIterations)
830  {
831  logger(VERBOSE, "Iteration: %\n", iteration);
832 
833  SmallMatrix<4, 4> jacobian = getJacobianOfContactDetectionObjective(contactPoint, p1, p2);
834  SmallVector<4> func = -computeResidualContactDetection(contactPoint, p1, p2);
835  jacobian.solve(func); //note: solve is in place
836 
837  SmallVector<4> newPoint = contactPoint + dampingFactor * func;
838  const auto residueOld = computeResidualContactDetection(contactPoint, p1, p2);
839  const auto residueNew = computeResidualContactDetection(newPoint, p1, p2);
840  logger(VERBOSE, "Old value: % (%) new value: % (%)",
841  computeResidualContactDetection(contactPoint, p1, p2),
842  computeResidualContactDetection(contactPoint, p1, p2).length(),
843  computeResidualContactDetection(newPoint, p1, p2),
844  computeResidualContactDetection(newPoint, p1, p2).length());
845 
846  logger(VERBOSE, "current contact point: %, new contact point: %\n", contactPoint, newPoint);
847  logger(VERBOSE, "damping factor: %", dampingFactor);
848 
849  if (residueNew.length()
850  < residueOld.length())
851  {
852  contactPoint = newPoint;
853  dampingFactor = 1;
854  }
855  else
856  {
857  dampingFactor /= 2;
858 
859  if (dampingFactor < 1e-10)
860  {
861  return false;
862  }
863  }
864 
865  iteration++;
866  }
867  return (iteration != maxIterations);
868 }

References computeResidualContactDetection(), getJacobianOfContactDetectionObjective(), logger, SmallMatrix< numberOfRows, numberOfColumns >::solve(), and VERBOSE.

Referenced by getContactPoint(), and getContactPointPlanB().

◆ computeHessianLabFixed()

SmallMatrix< 3, 3 > SuperQuadricParticle::computeHessianLabFixed ( const LabFixedCoordinates labFixedCoordinates) const

Compute and get the hessian ("second derivative") of the shape-function at the given (lab-fixed) position.

This function computes the Hessian ("second derivative") of the shape function in the body-fixed coordinate system. The expressions are provided in Eq. 15 of the article in Comp. Part. Mech. (2017) 4 : 101-118.

Returns
The hessian of the shape function at the given (lab-fixed) coordinates. Note, that this is the hessian to the lab-fixed coordinates, H_X (F)(X)
Todo:
Come up with good expression for when x = y = 0 and n1 < n2
528 {
529  SmallMatrix<3, 3> hessian;
530  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
531  getOrientation().rotateBack(bodyFixedCoordinates);
532 
533  const Mdouble n1 = 2.0 / eps1_;
534  const Mdouble n2 = 2.0 / eps2_;
535 
536  const Mdouble absXoa = std::abs(bodyFixedCoordinates.X / axes_.X);
537  const Mdouble absYob = std::abs(bodyFixedCoordinates.Y / axes_.Y);
538  const Mdouble absZoc = std::abs(bodyFixedCoordinates.Z / axes_.Z);
539 
540  const Mdouble nu = std::pow(absXoa, n2) + std::pow(absYob, n2);
541  const Mdouble help1 = std::pow(nu, n1 / n2 - 1.0);
542  const Mdouble help2 = std::pow(nu, n1 / n2 - 2.0);
543 
544  hessian(0, 0) = n1 * (n2 - 1) / axes_.X / axes_.X * std::pow(absXoa, n2 - 2.0) * help1
545  + ((n1 - n2) * n1 / axes_.X / axes_.X) * std::pow(absXoa, 2. * n2 - 2.0) * help2;
546  hessian(1, 1) = n1 * (n2 - 1) / axes_.Y / axes_.Y * std::pow(absYob, n2 - 2.0) * help1
547  + ((n1 - n2) * n1 / axes_.Y / axes_.Y) * std::pow(absYob, 2. * n2 - 2.0) * help2;
548  hessian(2, 2) = n1 * (n1 - 1) / axes_.Z / axes_.Z * std::pow(absZoc, n1 - 2.0);
549  hessian(1, 0) = n1 * (n1 - n2) / axes_.X / axes_.Y * std::pow(absXoa, n2 - 1.0) * std::pow(absYob, n2 - 1.0) *
550  help2 * mathsFunc::sign(bodyFixedCoordinates.X * bodyFixedCoordinates.Y);
551  hessian(0, 1) = hessian(1, 0);
554  return A * hessian * (A.transpose());
555 }

References A, axes_, eps1_, eps2_, BaseInteractable::getOrientation(), BaseInteractable::getPosition(), Quaternion::getRotationMatrix(), Quaternion::rotateBack(), mathsFunc::sign(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getCurvature(), getJacobianOfContactDetectionObjective(), ShapeGradientHessianTester::testCushion(), ShapeGradientHessianTester::testEllipsoid(), ShapeGradientHessianTester::testRoundedBeam(), and ShapeGradientHessianTester::testSphere().

◆ computeMass()

void SuperQuadricParticle::computeMass ( const ParticleSpecies s)
overridevirtual

Computes the particle's (inverse) mass and inertia.

Reimplemented from BaseParticle.

984  {
985  if (isFixed()) return;
986  if (getParticleDimensions()==3) {
987  Mdouble volume = getVolume();
988 
989  Mdouble help1 = mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_);
990  Mdouble help2 = mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0);
991  Mdouble help3 = mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_ + 1);
992  Mdouble help4 = mathsFunc::beta(1.5 * eps1_, eps1_ + 1.0);
993 
994  invMass_ = 1.0 / (volume * s.getDensity());
995  invInertia_.XX = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
996  (axes_.Y * axes_.Y * help1 * help2
997  + 4.0 * axes_.Z * axes_.Z * help3 * help4));
998  invInertia_.XY = 0.0;
999  invInertia_.XZ = 0.0;
1000  invInertia_.YY = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
1001  (axes_.X * axes_.X * help1 * help2
1002  + 4.0 * axes_.Z * axes_.Z * help3 * help4));
1003  invInertia_.YZ = 0.0;
1004  invInertia_.ZZ = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
1005  (axes_.X * axes_.X + axes_.Y * axes_.Y) * help1 * help2);
1006  } else {
1007  logger(ERROR, "[SuperQuadricParticle::computeMass] SuperQuadricParticle cannot be two-dimensional (yet)");
1008  }
1009 };

References axes_, mathsFunc::beta(), eps1_, eps2_, ERROR, ParticleSpecies::getDensity(), BaseParticle::getParticleDimensions(), getVolume(), BaseParticle::invInertia_, BaseParticle::invMass_, BaseParticle::isFixed(), logger, Vec3D::X, MatrixSymmetric3D::XX, MatrixSymmetric3D::XY, MatrixSymmetric3D::XZ, Vec3D::Y, MatrixSymmetric3D::YY, MatrixSymmetric3D::YZ, Vec3D::Z, and MatrixSymmetric3D::ZZ.

◆ computeResidualContactDetection()

SmallVector< 4 > SuperQuadricParticle::computeResidualContactDetection ( const SmallVector< 4 > &  position,
const SuperQuadricParticle p1,
const SuperQuadricParticle p2 
) const

Objective function for contact detection between the two given superquadrics. See Podlozhyuk et al. (2017) eq (22).

For the contact detection, we formulate the optimisation problem "minimise F1 + F2, s.t. F1 = F2", where F1 and F2 are the shape functions of particles 1 and 2. Define a Lagrange multiplier mu^2, then this is equivalent to solving "Gradient1 + mu^2 Gradient2 = 0, F1 - F2 = 0". The left-hand side of that function is computed in this function. See also Comp. Part. Mech. (2017) 4 : 101-118, equation 22.

Parameters
[in]positionCurrent guess for the contact-point, where the expression above should be evaluated
[in]p1First particle for which we are looking for a contact with.
[in]p2Second particle for which we are looking for a contact with.
Returns
Residual of the objective function.
570 {
571  LabFixedCoordinates labFixedCoordinates;
572  labFixedCoordinates.X = position[0];
573  labFixedCoordinates.Y = position[1];
574  labFixedCoordinates.Z = position[2];
575  const Mdouble lagrangeMultiplier = position[3];
576 
577  // First compute the gradient of the superquadric shape function and then rotate it to the
578  // global frame of reference
579  SmallVector<3> gradThis = p1->computeShapeGradientLabFixed(labFixedCoordinates);
580 
581  SmallVector<3> gradOther = p2->computeShapeGradientLabFixed(labFixedCoordinates);
582 
583  // Eqn. (23) in Podhlozhnyuk et al, Comp. Part. Mech. 2017.
584  // where contactPoint[3] = mu
585  SmallVector<4> toBecomeZero;
586  for (unsigned int i = 0; i < 3; ++i)
587  {
588  toBecomeZero[i] = gradThis[i] + lagrangeMultiplier * lagrangeMultiplier * gradOther[i];
589  }
590  toBecomeZero[3] = p1->computeShape(labFixedCoordinates) - p2->computeShape(labFixedCoordinates);
591 
592  return toBecomeZero;
593 }

References computeShape(), computeShapeGradientLabFixed(), constants::i, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by computeContactPoint().

◆ computeShape()

Mdouble SuperQuadricParticle::computeShape ( const LabFixedCoordinates labFixedCoordinates) const

Compute and get the shape-functiion at the given (lab-fixed) position.

This function computes the value of the shape function in the lab-fixed coordinate system. The expression is provided in Section 2.3 of the article in Comp. Part. Mech. (2017) 4 : 101-118.

Returns
The value of the shape-function at the given (lab-fixed) coordinates.
475 {
476  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
477  getOrientation().rotateBack(bodyFixedCoordinates);
478 
479  const Mdouble n1 = 2. / eps1_;
480  const Mdouble n2 = 2. / eps2_;
481 
482  return std::pow(std::pow(std::abs(bodyFixedCoordinates.X / axes_.X), n2)
483  + std::pow(std::abs(bodyFixedCoordinates.Y / axes_.Y), n2), n1 / n2)
484  + std::pow(std::abs(bodyFixedCoordinates.Z / axes_.Z), n1) - 1.0;
485 }

References axes_, eps1_, eps2_, BaseInteractable::getOrientation(), BaseInteractable::getPosition(), Quaternion::rotateBack(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by computeResidualContactDetection(), getInteractionWithSuperQuad(), isInContactWith(), overlapFromContactPoint(), ShapeGradientHessianTester::testCushion(), ShapeGradientHessianTester::testEllipsoid(), ShapeGradientHessianTester::testRoundedBeam(), and ShapeGradientHessianTester::testSphere().

◆ computeShapeGradientLabFixed()

SmallVector< 3 > SuperQuadricParticle::computeShapeGradientLabFixed ( const LabFixedCoordinates labFixedCoordinates) const

Compute and get the gradient of the shape-function at the given (lab-fixed) position.

This function computes the gradient ("first derivative") of the shape function in the lab-fixed coordinate system. The expressions are provided in Eq. 14 of the article in Comp. Part. Mech. (2017) 4 : 101-118.

Returns
The gradient of the shape function at the given (lab-fixed) coordinates. Note, that this is the gradient to the lab-fixed coordinates, \nabla_X F(X)
Todo:
Come up with good expression for when x = y = 0 and n1 < n2
495 {
496 
497  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
498  getOrientation().rotateBack(bodyFixedCoordinates);
499 
500  const Mdouble n1 = 2. / eps1_;
501  const Mdouble n2 = 2. / eps2_;
502 
503  const Mdouble absXoa = std::abs(bodyFixedCoordinates.X / axes_.X);
504  const Mdouble absYob = std::abs(bodyFixedCoordinates.Y / axes_.Y);
505  const Mdouble absZoc = std::abs(bodyFixedCoordinates.Z / axes_.Z);
506 
507  const Mdouble nu = std::pow(absXoa, n2) + std::pow(absYob, n2);
508 
509  const Mdouble help1 = std::pow(nu, n1 / n2 - 1.0);
510 
511  Vec3D gradientBodyFixed = {
512  (n1 / axes_.X) * std::pow(absXoa, n2 - 1.0) * help1 * mathsFunc::sign(bodyFixedCoordinates.X),
513  (n1 / axes_.Y) * std::pow(absYob, n2 - 1.0) * help1 * mathsFunc::sign(bodyFixedCoordinates.Y),
514  (n1 / axes_.Z) * std::pow(absZoc, n1 - 1.0) * mathsFunc::sign(bodyFixedCoordinates.Z)};
515  //we get the gradient in terms of the lab-fixed coordinates by \nabla_X F(X) = A \nabla_x f(x).
516  getOrientation().rotate(gradientBodyFixed);
517  return {gradientBodyFixed.X, gradientBodyFixed.Y, gradientBodyFixed.Z};
518 }

References axes_, eps1_, eps2_, BaseInteractable::getOrientation(), BaseInteractable::getPosition(), Quaternion::rotate(), Quaternion::rotateBack(), mathsFunc::sign(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by computeResidualContactDetection(), getCurvature(), getInteractionWithSuperQuad(), getJacobianOfContactDetectionObjective(), overlapFromContactPoint(), ShapeGradientHessianTester::testCushion(), ShapeGradientHessianTester::testEllipsoid(), ShapeGradientHessianTester::testRoundedBeam(), and ShapeGradientHessianTester::testSphere().

◆ copy()

SuperQuadricParticle * SuperQuadricParticle::copy ( ) const
overridevirtual

Copy method. It calls to copy constructor of this superquadric, useful for polymorphism.

Copy method. Uses copy constructor to create a copy on the heap. Useful for polymorphism.

Returns
pointer to the particle's copy

Implements BaseParticle.

91 {
92  return new SuperQuadricParticle(*this);
93 }

References SuperQuadricParticle().

Referenced by EllipsoidsBouncingOnWallDemo::setupInitialConditions(), EllipticalSuperQuadricCollision::setupInitialConditions(), SlidingSpheresUnitTest::setupInitialConditions(), SphericalSuperQuadricCollision::setupInitialConditions(), ContactDetectionNormalSpheresTest::setupInitialConditions(), ContactDetectionRotatedSpheresTest::setupInitialConditions(), and ContactDetectionTester::setupParticles().

◆ getAxes()

Vec3D SuperQuadricParticle::getAxes ( ) const
overridevirtual

Get the axes-lengths of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al.

Todo:

TW we could remove this function from the BaseParticle and use a dynamic_cast instead

ID Middle-term plan is to template the BaseParticle on shape-type, so that we won't have to cast etc.

Reimplemented from BaseParticle.

183 {
184  return axes_;
185 }

References axes_.

Referenced by getContactPointPlanB(), InfiniteWall::getDistanceNormalOverlapSuperquadric(), writeDebugMessageStep1(), and writeDebugMessageStep2().

◆ getContactPoint()

SmallVector< 4 > SuperQuadricParticle::getContactPoint ( const SuperQuadricParticle p,
BaseInteraction C 
) const

Compute the contact point between this and the given superquadric particle.

Computes the contact point between this and the given superquadric particle. It first gets a first guess for the contact point, and then uses that to compute the real contact-point with Newton-iterations.

Parameters
[in]pthe particle with which we want to know the contact point with.
[in]CThe contact between this and the give superquadric. Does not contain information of this time-step yet. This contact is used to see if there was an interaction between p and this superquadric in the time-step before, and if so, that contact point can serve as an initial point for contact-detection between the particles.
Returns
A SmallVector of length 4, with the contact point and the lagrange-multiplier needed for Newton's method for finding the contact point.
373 {
374  //Assuming contact between two spheres
375  SmallVector<4> approximateContactPoint = getInitialGuessForContact(p, C);
376  if (computeContactPoint(approximateContactPoint, this, p))
377  {
378  return approximateContactPoint;
379  }
380  else
381  {
382  return getContactPointPlanB(p, 4);
383  }
384 }

References computeContactPoint(), getContactPointPlanB(), and getInitialGuessForContact().

Referenced by getInteractionWithSuperQuad(), and isInContactWith().

◆ getContactPointPlanB()

SmallVector< 4 > SuperQuadricParticle::getContactPointPlanB ( const SuperQuadricParticle pOther,
unsigned  numberOfSteps 
) const

If the "normal" procedure fails to find a contact point, use an alternative approach that involves starting with two spheres to compute the interaction, and becoming less and less spherical.

If the "normal" method of finding a contact point diverges, then we start plan B. For this, we approximate the current and given particle by spheres, and move more and more towards the actual shape for the particles in small increments. This way, we have a relatively good starting point for each new optimisation problem, and therefore the chance that our contact-detection algorithm diverges is much smaller. This is based on equations (24) - (27) of Comp. Part. Mech. (2017) 4 : 101-118.

Parameters
[in]pOtherParticle we want to know if the contact-point with
Returns
The point where the shape-functions of both particles are minimised.
722 {
723  logger(VERBOSE, "Number of iterations: %", numberOfSteps);
724  // Step 1: Compute contact point for two volume equivalent spheres
725  const Mdouble interactionRadiusThis = getInteractionRadius(pOther);
726  const Mdouble interactionRadiusOther = pOther->getInteractionRadius(this);
727  const Vec3D positionThis = getPosition();
728  const Vec3D positionOther = pOther->getPosition();
729 
730  //analytical solution for contact point between bounding spheres
731  const Vec3D initialPoint = (interactionRadiusOther * positionThis + interactionRadiusThis * positionOther) /
732  (interactionRadiusThis + interactionRadiusOther);
733 
734  SmallVector<4> contactPointPlanB;
735  // Analytical solution
736  contactPointPlanB[0] = initialPoint.X;
737  contactPointPlanB[1] = initialPoint.Y;
738  contactPointPlanB[2] = initialPoint.Z;
739  contactPointPlanB[3] = 1.0; //Need to check.
740 
741  writeDebugMessageStep1(pOther, contactPointPlanB);
742 
743  //Step 2: Compute the deltas
744  const Vec3D dAxesThis = (getAxes() - interactionRadiusThis * Vec3D(1, 1, 1)) / numberOfSteps;
745  const Mdouble dn11 = (2.0 / getExponentEps1() - 2.0) / numberOfSteps;
746  const Mdouble dn12 = (2.0 / getExponentEps2() - 2.0) / numberOfSteps;
747 
748  const Vec3D dAxesOther = (pOther->getAxes() - interactionRadiusOther * Vec3D(1, 1, 1)) / numberOfSteps;
749  const Mdouble dn21 = (2.0 / pOther->getExponentEps1() - 2.0) / numberOfSteps;
750  const Mdouble dn22 = (2.0 / pOther->getExponentEps2() - 2.0) / numberOfSteps;
751 
752  writeDebugMessageStep2(pOther, dAxesThis, dn11, dn12, dAxesOther, dn21, dn22);
753 
754  // Create two superquadrics with the above parameters for incrementally computing the contact point
757 
758  p1.setPosition(getPosition());
759  p2.setPosition(pOther->getPosition());
760 
762  p2.setOrientation(pOther->getOrientation());
763  bool success = true;
764  const unsigned maxIterations = 20;
765  for (unsigned int counter = 1; counter <= numberOfSteps; ++counter)
766  {
767 
768  // Step 3: Increment the shape and blockiness parameters
769  const Vec3D axesThis = interactionRadiusThis * Vec3D(1, 1, 1) + counter * dAxesThis;
770  const Mdouble n11 = 2.0 + counter * dn11;
771  const Mdouble n12 = 2.0 + counter * dn12;
772 
773  Vec3D axesOther = interactionRadiusOther * Vec3D(1, 1, 1) + counter * dAxesOther;
774  const Mdouble n21 = 2.0 + counter * dn21;
775  const Mdouble n22 = 2.0 + counter * dn22;
776 
777  p1.setAxesAndExponents(axesThis, 2.0 / n11, 2.0 / n12);
778  p2.setAxesAndExponents(axesOther, 2.0 / n21, 2.0 / n22);
779 
780  writeDebugMessageStep3(axesThis, n11, n12, axesOther, n21, n22);
781 
782  writeDebugMessageMiddleOfLoop(p1, p2, contactPointPlanB, counter);
783 
784  //compute the contact point of the new particles
785  success = computeContactPoint(contactPointPlanB, &p1, &p2);
786  if (!success)
787  {
788  if (numberOfSteps > maxIterations)
789  {
790  write(std::cout);
791  pOther->write(std::cout);
792  logger(ERROR, "Plan B fails even with more than 20 intermediate steps");
793  }
794  return getContactPointPlanB(pOther, 2 * numberOfSteps);
795  }
796  }
797  logger.assert_debug(p1.getAxes().X == getAxes().X, "In getContactPointPlanB, final particle for contact-detection not "
798  "the same as original particle");
799 
800  logger(VERBOSE, "Plan B contact point: %", contactPointPlanB);
801 
802  return contactPointPlanB;
803 }

References computeContactPoint(), ERROR, getAxes(), getExponentEps1(), getExponentEps2(), getInteractionRadius(), BaseInteractable::getOrientation(), BaseInteractable::getPosition(), logger, setAxesAndExponents(), BaseInteractable::setOrientation(), BaseInteractable::setPosition(), VERBOSE, write(), writeDebugMessageMiddleOfLoop(), writeDebugMessageStep1(), writeDebugMessageStep2(), writeDebugMessageStep3(), Vec3D::X, X, Vec3D::Y, and Vec3D::Z.

Referenced by getContactPoint().

◆ getCurvature()

Mdouble SuperQuadricParticle::getCurvature ( const LabFixedCoordinates labFixedCoordinates) const
overridevirtual

Get the mean curvature of this superquadric at the given (lab-fixed) position, see Podlozhyuk et al. (2017) eq (39)

Compute the mean curvature, see Comp. Part. Mech. (2017) 4 : 101-118, eq (39)

Parameters
[in]labFixedCoordinatesposition in lab fixed coordinate system
Returns
mean curvature of particle at labFixedCoordinates
Todo:
Minus sign the wrong way around, check why

Reimplemented from BaseParticle.

666 {
667  SmallVector<3> gradientVec = computeShapeGradientLabFixed(labFixedCoordinates);
668  SmallMatrix<3, 1> gradient = gradientVec;
669  SmallMatrix<3, 3> hessian = computeHessianLabFixed(labFixedCoordinates);
670  Mdouble helper = ((gradient.transpose() * hessian) * gradient)(0, 0);
671  Mdouble gradientLength = gradientVec.length();
672  Mdouble helper2 = gradientLength * gradientLength * (hessian(0, 0) + hessian(1, 1) + hessian(2, 2));
673  Mdouble denominator = 2 * gradientLength * gradientLength * gradientLength;
674 
676  return (helper2 - helper) / denominator;
677 }

References computeHessianLabFixed(), computeShapeGradientLabFixed(), SmallVector< numberOfRows >::length(), and SmallMatrix< numberOfRows, numberOfColumns >::transpose().

◆ getExponentEps1()

Mdouble SuperQuadricParticle::getExponentEps1 ( ) const
overridevirtual

Get the first exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al.

Reimplemented from BaseParticle.

188 {
189  return eps1_;
190 }

References eps1_.

Referenced by getContactPointPlanB(), InfiniteWall::getDistanceNormalOverlapSuperquadric(), and writeDebugMessageStep2().

◆ getExponentEps2()

Mdouble SuperQuadricParticle::getExponentEps2 ( ) const
overridevirtual

Get the second exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al.

Reimplemented from BaseParticle.

193 {
194  return eps2_;
195 }

References eps2_.

Referenced by getContactPointPlanB(), InfiniteWall::getDistanceNormalOverlapSuperquadric(), and writeDebugMessageStep2().

◆ getInitialGuessForContact()

SmallVector< 4 > SuperQuadricParticle::getInitialGuessForContact ( const SuperQuadricParticle pQuad,
BaseInteraction C 
) const

Get an initial guess for the contact-point between this particle and the given particle.

For the contact-detection, we need an initial guess where the contact will be (Newton's method only converges if you start sufficiently close). If there was already a contact during the last time-step, the values of the contact-point of last time-step is taken, otherwise it is taken as the middle between both particle-positions.

Parameters
pQuadother superquadric for which we want to know if there is a contact
Ccontact between this particle and the given particle
Returns
A guess for the initial contact: contact point X,Y,Z, lagrange multiplier. Note that these are given in the LAB-FIXED coordinates.
437 {
438  SmallVector<4> approximateContactPoint;
439  if (C == nullptr || C->getOverlap() < 1e-15)
440  {
441  //this is a new interaction
442  const LabFixedCoordinates branchVector = pQuad->getPosition() - getPosition();
443  //Get the square of the distance between particle i and particle j
444  const Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
445  const Mdouble distance = sqrt(distanceSquared);
446  const LabFixedCoordinates normal = (branchVector / distance);
447  const Mdouble overlap = getSumOfInteractionRadii(pQuad) - distance;
448  const LabFixedCoordinates contactPoint = (pQuad->getPosition() -
449  (pQuad->getInteractionRadius(this) - 0.5 * overlap) * normal);
450 
451  approximateContactPoint[0] = contactPoint.X;
452  approximateContactPoint[1] = contactPoint.Y;
453  approximateContactPoint[2] = contactPoint.Z;
454  approximateContactPoint[3] = 1;
455  }
456  else
457  {
458  //this is an existing interaction
459  const LabFixedCoordinates contactPoint = C->getContactPoint();
460  const Mdouble multiplier = C->getLagrangeMultiplier();
461  approximateContactPoint[0] = contactPoint.X;
462  approximateContactPoint[1] = contactPoint.Y;
463  approximateContactPoint[2] = contactPoint.Z;
464  approximateContactPoint[3] = multiplier;
465  }
466  return approximateContactPoint;
467 }

References BaseInteraction::getContactPoint(), getInteractionRadius(), BaseInteraction::getLagrangeMultiplier(), Vec3D::getLengthSquared(), BaseInteraction::getOverlap(), BaseInteractable::getPosition(), BaseParticle::getSumOfInteractionRadii(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getContactPoint().

◆ getInteractionRadius()

Mdouble SuperQuadricParticle::getInteractionRadius ( const BaseParticle particle) const

returns the radius plus half the interactionDistance of the mixed species

978 {
979  const auto mixedSpecies = getSpecies()->getHandler()->getMixedObject(getSpecies(),particle->getSpecies());
980  return getRadius() + 0.5 * mixedSpecies->getInteractionDistance();
981 }

References BaseSpecies::getHandler(), SpeciesHandler::getMixedObject(), BaseParticle::getRadius(), and BaseInteractable::getSpecies().

Referenced by getContactPointPlanB(), getInitialGuessForContact(), and writeDebugMessageStep1().

◆ getInteractionWith()

BaseInteraction * SuperQuadricParticle::getInteractionWith ( BaseParticle P,
unsigned  timeStamp,
InteractionHandler interactionHandler 
)
overridevirtual

Checks if this superquadric is in interaction with the given particle, and if so, returns vector of pointer to the associated BaseInteraction object (else returns empty vector).

Overwrites BaseInteractable::getInteractionWith. First checks if the bounding radii overlap, and if so, changes the given particle into a superquadric and calls getInteractionWithSuperQuadric to obtain the relevant contact data.

Parameters
pPointer to the particle we want to get the interaction with, can be any type of particle
timeStampTime stamp to be assigned to the interaction object (i.e., the current time)
interactionHandlerBaseInteraction container from where the interaction is retrieved, and to which it is assigned (if it is a new interaction).
Returns
A vector with size 1 (if there is an interaction) or 0 (if there is no interaction).

Reimplemented from BaseParticle.

288 {
289  //get the normal (from P away from the contact)
290  const LabFixedCoordinates branchVector = p->getPosition() - getPosition();
291  //Get the square of the distance between particle i and particle j
292  const Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
293  auto mixedSpecies = getSpecies()->getHandler()->getMixedObject(getSpecies(),p->getSpecies());
294 
295  const Mdouble sumOfInteractionRadii = getMaxInteractionRadius()+p->getMaxInteractionRadius();
296  if (distanceSquared < (sumOfInteractionRadii * sumOfInteractionRadii))
297  {
298  if (p->isSphericalParticle())
299  {
301  pQuad->setAxes(p->getRadius(), p->getRadius(), p->getRadius());
302  BaseInteraction* contacts = getInteractionWithSuperQuad(pQuad, timeStamp, interactionHandler);
303  delete pQuad;
304  return contacts;
305  } else {
306  SuperQuadricParticle* pQuad = static_cast<SuperQuadricParticle*>(p);
307  return getInteractionWithSuperQuad(pQuad, timeStamp, interactionHandler);
308  }
309  }
310  return nullptr;
311 }

References BaseSpecies::getHandler(), getInteractionWithSuperQuad(), Vec3D::getLengthSquared(), BaseParticle::getMaxInteractionRadius(), SpeciesHandler::getMixedObject(), BaseInteractable::getPosition(), BaseParticle::getRadius(), BaseInteractable::getSpecies(), BaseParticle::isSphericalParticle(), setAxes(), and SuperQuadricParticle().

Referenced by ContactDetectionTester::testEllipsoidsContact(), and ContactDetectionTester::testSpheresContact().

◆ getInteractionWithSuperQuad()

BaseInteraction * SuperQuadricParticle::getInteractionWithSuperQuad ( SuperQuadricParticle p,
unsigned  timeStamp,
InteractionHandler interactionHandler 
)

Checks if this superquadric is in interaction with the given superquadric, and if so, returns vector of pointer to the associated BaseInteraction object (else returns empty vector).

Computes all relevant information for the contact between this and a superquadric particle. Returns nullptr if there is no (positive) overlap.

Parameters
[in]pSuperquadric particle we want to have the contact-information for
[in]timeStampTime stamp to be assigned to the interaction object (i.e., the current time)
[in]interactionHandlerBaseInteraction container from where the interaction is retrieved, and to which it is assigned (if it is a new interaction).
Returns
A vector of Interaction* with size 1 (if there is an interaction) or 0 (if there is no interaction).
Todo:
: find correct value for 'distance'
325 {
326  BaseInteraction* const C = interactionHandler->getInteraction(p, this, timeStamp);
327  const SmallVector<4> approximateContactPoint = getContactPoint(p, C);
328 
329  //set the new contact point:
330  LabFixedCoordinates contactPoint;
331  contactPoint.X = approximateContactPoint[0];
332  contactPoint.Y = approximateContactPoint[1];
333  contactPoint.Z = approximateContactPoint[2];
334 
335  //if the minimum is positive, there is no contact: return empty vector
336  if (computeShape(contactPoint) > -1e-10)
337  {
338  interactionHandler->removeObject(C->getIndex());
339  return nullptr;
340  }
341 
342  C->setContactPoint(contactPoint);
343  C->setLagrangeMultiplier(approximateContactPoint[3]);
344 
345  SmallVector<3> gradThis = computeShapeGradientLabFixed(contactPoint).getNormalised();
346  LabFixedCoordinates normal;
347  normal.X = gradThis[0];
348  normal.Y = gradThis[1];
349  normal.Z = gradThis[2];
350  C->setNormal(normal);
351 
352  const Mdouble alphaI = overlapFromContactPoint(contactPoint, normal);
353  const Mdouble alphaJ = p->overlapFromContactPoint(contactPoint, -normal);
354  C->setOverlap(alphaI + alphaJ);
356  C->setDistance((getPosition() - p->getPosition()).getLength() - C->getOverlap());
357 
358  return C;
359 }

References computeShape(), computeShapeGradientLabFixed(), getContactPoint(), BaseObject::getIndex(), InteractionHandler::getInteraction(), BaseInteraction::getOverlap(), BaseInteractable::getPosition(), overlapFromContactPoint(), BaseHandler< T >::removeObject(), BaseInteraction::setContactPoint(), BaseInteraction::setDistance(), BaseInteraction::setLagrangeMultiplier(), BaseInteraction::setNormal(), BaseInteraction::setOverlap(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getInteractionWith().

◆ getJacobianOfContactDetectionObjective()

SmallMatrix< 4, 4 > SuperQuadricParticle::getJacobianOfContactDetectionObjective ( const SmallVector< 4 > &  contactPoint,
const SuperQuadricParticle p1,
const SuperQuadricParticle p2 
) const

Compute and return the derivative of functionThatShouldBecomeZeroForContactDetection, both to the position and the Lagrange multiplier, and evaluated at the contact point.

Compute and return the derivative of computeResidualContactDetection, both to the position and the Lagrange multiplier, and evaluated at the contact point. This is done in order to use Newton's method on computeResidualContactDetection, which comes from Comp. Part. Mech. (2017) 4 : 101-118.

The result is a 4x4 matrix, with 4 parts: [0,2]x[0,2] (upper-left corner): Hessian1 + mu^2 Hessian2, where Hessian is the "second derivative" of the shape function and mu^2 the Lagrange multiplier. This is the derivative of (22a) to x in Comp. Part. Mech. (2017) 4 : 101-118. [0,2]x[3] (upper-right corner): 2 * mu * Gradient2, where Gradient is the gradient of the shape function and mu is the Lagrange multiplier. This is the derivative of (22a) to mu in Comp. Part. Mech. (2017) 4 : 101-118. [3]x[0,2] (lower-left corner): Gradient1 - Gradient2, where Gradient is the gradient of the shape function. This is the derivative of (22b) to x in Comp. Part. Mech. (2017) 4 : 101-118. [3]x[3] (lower-right corner): 0. This is the derivative of (22b) to mu in Comp. Part. Mech. (2017) 4 : 101-118.

In order to compute these parts, first compute the first and second derivative of the shape function (the gradient and hessian) in the lab-fixed coordinate system for both particles. To do so, we apply the rotation matrix A on the body-fixed first and second derivative of the shape function, which are computed in the functions computeShapeGradientBodyFixed() and computeHessian(). The expressions for the purpose of mapping are provided in Eq. 18 of the article in Comp. Part. Mech. (2017) 4 : 101-118. Then fill in the respective contributions.

Parameters
[in]
627 {
628  LabFixedCoordinates labFixedCoordinates;
629  labFixedCoordinates.X = contactPoint[0];
630  labFixedCoordinates.Y = contactPoint[1];
631  labFixedCoordinates.Z = contactPoint[2];
632  const Mdouble lagrangeMultiplier = contactPoint[3];
633 
634  SmallVector<3> gradThis = p1->computeShapeGradientLabFixed(labFixedCoordinates);
635  SmallVector<3> gradOther = p2->computeShapeGradientLabFixed(labFixedCoordinates);
636 
637  SmallMatrix<3, 3> hessianThis = p1->computeHessianLabFixed(labFixedCoordinates);
638  SmallMatrix<3, 3> hessianOther = p2->computeHessianLabFixed(labFixedCoordinates);
639 
640  SmallMatrix<3, 3> upperLeftCorner = hessianThis + lagrangeMultiplier * lagrangeMultiplier * hessianOther;
641 
642  SmallMatrix<4, 4> jacobian;
643  for (unsigned int i = 0; i < 3; ++i)
644  {
645  for (unsigned int j = 0; j < 3; ++j)
646  {
647  jacobian(i, j) = upperLeftCorner(i, j);
648  }
649  jacobian(i, 3) = 2 * lagrangeMultiplier * gradOther[i];
650  }
651 
652  for (unsigned int j = 0; j < 3; ++j)
653  {
654  jacobian(3, j) = gradThis[j] - gradOther[j];
655  }
656 
657  return jacobian;
658 }

References computeHessianLabFixed(), computeShapeGradientLabFixed(), constants::i, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by computeContactPoint().

◆ getName()

std::string SuperQuadricParticle::getName ( ) const
overridevirtual

Returns the name of the class, here "SuperQuadricParticle".

Returns the name of the object; in this case "SuperQuadricParticle".

Returns
The std::string "SuperQuadricParticle".

Reimplemented from BaseParticle.

113 {
114  return "SuperQuadricParticle";
115 }

◆ getVolume()

Mdouble SuperQuadricParticle::getVolume ( ) const
overridevirtual

Get the volume of this superquadric.

Returns the volume of the SuperEllipsoid, which is calculated using its principal axis and the exponents. The analytical expressions are taken from Chapter 2 of the book "Segmentation and Recovery of Superquadrics" by Jaklic et al . However, the beta functions that are part of these expressions are approximations, see ExtendedMath.cc

Returns
The actual volume of this superquadric.

Reimplemented from BaseParticle.

205 {
206  logger.assert_debug(getHandler() != nullptr, "[SuperQuadricParticle::getVolume] no particle handler specified");
207 
208  return (2.0 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) * mathsFunc::beta(0.5 * eps1_ + 1.0, eps1_) *
209  mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_);
210 }

References axes_, mathsFunc::beta(), eps1_, eps2_, BaseParticle::getHandler(), logger, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by computeMass(), and VolumeTest::test().

◆ isInContactWith()

bool SuperQuadricParticle::isInContactWith ( const BaseParticle p) const
overridevirtual

Get whether or not this superquadric is in contact with the given particle.

Get whether or not this superquadric is in contact with the given particle: first transform the particle to a superquadric if necessary, then compute the contact-point using the function getContactPoint. If the shape-function at the contact point is negative, then there is a contact with the given particle, and otherwise not

Parameters
[in]pThe particle for which we want to know if there is a contact
Returns
True if there is a contact, false otherwise.

Reimplemented from BaseParticle.

687 {
688  SmallVector<4> approximateContactPoint;
689 
690  if (p->isSphericalParticle())
691  {
693  pQuad->setAxes(p->getRadius(), p->getRadius(), p->getRadius());
694  approximateContactPoint = getContactPoint(pQuad, nullptr);
695  delete pQuad;
696  } else {
697  const SuperQuadricParticle* pQuad = static_cast<const SuperQuadricParticle*>(p);
698  approximateContactPoint = getContactPoint(pQuad, nullptr);
699  }
700 
701  //set the new contact point:
702  LabFixedCoordinates contactPoint;
703  contactPoint.X = approximateContactPoint[0];
704  contactPoint.Y = approximateContactPoint[1];
705  contactPoint.Z = approximateContactPoint[2];
706 
707  //if the minimum is positive, there is no contact: return false
708  return (computeShape(contactPoint) < 0);
709 
710 }

References computeShape(), getContactPoint(), BaseParticle::getRadius(), BaseParticle::isSphericalParticle(), setAxes(), SuperQuadricParticle(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by ContactDetectionTester::testEllipsoidsContact().

◆ overlapFromContactPoint()

Mdouble SuperQuadricParticle::overlapFromContactPoint ( const LabFixedCoordinates contactPoint,
const LabFixedCoordinates normal 
) const

Compute the distance between the contact-point and surface of this superquadric particle.

Compute the distance between the contact-point and surface of this superquadric particle with Newton-iterations. This is the procedure as described by Comp. Part. Mech. (2017) 4 : 101-118, eq 37.

Parameters
[in]contactPointThe contact point between this particle with another BaseInteractable
[in]normalThe normal direction of the contact.
Returns
Distance between contact-point and surface of this particle
395 {
396  Mdouble alphaI = 0;
397  LabFixedCoordinates xEdge = contactPoint + alphaI * normal;
398  Mdouble dampingCoefficient = 1;
399  while (std::abs(computeShape(xEdge)) > 1e-10)
400  {
401  SmallVector<3> gradientShape = computeShapeGradientLabFixed(xEdge);
402  LabFixedCoordinates gradient;
403  gradient.X = gradientShape(0);
404  gradient.Y = gradientShape(1);
405  gradient.Z = gradientShape(2);
406  Mdouble newValue = alphaI - dampingCoefficient * computeShape(xEdge) / Vec3D::dot(gradient, normal);
407  const LabFixedCoordinates xEdgeNew = contactPoint + newValue * normal;
408 
409  if (std::abs(computeShape(xEdgeNew)) < std::abs(computeShape(xEdge)))
410  {
411  alphaI = newValue;
412  xEdge = xEdgeNew;
413  dampingCoefficient = 1;
414  }
415  else
416  {
417  dampingCoefficient /= 2;
418  if (dampingCoefficient < 1e-10)
419  {
420  logger(ERROR, "overlap detection does not converge");
421  }
422  }
423  }
424  return alphaI;
425 }

References computeShape(), computeShapeGradientLabFixed(), Vec3D::dot(), ERROR, logger, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getInteractionWithSuperQuad().

◆ read()

void SuperQuadricParticle::read ( std::istream &  is)
overridevirtual

Read function: read in the information for this superquadric from the given input-stream, for example a restart file.

Particle read function, which reads the axes_ and both epsilons from the given input-stream.

Parameters
[in,out]isinput stream with particle properties, e.g. a restart-file.

Reimplemented from BaseParticle.

123 {
124  BaseParticle::read(is);
125  std::string dummy;
126  is >> dummy >> axes_ >> dummy >> eps1_ >> dummy >> eps2_;
127  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
128  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
129 }

References axes_, eps1_, eps2_, logger, and BaseParticle::read().

◆ setAxes() [1/2]

void SuperQuadricParticle::setAxes ( const Mdouble a1,
const Mdouble a2,
const Mdouble a3 
)

Set the axes-lengths to a1, a2 and a3 for this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al.

165 {
166  setAxes({a1,a2,a3});
167 }

Referenced by getInteractionWith(), isInContactWith(), main(), setAxesAndExponents(), ShapesDemo::setupInitialConditions(), ContactDetectionTester::testEllipsoidsContact(), and ContactDetectionWithWallTester::testEllipsoidsContact().

◆ setAxes() [2/2]

void SuperQuadricParticle::setAxes ( const Vec3D axes)
overridevirtual

Set the axes-lengths to axes for this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al.

Reimplemented from BaseParticle.

155 {
156  axes_ = axes;
158  if (getSpecies() != nullptr)
159  {
160  getSpecies()->computeMass(this);
161  }
162 }

References axes_, ParticleSpecies::computeMass(), BaseInteractable::getSpecies(), and setBoundingRadius().

◆ setAxesAndExponents() [1/2]

void SuperQuadricParticle::setAxesAndExponents ( const Mdouble a1,
const Mdouble a2,
const Mdouble a3,
const Mdouble eps1,
const Mdouble eps2 
)

◆ setAxesAndExponents() [2/2]

void SuperQuadricParticle::setAxesAndExponents ( const Vec3D axes,
const Mdouble eps1,
const Mdouble eps2 
)

Set the geometrical properties of the superquadrics, namely the axes-lengths axes, and the exponents epsilon1 and epsilon2. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al.

145 {
146  eps1_ = eps1;
147  eps2_ = eps2;
148  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
149  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
150 
151  setAxes(axes);
152 }

References eps1_, eps2_, logger, and setAxes().

◆ setBoundingRadius()

void SuperQuadricParticle::setBoundingRadius ( )
private

Get the radius of the sphere that fits precisely around the particle.

Todo:
Currently only implemented for ellipsoids
245 {
246 
247  if(mathsFunc::isEqual(eps2_,1,std::numeric_limits<Mdouble>::epsilon()) && mathsFunc::isEqual(eps1_,1,std::numeric_limits<Mdouble>::epsilon()))
248  {
249  BaseParticle::setRadius(std::max(std::max(axes_.Y,axes_.X),axes_.Z));
250  return;
251  }else
252  {
253 
254  const Mdouble axesX = std::max(axes_.X, axes_.Y);
255  const Mdouble axesY = std::min(axes_.X, axes_.Y);
256 
257  Mdouble alpha;
258 
259  const Mdouble eps1 = std::min(.96, eps1_);
260  const Mdouble eps2 = std::min(.96, eps2_);
261 
262  alpha = std::pow(axesY / axesX, 2.0 / (2.0 / eps2 - 2.0));
263 
264  const Mdouble help1 = std::pow(alpha, 2.0 / eps2);
265  const Mdouble gamma = std::pow(1.0 + help1, eps2 / eps1 - 1.0);
266  const Mdouble beta = std::pow(gamma * axes_.Z * axes_.Z / (axesX * axesX), 1.0 / (2.0 / eps1 - 2.0));
267  const Mdouble xTilde = std::pow(std::pow(1 + help1, eps2 / eps1) + std::pow(beta, 2.0 / eps1),
268  -eps1 / 2.0);
269  BaseParticle::setRadius(std::sqrt(mathsFunc::square(axesX * xTilde)
270  + mathsFunc::square(alpha * axesY * xTilde)
271  + mathsFunc::square(beta * axes_.Z * xTilde)));
272 
273  }
274 }

References axes_, mathsFunc::beta(), eps1_, eps2_, mathsFunc::gamma(), mathsFunc::isEqual(), BaseParticle::setRadius(), mathsFunc::square(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by setAxes(), and setExponents().

◆ setExponents()

void SuperQuadricParticle::setExponents ( const Mdouble eps1,
const Mdouble eps2 
)
overridevirtual

Set the exponents to eps1 and eps2 for this superquadric. We use the super-ellipsoid definition stated in Chapter 2 of the book "Segmentation and recovery of superquadrics" by Jaklic et al.

Reimplemented from BaseParticle.

170 {
171  eps1_ = eps1;
172  eps2_ = eps2;
174  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
175  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
176  if (getSpecies() != nullptr)
177  {
178  getSpecies()->computeMass(this);
179  }
180 }

References ParticleSpecies::computeMass(), eps1_, eps2_, BaseInteractable::getSpecies(), logger, and setBoundingRadius().

Referenced by main(), and ShapesDemo::setupInitialConditions().

◆ setInertia()

void SuperQuadricParticle::setInertia ( )
overridevirtual

Compute and set the inertia-tensor for this superquadric. For internal use only.

This function computes the principal moments of inertia for the superellipsoids. Again, the analytical expressions are taken from Chapter 2 (pg. 36) of the book "Segmentation and Recovery of Superquadrics" by by Jaklic et al .

Reimplemented from BaseParticle.

217 {
218  MatrixSymmetric3D inertia;
219 
220  inertia.XX = getSpecies()->getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
221  (axes_.Y * axes_.Y * mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_) *
222  mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0)
223  + 4.0 * axes_.Z * axes_.Z * mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_ + 1) *
224  mathsFunc::beta(1.5 * eps1_, eps1_ + 1.0));
225 
226  inertia.YY = getSpecies()->getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
227  (axes_.X * axes_.X * mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_) *
228  mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0)
229  + 4.0 * axes_.Z * axes_.Z * mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_ + 1) *
230  mathsFunc::beta(1.5 * eps1_, eps1_ + 1.0));
231 
232  inertia.ZZ = getSpecies()->getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
233  (axes_.X * axes_.X + axes_.Y * axes_.Y) * mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_) *
234  mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0);
235 
236  BaseParticle::setInertia(inertia);
237 }

References axes_, mathsFunc::beta(), eps1_, eps2_, ParticleSpecies::getDensity(), BaseInteractable::getSpecies(), BaseParticle::setInertia(), Vec3D::X, MatrixSymmetric3D::XX, Vec3D::Y, MatrixSymmetric3D::YY, Vec3D::Z, and MatrixSymmetric3D::ZZ.

Referenced by BouncingSuperQuadric::setupInitialConditions(), EllipsoidsBouncingOnWallDemo::setupInitialConditions(), EllipticalSuperQuadricCollision::setupInitialConditions(), GranularCollapse::setupInitialConditions(), SlidingSpheresUnitTest::setupInitialConditions(), SphericalSuperQuadricCollision::setupInitialConditions(), ContactDetectionNormalSpheresTest::setupInitialConditions(), ContactDetectionRotatedSpheresTest::setupInitialConditions(), VisualisationTest::setupInitialConditions(), ContactDetectionWithWallTester::setupParticleAndWall(), and ContactDetectionTester::setupParticles().

◆ setRadius()

void SuperQuadricParticle::setRadius ( const Mdouble  radius)
overridevirtual

Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species)

Sets the radius of the particle, and from that computes the new mass (using its species) and checks whether it now is either the smallest or biggest particle in its ParticleHandler.

Parameters
[in]radiusthe new radius

Reimplemented from BaseParticle.

240 {
241  logger(ERROR,"This function should not be used");
242 }

References ERROR, and logger.

◆ write()

void SuperQuadricParticle::write ( std::ostream &  os) const
overridevirtual

Write function: write this superquadric to the given output-stream, for example a restart-file.

SuperQuadricParticle print method, which accepts an std::ostream as input. It prints human readable SuperQuadricParticle information to the given output-stream.

Parameters
[in,out]osstream to which the info is written, e.g. a restart-file or std::cout.

Reimplemented from BaseParticle.

101 {
103  os << " axes " << axes_
104  << " exp1 " << eps1_
105  << " exp2 " << eps2_;
106 }

References axes_, eps1_, eps2_, and BaseParticle::write().

Referenced by getContactPointPlanB().

◆ writeDebugMessageMiddleOfLoop()

void SuperQuadricParticle::writeDebugMessageMiddleOfLoop ( const SuperQuadricParticle p1,
const SuperQuadricParticle p2,
SmallVector< 4 > &  contactPointPlanB,
const unsigned int counter 
) const
872 {
873  logger(VERBOSE, "Position of particle 1 (p1): % \nPosition of particle 2 (p2): % \n",
874  p1.getPosition(), p2.getPosition());
875  logger(VERBOSE, "Orientation of particle 1 (p1): % \nOrientation of particle 2 (p2): %\n",
876  p1.getOrientation(), p2.getOrientation());
877 
878  // Step 4: Calculate the contact point for the updated shape parameters
879 
880  logger(VERBOSE, "----------------------------------------");
881  logger(VERBOSE, "STEP 4: Compute contact point");
882  logger(VERBOSE, "----------------------------------------");
883  logger(VERBOSE, " ");
884  logger(VERBOSE, "Counter: %", counter);
885  logger(VERBOSE, " ");
886 
887  logger(VERBOSE, "Analytical solution Contact Point: % ", contactPointPlanB);
888  logger(VERBOSE, " ");
889 }

References BaseInteractable::getOrientation(), BaseInteractable::getPosition(), logger, and VERBOSE.

Referenced by getContactPointPlanB().

◆ writeDebugMessageStep1()

void SuperQuadricParticle::writeDebugMessageStep1 ( const SuperQuadricParticle pQuad,
const SmallVector< 4 > &  contactPointPlanB 
) const
955 {
956  logger(VERBOSE, "---------------------");
957  logger(VERBOSE, "STEP 1");
958  logger(VERBOSE, "---------------------\n");
959 
960  logger(VERBOSE, "Position of particle 1: % ", getPosition());
961  logger(VERBOSE, "Position of particle 2 (pQuad): % ", pQuad->getPosition());
962  logger(VERBOSE, " ");
963 
964  logger(VERBOSE, "Radius of particle 1: % ", getInteractionRadius(pQuad));
965  logger(VERBOSE, "Radius of particle 2 (pQuad): % \n", pQuad->getInteractionRadius(this));
966 
967  logger(VERBOSE, "Orientation of particle 1: % ", getOrientation());
968  logger(VERBOSE, "Orientation of particle 2 (pQuad): % ", pQuad->getOrientation());
969  logger(VERBOSE, " ");
970 
971  logger(VERBOSE, "Particle 1 axes: %", getAxes());
972  logger(VERBOSE, "Particle 2 axes (pQuad): % \n", pQuad->getAxes());
973 
974  logger(VERBOSE, "Analytical solution for two equivalent spheres in contact: % \n", contactPointPlanB);
975 }

References getAxes(), getInteractionRadius(), BaseInteractable::getOrientation(), BaseInteractable::getPosition(), logger, and VERBOSE.

Referenced by getContactPointPlanB().

◆ writeDebugMessageStep2()

void SuperQuadricParticle::writeDebugMessageStep2 ( const SuperQuadricParticle pQuad,
const Vec3D dAxesThis,
const Mdouble dn11,
const Mdouble dn12,
const Vec3D dAxesOther,
const Mdouble dn21,
const Mdouble dn22 
) const
919 {
920  logger(VERBOSE, "---------------------");
921  logger(VERBOSE, "STEP 2");
922  logger(VERBOSE, "---------------------");
923  logger(VERBOSE, " ");
924 
925  logger(VERBOSE, "Particle 1 x-scale: %", pQuad->getAxes().X);
926  logger(VERBOSE, "Particle 1 y-scale: %", pQuad->getAxes().Y);
927  logger(VERBOSE, "Particle 1 z-scale: %", pQuad->getAxes().Z);
928  logger(VERBOSE, "Particle 1 n1: %", 2. / pQuad->getExponentEps1());
929  logger(VERBOSE, "Particle 1 n2: %", 2. / pQuad->getExponentEps2());
930  logger(VERBOSE, " ");
931 
932  logger(VERBOSE, "Particle 1 x-scale increment: %", dAxesThis.X);
933  logger(VERBOSE, "Particle 1 y-scale increment: %", dAxesThis.Y);
934  logger(VERBOSE, "Particle 1 z-scale increment: %", dAxesThis.Z);
935  logger(VERBOSE, "Particle 1 n1 increment: %", dn11);
936  logger(VERBOSE, "Particle 1 n2 increment: %", dn12);
937  logger(VERBOSE, " ");
938 
939  logger(VERBOSE, "Particle 2 x-scale: %", getAxes().X);
940  logger(VERBOSE, "Particle 2 y-scale: %", getAxes().Y);
941  logger(VERBOSE, "Particle 2 z-scale: %", getAxes().Z);
942  logger(VERBOSE, "Particle 2 n1: %", 2. / getExponentEps1());
943  logger(VERBOSE, "Particle 2 n2: %", 2. / getExponentEps2());
944  logger(VERBOSE, " ");
945 
946  logger(VERBOSE, "Particle 2 x-scale increment: %", dAxesOther.X);
947  logger(VERBOSE, "Particle 2 y-scale increment: %", dAxesOther.Y);
948  logger(VERBOSE, "Particle 2 z-scale increment: %", dAxesOther.Z);
949  logger(VERBOSE, "Particle 2 n1 increment: %", dn21);
950  logger(VERBOSE, "Particle 2 n2 increment: %", dn22);
951  logger(VERBOSE, " ");
952 }

References getAxes(), getExponentEps1(), getExponentEps2(), logger, VERBOSE, Vec3D::X, X, Vec3D::Y, Y, Vec3D::Z, and Z.

Referenced by getContactPointPlanB().

◆ writeDebugMessageStep3()

void SuperQuadricParticle::writeDebugMessageStep3 ( const Vec3D axesThis,
const Mdouble n11,
const Mdouble n12,
const Vec3D axesOther,
const Mdouble n21,
const Mdouble n22 
) const
894 {
895  logger(VERBOSE, "-----------------------------------");
896  logger(VERBOSE, "STEP 3: First increment");
897  logger(VERBOSE, "-----------------------------------");
898  logger(VERBOSE, " ");
899 
900  logger(VERBOSE, "Particle 1 x-scale after increment: %", axesThis.X);
901  logger(VERBOSE, "Particle 1 y-scale after increment: %", axesThis.Y);
902  logger(VERBOSE, "Particle 1 z-scale after increment: %", axesThis.Z);
903  logger(VERBOSE, "Particle 1 n1 after increment: %", n11);
904  logger(VERBOSE, "Particle 1 n2 after increment: %", n12);
905  logger(VERBOSE, " ");
906 
907  logger(VERBOSE, "Particle 2 x-scale after increment: %", axesOther.X);
908  logger(VERBOSE, "Particle 2 y-scale after increment: %", axesOther.Y);
909  logger(VERBOSE, "Particle 2 z-scale after increment: %", axesOther.Z);
910  logger(VERBOSE, "Particle 2 n1 after increment: %", n21);
911  logger(VERBOSE, "Particle 2 n2 after increment: %", n22);
912  logger(VERBOSE, " ");
913 }

References logger, VERBOSE, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getContactPointPlanB().

Member Data Documentation

◆ axes_

◆ eps1_

Mdouble SuperQuadricParticle::eps1_
private

Blockiness parameters.

Blockiness parameters should be in the range (0,1], where a sphere or ellipsoid is represented by eps1_ = eps2_ = 1.

Referenced by computeHessianLabFixed(), computeMass(), computeShape(), computeShapeGradientLabFixed(), getExponentEps1(), getVolume(), read(), setAxesAndExponents(), setBoundingRadius(), setExponents(), setInertia(), SuperQuadricParticle(), and write().

◆ eps2_


The documentation for this class was generated from the following files:
mathsFunc::square
T square(const T val)
squares a number
Definition: ExtendedMath.h:106
SuperQuadricParticle::getExponentEps2
Mdouble getExponentEps2() const override
Get the second exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter...
Definition: SuperQuadricParticle.cc:192
BaseInteractable::getSpecies
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
Definition: BaseInteractable.h:108
BaseParticle::invInertia_
MatrixSymmetric3D invInertia_
Inverse Particle mass (for computation optimization)
Definition: BaseParticle.h:657
Y
@ Y
Definition: StatisticsVector.h:42
SuperQuadricParticle::getJacobianOfContactDetectionObjective
SmallMatrix< 4, 4 > getJacobianOfContactDetectionObjective(const SmallVector< 4 > &contactPoint, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
Compute and return the derivative of functionThatShouldBecomeZeroForContactDetection,...
Definition: SuperQuadricParticle.cc:624
BaseSpecies::getHandler
SpeciesHandler * getHandler() const
Returns the pointer to the handler to which this species belongs.
Definition: BaseSpecies.cc:99
SuperQuadricParticle::computeResidualContactDetection
SmallVector< 4 > computeResidualContactDetection(const SmallVector< 4 > &position, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
Objective function for contact detection between the two given superquadrics. See Podlozhyuk et al....
Definition: SuperQuadricParticle.cc:567
BaseInteractable::setPosition
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Definition: BaseInteractable.h:239
SuperQuadricParticle::writeDebugMessageStep3
void writeDebugMessageStep3(const Vec3D &axesThis, const Mdouble &n11, const Mdouble &n12, const Vec3D &axesOther, const Mdouble &n21, const Mdouble &n22) const
Definition: SuperQuadricParticle.cc:892
logger
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
SmallMatrix::solve
void solve(SmallMatrix< numberOfRows, numberOfRightHandSideColumns > &B) const
solves Ax=B where A is the current matrix and B is passed in. The result is returned in B.
Definition: SmallMatrix_impl.h:329
SmallMatrix::transpose
SmallMatrix< numberOfColumns, numberOfRows > transpose() const
Definition: SmallMatrix.h:330
Vec3D::X
Mdouble X
the vector components
Definition: Vector.h:65
Vec3D::dot
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:76
SpeciesHandler::getMixedObject
std::enable_if<!std::is_pointer< typename U::MixedSpeciesType >::value, typename U::MixedSpeciesType * >::type getMixedObject(const U *S, const U *T)
Definition: SpeciesHandler.h:74
BaseInteractable::getOrientation
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Definition: BaseInteractable.h:230
BaseParticle::setInertia
virtual void setInertia()
Definition: BaseParticle.cc:497
SuperQuadricParticle::writeDebugMessageMiddleOfLoop
void writeDebugMessageMiddleOfLoop(const SuperQuadricParticle &p1, const SuperQuadricParticle &p2, SmallVector< 4 > &contactPointPlanB, const unsigned int &counter) const
Definition: SuperQuadricParticle.cc:870
BaseParticle::getParticleDimensions
unsigned int getParticleDimensions() const
Returns the particle's dimensions (either 2 or 3).
Definition: BaseParticle.cc:770
SuperQuadricParticle::axes_
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Definition: SuperQuadricParticle.h:297
BaseParticle::setRadius
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species)
Definition: BaseParticle.cc:542
Vec3D
Definition: Vector.h:50
ParticleSpecies::getDensity
Mdouble getDensity() const
Allows density_ to be accessed.
Definition: ParticleSpecies.cc:117
BaseInteraction::setDistance
void setDistance(Mdouble distance)
Sets the interaction distance between the two interacting objects.
Definition: BaseInteraction.cc:219
BaseInteractable::setOrientation
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
Definition: BaseInteractable.h:260
BaseInteraction::getContactPoint
const Vec3D & getContactPoint() const
Gets constant reference to contact point (vector).
Definition: BaseInteraction.h:234
A
@ A
Definition: StatisticsVector.h:42
SuperQuadricParticle::writeDebugMessageStep2
void writeDebugMessageStep2(const SuperQuadricParticle *pQuad, const Vec3D &dAxesThis, const Mdouble &dn11, const Mdouble &dn12, const Vec3D &dAxesOther, const Mdouble &dn21, const Mdouble &dn22) const
Definition: SuperQuadricParticle.cc:916
BaseInteraction
Stores information about interactions between two interactable objects; often particles but could be ...
Definition: BaseInteraction.h:60
SuperQuadricParticle::computeHessianLabFixed
SmallMatrix< 3, 3 > computeHessianLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the hessian ("second derivative") of the shape-function at the given (lab-fixed) posi...
Definition: SuperQuadricParticle.cc:527
Mdouble
double Mdouble
Definition: GeneralDefine.h:34
BaseParticle::getRadius
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:348
SmallMatrix
Data type for small dense matrix.
Definition: SmallMatrix.h:68
SuperQuadricParticle::eps1_
Mdouble eps1_
Blockiness parameters.
Definition: SuperQuadricParticle.h:292
BaseParticle::invMass_
Mdouble invMass_
Particle radius_.
Definition: BaseParticle.h:656
SuperQuadricParticle::getInitialGuessForContact
SmallVector< 4 > getInitialGuessForContact(const SuperQuadricParticle *pQuad, BaseInteraction *C) const
Get an initial guess for the contact-point between this particle and the given particle.
Definition: SuperQuadricParticle.cc:436
SuperQuadricParticle::getAxes
Vec3D getAxes() const override
Get the axes-lengths of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 ...
Definition: SuperQuadricParticle.cc:182
mathsFunc::sign
int sign(T val)
This is a sign function, it returns -1 for negative numbers, 1 for positive numbers and 0 for 0.
Definition: ExtendedMath.h:97
SuperQuadricParticle::writeDebugMessageStep1
void writeDebugMessageStep1(const SuperQuadricParticle *pQuad, const SmallVector< 4 > &contactPointPlanB) const
Definition: SuperQuadricParticle.cc:954
SuperQuadricParticle::SuperQuadricParticle
SuperQuadricParticle()
Basic Particle constructor, creates a superquadric with axes (1,1,1) and exponents (2,...
Definition: SuperQuadricParticle.cc:38
VERBOSE
LL< Log::VERBOSE > VERBOSE
Verbose information.
Definition: Logger.cc:57
ERROR
LL< Log::ERROR > ERROR
Error log level.
Definition: Logger.cc:53
BaseHandler::removeObject
virtual void removeObject(unsigned const int index)
Removes an Object from the BaseHandler.
Definition: BaseHandler.h:472
BaseInteraction::getOverlap
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
Definition: BaseInteraction.h:240
MatrixSymmetric3D::XZ
Mdouble XZ
Definition: MatrixSymmetric.h:42
SuperQuadricParticle::write
void write(std::ostream &os) const override
Write function: write this superquadric to the given output-stream, for example a restart-file.
Definition: SuperQuadricParticle.cc:100
BaseInteraction::getLagrangeMultiplier
Mdouble getLagrangeMultiplier()
Definition: BaseInteraction.h:190
BaseInteraction::setOverlap
void setOverlap(Mdouble overlap)
Set the overlap between the two interacting object.
Definition: BaseInteraction.cc:229
Log::FATAL
@ FATAL
MatrixSymmetric3D::ZZ
Mdouble ZZ
Definition: MatrixSymmetric.h:42
BaseObject::getIndex
unsigned int getIndex() const
Returns the index of the object in the handler.
Definition: BaseObject.h:118
MatrixSymmetric3D::XY
Mdouble XY
Definition: MatrixSymmetric.h:42
SuperQuadricParticle::getExponentEps1
Mdouble getExponentEps1() const override
Get the first exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter ...
Definition: SuperQuadricParticle.cc:187
SuperQuadricParticle::getContactPoint
SmallVector< 4 > getContactPoint(const SuperQuadricParticle *p, BaseInteraction *C) const
Compute the contact point between this and the given superquadric particle.
Definition: SuperQuadricParticle.cc:372
Vec3D::Y
Mdouble Y
Definition: Vector.h:65
InteractionHandler::getInteraction
BaseInteraction * getInteraction(BaseInteractable *P, BaseInteractable *I, unsigned timeStamp)
Returns the Interaction between the BaseInteractable's P and I.
Definition: InteractionHandler.cc:149
BaseInteraction::setContactPoint
void setContactPoint(Vec3D contactPoint)
Set the location of the contact point between the two interacting objects.
Definition: BaseInteraction.cc:238
constants::i
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
SuperQuadricParticle::setBoundingRadius
void setBoundingRadius()
Get the radius of the sphere that fits precisely around the particle.
Definition: SuperQuadricParticle.cc:244
BaseInteraction::setNormal
void setNormal(Vec3D normal)
Sets the normal vector between the two interacting objects.
Definition: BaseInteraction.cc:210
SuperQuadricParticle::computeShape
Mdouble computeShape(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the shape-functiion at the given (lab-fixed) position.
Definition: SuperQuadricParticle.cc:474
SuperQuadricParticle::overlapFromContactPoint
Mdouble overlapFromContactPoint(const LabFixedCoordinates &contactPoint, const LabFixedCoordinates &normal) const
Compute the distance between the contact-point and surface of this superquadric particle.
Definition: SuperQuadricParticle.cc:393
MatrixSymmetric3D::XX
Mdouble XX
The six distinctive matrix elements.
Definition: MatrixSymmetric.h:42
BaseObject::getId
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:125
SuperQuadricParticle::computeContactPoint
bool computeContactPoint(SmallVector< 4 > &contactPoint, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
Perform the actual Newton-iterations to find the contact point. Note, that it is given back as a para...
Definition: SuperQuadricParticle.cc:817
X
@ X
Definition: StatisticsVector.h:42
Quaternion::getRotationMatrix
void getRotationMatrix(SmallMatrix< 3, 3 > &A) const
Definition: Quaternion.cc:508
Quaternion::rotateBack
void rotateBack(Vec3D &position) const
Definition: Quaternion.cc:592
BaseParticle::write
void write(std::ostream &os) const override
Particle print function, which accepts an std::ostream as input.
Definition: BaseParticle.cc:330
SuperQuadricParticle::setAxes
void setAxes(const Mdouble &a1, const Mdouble &a2, const Mdouble &a3)
Set the axes-lengths to a1, a2 and a3 for this superquadric. We use the super-ellipsoid definition st...
Definition: SuperQuadricParticle.cc:164
BaseParticle::getMaxInteractionRadius
Mdouble getMaxInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e....
Definition: BaseParticle.h:362
SuperQuadricParticle::getVolume
Mdouble getVolume() const override
Definition: SuperQuadricParticle.cc:204
BaseParticle::getHandler
ParticleHandler * getHandler() const
Returns pointer to the particle's ParticleHandler.
Definition: BaseParticle.cc:661
SuperQuadricParticle
Definition: SuperQuadricParticle.h:56
Quaternion::rotate
void rotate(Vec3D &position) const
Definition: Quaternion.cc:563
BaseParticle::read
void read(std::istream &is) override
Particle read function, which accepts an std::istream as input.
Definition: BaseParticle.cc:368
MatrixSymmetric3D::YZ
Mdouble YZ
Definition: MatrixSymmetric.h:42
SuperQuadricParticle::getInteractionWithSuperQuad
BaseInteraction * getInteractionWithSuperQuad(SuperQuadricParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler)
Checks if this superquadric is in interaction with the given superquadric, and if so,...
Definition: SuperQuadricParticle.cc:323
mathsFunc::beta
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma)
Definition: ExtendedMath.cc:164
BaseParticle::isFixed
bool isFixed() const override
Is fixed Particle function. It returns whether a Particle is fixed or not, by checking its inverse Ma...
Definition: BaseParticle.h:93
BaseParticle::isSphericalParticle
virtual bool isSphericalParticle() const
Definition: BaseParticle.h:645
MatrixSymmetric3D
Implementation of a 3D symmetric matrix.
Definition: MatrixSymmetric.h:37
Vec3D::Z
Mdouble Z
Definition: Vector.h:65
SuperQuadricParticle::setAxesAndExponents
void setAxesAndExponents(const Mdouble &a1, const Mdouble &a2, const Mdouble &a3, const Mdouble &eps1, const Mdouble &eps2)
Set the geometrical properties of the superquadrics, namely the axes-lengths a1, a2 and a3,...
Definition: SuperQuadricParticle.cc:133
ParticleHandler::checkExtremaOnDelete
void checkExtremaOnDelete(BaseParticle *P)
Checks if the extrema of this ParticleHandler needs updating when a particle is deleted.
Definition: ParticleHandler.cc:1189
SuperQuadricParticle::getInteractionRadius
Mdouble getInteractionRadius(const BaseParticle *particle) const
returns the radius plus half the interactionDistance of the mixed species
Definition: SuperQuadricParticle.cc:977
BaseInteraction::setLagrangeMultiplier
void setLagrangeMultiplier(Mdouble multiplier)
Definition: BaseInteraction.h:185
BaseParticle::BaseParticle
BaseParticle()
Basic Particle constructor, creates an Particle at (0,0,0) with radius, mass and inertia equal to 1.
Definition: BaseParticle.cc:33
SuperQuadricParticle::eps2_
Mdouble eps2_
Definition: SuperQuadricParticle.h:292
SmallVector::length
Mdouble length() const
Definition: SmallVector.h:253
BaseParticle::getSumOfInteractionRadii
Mdouble getSumOfInteractionRadii(const BaseParticle *particle) const
returns the sum of the radii plus the interactionDistance
Definition: BaseParticle.h:379
SuperQuadricParticle::computeShapeGradientLabFixed
SmallVector< 3 > computeShapeGradientLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the gradient of the shape-function at the given (lab-fixed) position.
Definition: SuperQuadricParticle.cc:494
mathsFunc::isEqual
bool isEqual(Mdouble v1, Mdouble v2, Mdouble absError)
Compares the difference of two Mdouble with an absolute error, useful in UnitTests.
Definition: ExtendedMath.cc:251
mathsFunc::gamma
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
Definition: ExtendedMath.cc:137
SmallVector
Definition: SmallVector.h:62
Vec3D::getLengthSquared
Mdouble getLengthSquared() const
Calculates the squared length of this Vec3D: .
Definition: Vector.cc:184
MatrixSymmetric3D::YY
Mdouble YY
Definition: MatrixSymmetric.h:42
Z
@ Z
Definition: StatisticsVector.h:42
BaseInteractable::getPosition
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Definition: BaseInteractable.h:218
SuperQuadricParticle::getContactPointPlanB
SmallVector< 4 > getContactPointPlanB(const SuperQuadricParticle *pOther, unsigned numberOfSteps) const
If the "normal" procedure fails to find a contact point, use an alternative approach that involves st...
Definition: SuperQuadricParticle.cc:721
ParticleSpecies::computeMass
void computeMass(BaseParticle *p) const
Compute Particle mass function, which required a reference to the Species vector. It computes the Par...
Definition: ParticleSpecies.cc:166