MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 polymorfism. 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...
 
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...
 

Detailed Description

Definition at line 55 of file SuperQuadricParticle.h.

Constructor & Destructor Documentation

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.

Definition at line 38 of file SuperQuadricParticle.cc.

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

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

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 }
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
BaseParticle()
Basic Particle constructor, creates an Particle at (0,0,0) with radius, mass and inertia equal to 1...
Definition: BaseParticle.cc:33
Definition: Vector.h:49
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.

Definition at line 55 of file SuperQuadricParticle.cc.

References axes_, eps1_, and eps2_.

56  : BaseParticle(p)
57 {
58  axes_ = p.axes_;
59  eps1_ = p.eps1_;
60  eps2_ = p.eps2_;
61 }
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
BaseParticle()
Basic Particle constructor, creates an Particle at (0,0,0) with radius, mass and inertia equal to 1...
Definition: BaseParticle.cc:33
SuperQuadricParticle::SuperQuadricParticle ( const BaseParticle p)

Definition at line 64 of file SuperQuadricParticle.cc.

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

64  : BaseParticle(p)
65 {
66  Mdouble radius = p.getRadius();
67  axes_ = Vec3D(radius, radius, radius);
68  eps1_ = 1.0;
69  eps2_ = 1.0;
70 }
double Mdouble
Definition: GeneralDefine.h:34
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
BaseParticle()
Basic Particle constructor, creates an Particle at (0,0,0) with radius, mass and inertia equal to 1...
Definition: BaseParticle.cc:33
Definition: Vector.h:49
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.

Definition at line 76 of file SuperQuadricParticle.cc.

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

77 {
78  if (getHandler() != nullptr)
79  {
81  }
82  logger(DEBUG, "SuperQuadricParticle::SuperQuadricParticleParticle() of particle % finished.", getId());
83 
84 }
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:125
void checkExtremaOnDelete(BaseParticle *P)
Checks if the extrema of this ParticleHandler needs updating when a particle is deleted.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
ParticleHandler * getHandler() const
Returns pointer to the particle's ParticleHandler.

Member Function Documentation

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.

Definition at line 816 of file SuperQuadricParticle.cc.

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

Referenced by getContactPoint(), and getContactPointPlanB().

818 {
819  // Damped Newton's method: (dampingFactor 1 is undamped)
820  Mdouble dampingFactor = 1;
821  int iteration = 1;
822  const Mdouble tolerance = 1e-5;
823  const unsigned maxIterations = 100;
824 
825  logger(VERBOSE, "func to be zero value: % \n",
826  computeResidualContactDetection(contactPoint, p1, p2));
827 
828  while (computeResidualContactDetection(contactPoint, p1, p2).length() > tolerance && iteration < maxIterations)
829  {
830  logger(VERBOSE, "Iteration: %\n", iteration);
831 
832  SmallMatrix<4, 4> jacobian = getJacobianOfContactDetectionObjective(contactPoint, p1, p2);
833  SmallVector<4> func = -computeResidualContactDetection(contactPoint, p1, p2);
834  jacobian.solve(func); //note: solve is in place
835 
836  SmallVector<4> newPoint = contactPoint + dampingFactor * func;
837  const auto residueOld = computeResidualContactDetection(contactPoint, p1, p2);
838  const auto residueNew = computeResidualContactDetection(newPoint, p1, p2);
839  logger(VERBOSE, "Old value: % (%) new value: % (%)",
840  computeResidualContactDetection(contactPoint, p1, p2),
841  computeResidualContactDetection(contactPoint, p1, p2).length(),
842  computeResidualContactDetection(newPoint, p1, p2),
843  computeResidualContactDetection(newPoint, p1, p2).length());
844 
845  logger(VERBOSE, "current contact point: %, new contact point: %\n", contactPoint, newPoint);
846  logger(VERBOSE, "damping factor: %", dampingFactor);
847 
848  if (residueNew.length()
849  < residueOld.length())
850  {
851  contactPoint = newPoint;
852  dampingFactor = 1;
853  }
854  else
855  {
856  dampingFactor /= 2;
857 
858  if (dampingFactor < 1e-10)
859  {
860  return false;
861  }
862  }
863 
864  iteration++;
865  }
866  return (iteration != maxIterations);
867 }
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
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...
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...
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.
Data type for small dense matrix.
Definition: SmallMatrix.h:67
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

Definition at line 526 of file SuperQuadricParticle.cc.

References A, axes_, eps1_, eps2_, BaseInteractable::getOrientation(), BaseInteractable::getPosition(), Quaternion::getRotationMatrix(), Quaternion::rotateBack(), mathsFunc::sign(), SmallMatrix< numberOfRows, numberOfColumns >::transpose(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getCurvature(), and getJacobianOfContactDetectionObjective().

527 {
528  SmallMatrix<3, 3> hessian;
529  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
530  getOrientation().rotateBack(bodyFixedCoordinates);
531 
532  const Mdouble n1 = 2.0 / eps1_;
533  const Mdouble n2 = 2.0 / eps2_;
534 
535  const Mdouble absXoa = std::abs(bodyFixedCoordinates.X / axes_.X);
536  const Mdouble absYob = std::abs(bodyFixedCoordinates.Y / axes_.Y);
537  const Mdouble absZoc = std::abs(bodyFixedCoordinates.Z / axes_.Z);
538 
539  const Mdouble nu = std::pow(absXoa, n2) + std::pow(absYob, n2);
540  const Mdouble help1 = std::pow(nu, n1 / n2 - 1.0);
541  const Mdouble help2 = std::pow(nu, n1 / n2 - 2.0);
542 
543  hessian(0, 0) = n1 * (n2 - 1) / axes_.X / axes_.X * std::pow(absXoa, n2 - 2.0) * help1
544  + ((n1 - n2) * n1 / axes_.X / axes_.X) * std::pow(absXoa, 2. * n2 - 2.0) * help2;
545  hessian(1, 1) = n1 * (n2 - 1) / axes_.Y / axes_.Y * std::pow(absYob, n2 - 2.0) * help1
546  + ((n1 - n2) * n1 / axes_.Y / axes_.Y) * std::pow(absYob, 2. * n2 - 2.0) * help2;
547  hessian(2, 2) = n1 * (n1 - 1) / axes_.Z / axes_.Z * std::pow(absZoc, n1 - 2.0);
548  hessian(1, 0) = n1 * (n1 - n2) / axes_.X / axes_.Y * std::pow(absXoa, n2 - 1.0) * std::pow(absYob, n2 - 1.0) *
549  help2 * mathsFunc::sign(bodyFixedCoordinates.X * bodyFixedCoordinates.Y);
550  hessian(0, 1) = hessian(1, 0);
553  return A * hessian * (A.transpose());
554 }
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
void rotateBack(Vec3D &position) const
Applies the inverse rotation to a position.
Definition: Quaternion.cc:592
double Mdouble
Definition: GeneralDefine.h:34
SmallMatrix< numberOfColumns, numberOfRows > transpose() const
Definition: SmallMatrix.h:330
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:95
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
Mdouble Y
Definition: Vector.h:65
void getRotationMatrix(SmallMatrix< 3, 3 > &A) const
Retrieves the rotation matrix.
Definition: Quaternion.cc:508
Definition: Vector.h:49
Data type for small dense matrix.
Definition: SmallMatrix.h:67
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65
void SuperQuadricParticle::computeMass ( const ParticleSpecies s)
overridevirtual

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

Reimplemented from BaseParticle.

Definition at line 983 of file SuperQuadricParticle.cc.

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.

983  {
984  if (isFixed()) return;
985  if (getParticleDimensions()==3) {
986  Mdouble volume = getVolume();
987 
988  Mdouble help1 = mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_);
989  Mdouble help2 = mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0);
990  Mdouble help3 = mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_ + 1);
991  Mdouble help4 = mathsFunc::beta(1.5 * eps1_, eps1_ + 1.0);
992 
993  invMass_ = 1.0 / (volume * s.getDensity());
994  invInertia_.XX = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
995  (axes_.Y * axes_.Y * help1 * help2
996  + 4.0 * axes_.Z * axes_.Z * help3 * help4));
997  invInertia_.XY = 0.0;
998  invInertia_.XZ = 0.0;
999  invInertia_.YY = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
1000  (axes_.X * axes_.X * help1 * help2
1001  + 4.0 * axes_.Z * axes_.Z * help3 * help4));
1002  invInertia_.YZ = 0.0;
1003  invInertia_.ZZ = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
1004  (axes_.X * axes_.X + axes_.Y * axes_.Y) * help1 * help2);
1005  } else {
1006  logger(ERROR, "[SuperQuadricParticle::computeMass] SuperQuadricParticle cannot be two-dimensional (yet)");
1007  }
1008 };
Mdouble X
the vector components
Definition: Vector.h:65
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
Mdouble invMass_
Particle radius_.
Definition: BaseParticle.h:653
MatrixSymmetric3D invInertia_
Inverse Particle mass (for computation optimization)
Definition: BaseParticle.h:654
Mdouble getVolume() const override
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
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
Mdouble Y
Definition: Vector.h:65
Mdouble getDensity() const
Allows density_ to be accessed.
unsigned int getParticleDimensions() const
Returns the particle's dimensions (either 2 or 3).
Mdouble Z
Definition: Vector.h:65
Mdouble XX
The six distinctive matrix elements.
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.

Definition at line 566 of file SuperQuadricParticle.cc.

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

Referenced by computeContactPoint().

569 {
570  LabFixedCoordinates labFixedCoordinates;
571  labFixedCoordinates.X = position[0];
572  labFixedCoordinates.Y = position[1];
573  labFixedCoordinates.Z = position[2];
574  const Mdouble lagrangeMultiplier = position[3];
575 
576  // First compute the gradient of the superquadric shape function and then rotate it to the
577  // global frame of reference
578  SmallVector<3> gradThis = p1->computeShapeGradientLabFixed(labFixedCoordinates);
579 
580  SmallVector<3> gradOther = p2->computeShapeGradientLabFixed(labFixedCoordinates);
581 
582  // Eqn. (23) in Podhlozhnyuk et al, Comp. Part. Mech. 2017.
583  // where contactPoint[3] = mu
584  SmallVector<4> toBecomeZero;
585  for (unsigned int i = 0; i < 3; ++i)
586  {
587  toBecomeZero[i] = gradThis[i] + lagrangeMultiplier * lagrangeMultiplier * gradOther[i];
588  }
589  toBecomeZero[3] = p1->computeShape(labFixedCoordinates) - p2->computeShape(labFixedCoordinates);
590 
591  return toBecomeZero;
592 }
Mdouble X
the vector components
Definition: Vector.h:65
double Mdouble
Definition: GeneralDefine.h:34
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
Mdouble computeShape(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the shape-functiion at the given (lab-fixed) position.
Mdouble Y
Definition: Vector.h:65
Definition: Vector.h:49
Mdouble Z
Definition: Vector.h:65
SmallVector< 3 > computeShapeGradientLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the gradient of the shape-function at the given (lab-fixed) position.
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.

Definition at line 473 of file SuperQuadricParticle.cc.

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

Referenced by computeResidualContactDetection(), getInteractionWithSuperQuad(), isInContactWith(), and overlapFromContactPoint().

474 {
475  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
476  getOrientation().rotateBack(bodyFixedCoordinates);
477 
478  const Mdouble n1 = 2. / eps1_;
479  const Mdouble n2 = 2. / eps2_;
480 
481  return std::pow(std::pow(std::abs(bodyFixedCoordinates.X / axes_.X), n2)
482  + std::pow(std::abs(bodyFixedCoordinates.Y / axes_.Y), n2), n1 / n2)
483  + std::pow(std::abs(bodyFixedCoordinates.Z / axes_.Z), n1) - 1.0;
484 }
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
void rotateBack(Vec3D &position) const
Applies the inverse rotation to a position.
Definition: Quaternion.cc:592
double Mdouble
Definition: GeneralDefine.h:34
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
Mdouble Y
Definition: Vector.h:65
Definition: Vector.h:49
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65
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, F(X)
Todo:
Come up with good expression for when x = y = 0 and n1 < n2

Definition at line 493 of file SuperQuadricParticle.cc.

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(), and overlapFromContactPoint().

494 {
495 
496  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
497  getOrientation().rotateBack(bodyFixedCoordinates);
498 
499  const Mdouble n1 = 2. / eps1_;
500  const Mdouble n2 = 2. / eps2_;
501 
502  const Mdouble absXoa = std::abs(bodyFixedCoordinates.X / axes_.X);
503  const Mdouble absYob = std::abs(bodyFixedCoordinates.Y / axes_.Y);
504  const Mdouble absZoc = std::abs(bodyFixedCoordinates.Z / axes_.Z);
505 
506  const Mdouble nu = std::pow(absXoa, n2) + std::pow(absYob, n2);
507 
508  const Mdouble help1 = std::pow(nu, n1 / n2 - 1.0);
509 
510  Vec3D gradientBodyFixed = {
511  (n1 / axes_.X) * std::pow(absXoa, n2 - 1.0) * help1 * mathsFunc::sign(bodyFixedCoordinates.X),
512  (n1 / axes_.Y) * std::pow(absYob, n2 - 1.0) * help1 * mathsFunc::sign(bodyFixedCoordinates.Y),
513  (n1 / axes_.Z) * std::pow(absZoc, n1 - 1.0) * mathsFunc::sign(bodyFixedCoordinates.Z)};
514  //we get the gradient in terms of the lab-fixed coordinates by \nabla_X F(X) = A \nabla_x f(x).
515  getOrientation().rotate(gradientBodyFixed);
516  return {gradientBodyFixed.X, gradientBodyFixed.Y, gradientBodyFixed.Z};
517 }
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
void rotateBack(Vec3D &position) const
Applies the inverse rotation to a position.
Definition: Quaternion.cc:592
double Mdouble
Definition: GeneralDefine.h:34
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:95
void rotate(Vec3D &position) const
Applies the rotation to a position.
Definition: Quaternion.cc:563
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
Mdouble Y
Definition: Vector.h:65
Definition: Vector.h:49
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65
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.

Definition at line 90 of file SuperQuadricParticle.cc.

References SuperQuadricParticle().

91 {
92  return new SuperQuadricParticle(*this);
93 }
SuperQuadricParticle()
Basic Particle constructor, creates a superquadric with axes (1,1,1) and exponents (2...
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.

Definition at line 182 of file SuperQuadricParticle.cc.

References axes_.

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

183 {
184  return axes_;
185 }
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
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.

Definition at line 371 of file SuperQuadricParticle.cc.

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

Referenced by getInteractionWithSuperQuad(), and isInContactWith().

372 {
373  //Assuming contact between two spheres
374  SmallVector<4> approximateContactPoint = getInitialGuessForContact(p, C);
375  if (computeContactPoint(approximateContactPoint, this, p))
376  {
377  return approximateContactPoint;
378  }
379  else
380  {
381  return getContactPointPlanB(p, 4);
382  }
383 }
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...
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
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...
SmallVector< 4 > getInitialGuessForContact(const SuperQuadricParticle *pQuad, BaseInteraction *C) const
Get an initial guess for the contact-point between this particle and the given particle.
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.

Definition at line 720 of file SuperQuadricParticle.cc.

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, Vec3D::Y, and Vec3D::Z.

Referenced by getContactPoint().

721 {
722  logger(VERBOSE, "Number of iterations: %", numberOfSteps);
723  // Step 1: Compute contact point for two volume equivalent spheres
724  const Mdouble interactionRadiusThis = getInteractionRadius(pOther);
725  const Mdouble interactionRadiusOther = pOther->getInteractionRadius(this);
726  const Vec3D positionThis = getPosition();
727  const Vec3D positionOther = pOther->getPosition();
728 
729  //analytical solution for contact point between bounding spheres
730  const Vec3D initialPoint = (interactionRadiusOther * positionThis + interactionRadiusThis * positionOther) /
731  (interactionRadiusThis + interactionRadiusOther);
732 
733  SmallVector<4> contactPointPlanB;
734  // Analytical solution
735  contactPointPlanB[0] = initialPoint.X;
736  contactPointPlanB[1] = initialPoint.Y;
737  contactPointPlanB[2] = initialPoint.Z;
738  contactPointPlanB[3] = 1.0; //Need to check.
739 
740  writeDebugMessageStep1(pOther, contactPointPlanB);
741 
742  //Step 2: Compute the deltas
743  const Vec3D dAxesThis = (getAxes() - interactionRadiusThis * Vec3D(1, 1, 1)) / numberOfSteps;
744  const Mdouble dn11 = (2.0 / getExponentEps1() - 2.0) / numberOfSteps;
745  const Mdouble dn12 = (2.0 / getExponentEps2() - 2.0) / numberOfSteps;
746 
747  const Vec3D dAxesOther = (pOther->getAxes() - interactionRadiusOther * Vec3D(1, 1, 1)) / numberOfSteps;
748  const Mdouble dn21 = (2.0 / pOther->getExponentEps1() - 2.0) / numberOfSteps;
749  const Mdouble dn22 = (2.0 / pOther->getExponentEps2() - 2.0) / numberOfSteps;
750 
751  writeDebugMessageStep2(pOther, dAxesThis, dn11, dn12, dAxesOther, dn21, dn22);
752 
753  // Create two superquadrics with the above parameters for incrementally computing the contact point
756 
757  p1.setPosition(getPosition());
758  p2.setPosition(pOther->getPosition());
759 
761  p2.setOrientation(pOther->getOrientation());
762  bool success = true;
763  const unsigned maxIterations = 20;
764  for (unsigned int counter = 1; counter <= numberOfSteps; ++counter)
765  {
766 
767  // Step 3: Increment the shape and blockiness parameters
768  const Vec3D axesThis = interactionRadiusThis * Vec3D(1, 1, 1) + counter * dAxesThis;
769  const Mdouble n11 = 2.0 + counter * dn11;
770  const Mdouble n12 = 2.0 + counter * dn12;
771 
772  Vec3D axesOther = interactionRadiusOther * Vec3D(1, 1, 1) + counter * dAxesOther;
773  const Mdouble n21 = 2.0 + counter * dn21;
774  const Mdouble n22 = 2.0 + counter * dn22;
775 
776  p1.setAxesAndExponents(axesThis, 2.0 / n11, 2.0 / n12);
777  p2.setAxesAndExponents(axesOther, 2.0 / n21, 2.0 / n22);
778 
779  writeDebugMessageStep3(axesThis, n11, n12, axesOther, n21, n22);
780 
781  writeDebugMessageMiddleOfLoop(p1, p2, contactPointPlanB, counter);
782 
783  //compute the contact point of the new particles
784  success = computeContactPoint(contactPointPlanB, &p1, &p2);
785  if (!success)
786  {
787  if (numberOfSteps > maxIterations)
788  {
789  write(std::cout);
790  pOther->write(std::cout);
791  logger(ERROR, "Plan B fails even with more than 20 intermediate steps");
792  }
793  return getContactPointPlanB(pOther, 2 * numberOfSteps);
794  }
795  }
796  logger.assert(p1.getAxes().X == getAxes().X, "In getContactPointPlanB, final particle for contact-detection not "
797  "the same as original particle");
798 
799  logger(VERBOSE, "Plan B contact point: %", contactPointPlanB);
800 
801  return contactPointPlanB;
802 }
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...
Mdouble getInteractionRadius(const BaseParticle *particle) const
returns the radius plus half the interactionDistance of the mixed species
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
Vec3D getAxes() const override
Get the axes-lengths of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 ...
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...
Mdouble getExponentEps1() const override
Get the first exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter ...
void writeDebugMessageStep1(const SuperQuadricParticle *pQuad, const SmallVector< 4 > &contactPointPlanB) const
Mdouble Y
Definition: Vector.h:65
void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
void writeDebugMessageStep3(const Vec3D &axesThis, const Mdouble &n11, const Mdouble &n12, const Vec3D &axesOther, const Mdouble &n21, const Mdouble &n22) const
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...
void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
Definition: Vector.h:49
Mdouble getExponentEps2() const override
Get the second exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter...
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65
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 writeDebugMessageMiddleOfLoop(const SuperQuadricParticle &p1, const SuperQuadricParticle &p2, SmallVector< 4 > &contactPointPlanB, const unsigned int &counter) const
void write(std::ostream &os) const override
Write function: write this superquadric to the given output-stream, for example a restart-file...
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.

Definition at line 664 of file SuperQuadricParticle.cc.

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

665 {
666  SmallVector<3> gradientVec = computeShapeGradientLabFixed(labFixedCoordinates);
667  SmallMatrix<3, 1> gradient = gradientVec;
668  SmallMatrix<3, 3> hessian = computeHessianLabFixed(labFixedCoordinates);
669  Mdouble helper = ((gradient.transpose() * hessian) * gradient)(0, 0);
670  Mdouble gradientLength = gradientVec.length();
671  Mdouble helper2 = gradientLength * gradientLength * (hessian(0, 0) + hessian(1, 1) + hessian(2, 2));
672  Mdouble denominator = 2 * gradientLength * gradientLength * gradientLength;
673 
675  return (helper2 - helper) / denominator;
676 }
double Mdouble
Definition: GeneralDefine.h:34
SmallMatrix< numberOfColumns, numberOfRows > transpose() const
Definition: SmallMatrix.h:330
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
Mdouble length() const
Definition: SmallVector.h:253
Data type for small dense matrix.
Definition: SmallMatrix.h:67
SmallVector< 3 > computeShapeGradientLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the gradient of the shape-function at the given (lab-fixed) position.
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...
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.

Definition at line 187 of file SuperQuadricParticle.cc.

References eps1_.

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

188 {
189  return eps1_;
190 }
Mdouble eps1_
Blockiness parameters.
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.

Definition at line 192 of file SuperQuadricParticle.cc.

References eps2_.

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

193 {
194  return eps2_;
195 }
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.

Definition at line 435 of file SuperQuadricParticle.cc.

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

Referenced by getContactPoint().

436 {
437  SmallVector<4> approximateContactPoint;
438  if (C == nullptr || C->getOverlap() < 1e-15)
439  {
440  //this is a new interaction
441  const LabFixedCoordinates branchVector = pQuad->getPosition() - getPosition();
442  //Get the square of the distance between particle i and particle j
443  const Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
444  const Mdouble distance = sqrt(distanceSquared);
445  const LabFixedCoordinates normal = (branchVector / distance);
446  const Mdouble overlap = getSumOfInteractionRadii(pQuad) - distance;
447  const LabFixedCoordinates contactPoint = (pQuad->getPosition() -
448  (pQuad->getInteractionRadius(this) - 0.5 * overlap) * normal);
449 
450  approximateContactPoint[0] = contactPoint.X;
451  approximateContactPoint[1] = contactPoint.Y;
452  approximateContactPoint[2] = contactPoint.Z;
453  approximateContactPoint[3] = 1;
454  }
455  else
456  {
457  //this is an existing interaction
458  const LabFixedCoordinates contactPoint = C->getContactPoint();
459  const Mdouble multiplier = C->getLagrangeMultiplier();
460  approximateContactPoint[0] = contactPoint.X;
461  approximateContactPoint[1] = contactPoint.Y;
462  approximateContactPoint[2] = contactPoint.Z;
463  approximateContactPoint[3] = multiplier;
464  }
465  return approximateContactPoint;
466 }
Mdouble getInteractionRadius(const BaseParticle *particle) const
returns the radius plus half the interactionDistance of the mixed species
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
double Mdouble
Definition: GeneralDefine.h:34
Mdouble getSumOfInteractionRadii(const BaseParticle *particle) const
returns the sum of the radii plus the interactionDistance
Definition: BaseParticle.h:376
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
const Vec3D & getContactPoint() const
Gets constant reference to contact point (vector).
Mdouble getLengthSquared() const
Calculates the squared length of this Vec3D: .
Definition: Vector.cc:184
Mdouble getLagrangeMultiplier()
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
Mdouble Y
Definition: Vector.h:65
Definition: Vector.h:49
Mdouble Z
Definition: Vector.h:65
Mdouble SuperQuadricParticle::getInteractionRadius ( const BaseParticle particle) const

returns the radius plus half the interactionDistance of the mixed species

Definition at line 976 of file SuperQuadricParticle.cc.

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

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

977 {
978  const auto mixedSpecies = getSpecies()->getHandler()->getMixedObject(getSpecies(),particle->getSpecies());
979  return getRadius() + 0.5 * mixedSpecies->getInteractionDistance();
980 }
std::enable_if<!std::is_pointer< typename U::MixedSpeciesType >::value, typename U::MixedSpeciesType * >::type getMixedObject(const U *S, const U *T)
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
SpeciesHandler * getHandler() const
Returns the pointer to the handler to which this species belongs.
Definition: BaseSpecies.cc:99
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
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.

Definition at line 286 of file SuperQuadricParticle.cc.

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

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 }
std::enable_if<!std::is_pointer< typename U::MixedSpeciesType >::value, typename U::MixedSpeciesType * >::type getMixedObject(const U *S, const U *T)
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
double Mdouble
Definition: GeneralDefine.h:34
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
SuperQuadricParticle()
Basic Particle constructor, creates a superquadric with axes (1,1,1) and exponents (2...
Mdouble getMaxInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles)
Definition: BaseParticle.h:359
BaseInteraction * getInteractionWithSuperQuad(SuperQuadricParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler)
Checks if this superquadric is in interaction with the given superquadric, and if so...
Mdouble getLengthSquared() const
Calculates the squared length of this Vec3D: .
Definition: Vector.cc:184
Stores information about interactions between two interactable objects; often particles but could be ...
SpeciesHandler * getHandler() const
Returns the pointer to the handler to which this species belongs.
Definition: BaseSpecies.cc:99
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: Vector.h:49
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'

Definition at line 323 of file SuperQuadricParticle.cc.

References computeShape(), computeShapeGradientLabFixed(), getContactPoint(), InteractionHandler::getInteraction(), BaseInteraction::getOverlap(), BaseInteractable::getPosition(), overlapFromContactPoint(), BaseInteraction::setContactPoint(), BaseInteraction::setDistance(), BaseInteraction::setLagrangeMultiplier(), BaseInteraction::setNormal(), BaseInteraction::setOverlap(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getInteractionWith().

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  return nullptr;
339  }
340 
341  C->setContactPoint(contactPoint);
342  C->setLagrangeMultiplier(approximateContactPoint[3]);
343 
344  SmallVector<3> gradThis = computeShapeGradientLabFixed(contactPoint).getNormalised();
345  LabFixedCoordinates normal;
346  normal.X = gradThis[0];
347  normal.Y = gradThis[1];
348  normal.Z = gradThis[2];
349  C->setNormal(normal);
350 
351  const Mdouble alphaI = overlapFromContactPoint(contactPoint, normal);
352  const Mdouble alphaJ = p->overlapFromContactPoint(contactPoint, -normal);
353  C->setOverlap(alphaI + alphaJ);
355  C->setDistance((getPosition() - p->getPosition()).getLength() - C->getOverlap());
356 
357  return C;
358 }
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
void setNormal(Vec3D normal)
Sets the normal vector between the two interacting objects.
Mdouble X
the vector components
Definition: Vector.h:65
void setOverlap(Mdouble overlap)
Set the overlap between the two interacting object.
double Mdouble
Definition: GeneralDefine.h:34
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
void setContactPoint(Vec3D contactPoint)
Set the location of the contact point between the two interacting objects.
void setLagrangeMultiplier(Mdouble multiplier)
void setDistance(Mdouble distance)
Sets the interaction distance between the two interacting objects.
SmallVector< 4 > getContactPoint(const SuperQuadricParticle *p, BaseInteraction *C) const
Compute the contact point between this and the given superquadric particle.
Stores information about interactions between two interactable objects; often particles but could be ...
BaseInteraction * getInteraction(BaseInteractable *P, BaseInteractable *I, unsigned timeStamp)
Returns the Interaction between the BaseInteractable's P and I.
Mdouble computeShape(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the shape-functiion at the given (lab-fixed) position.
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
Mdouble Y
Definition: Vector.h:65
Definition: Vector.h:49
Mdouble overlapFromContactPoint(const LabFixedCoordinates &contactPoint, const LabFixedCoordinates &normal) const
Compute the distance between the contact-point and surface of this superquadric particle.
Mdouble Z
Definition: Vector.h:65
SmallVector< 3 > computeShapeGradientLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the gradient of the shape-function at the given (lab-fixed) position.
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]The contact point where this Jacobian should be evaluated
[in]p1First particle in the contact for which this Jacobian should be evaluated
[in]p2Second particle in the contact for which this Jacobian should be evaluated
Returns
A 4x4 matrix with the Jacobian of computeResidualContactDetection.

Definition at line 623 of file SuperQuadricParticle.cc.

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

Referenced by computeContactPoint().

626 {
627  LabFixedCoordinates labFixedCoordinates;
628  labFixedCoordinates.X = contactPoint[0];
629  labFixedCoordinates.Y = contactPoint[1];
630  labFixedCoordinates.Z = contactPoint[2];
631  const Mdouble lagrangeMultiplier = contactPoint[3];
632 
633  SmallVector<3> gradThis = p1->computeShapeGradientLabFixed(labFixedCoordinates);
634  SmallVector<3> gradOther = p2->computeShapeGradientLabFixed(labFixedCoordinates);
635 
636  SmallMatrix<3, 3> hessianThis = p1->computeHessianLabFixed(labFixedCoordinates);
637  SmallMatrix<3, 3> hessianOther = p2->computeHessianLabFixed(labFixedCoordinates);
638 
639  SmallMatrix<3, 3> upperLeftCorner = hessianThis + lagrangeMultiplier * lagrangeMultiplier * hessianOther;
640 
641  SmallMatrix<4, 4> jacobian;
642  for (unsigned int i = 0; i < 3; ++i)
643  {
644  for (unsigned int j = 0; j < 3; ++j)
645  {
646  jacobian(i, j) = upperLeftCorner(i, j);
647  }
648  jacobian(i, 3) = 2 * lagrangeMultiplier * gradOther[i];
649  }
650 
651  for (unsigned int j = 0; j < 3; ++j)
652  {
653  jacobian(3, j) = gradThis[j] - gradOther[j];
654  }
655 
656  return jacobian;
657 }
Mdouble X
the vector components
Definition: Vector.h:65
double Mdouble
Definition: GeneralDefine.h:34
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
Mdouble Y
Definition: Vector.h:65
Definition: Vector.h:49
Data type for small dense matrix.
Definition: SmallMatrix.h:67
Mdouble Z
Definition: Vector.h:65
SmallVector< 3 > computeShapeGradientLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the gradient of the shape-function at the given (lab-fixed) position.
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...
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.

Definition at line 112 of file SuperQuadricParticle.cc.

113 {
114  return "SuperQuadricParticle";
115 }
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.

Definition at line 204 of file SuperQuadricParticle.cc.

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

Referenced by computeMass().

205 {
206  logger.assert(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 }
Mdouble X
the vector components
Definition: Vector.h:65
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
ParticleHandler * getHandler() const
Returns pointer to the particle's ParticleHandler.
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
Mdouble Y
Definition: Vector.h:65
Mdouble Z
Definition: Vector.h:65
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.

Definition at line 685 of file SuperQuadricParticle.cc.

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

686 {
687  SmallVector<4> approximateContactPoint;
688 
689  if (p->isSphericalParticle())
690  {
692  pQuad->setAxes(p->getRadius(), p->getRadius(), p->getRadius());
693  approximateContactPoint = getContactPoint(pQuad, nullptr);
694  delete pQuad;
695  } else {
696  const SuperQuadricParticle* pQuad = static_cast<const SuperQuadricParticle*>(p);
697  approximateContactPoint = getContactPoint(pQuad, nullptr);
698  }
699 
700  //set the new contact point:
701  LabFixedCoordinates contactPoint;
702  contactPoint.X = approximateContactPoint[0];
703  contactPoint.Y = approximateContactPoint[1];
704  contactPoint.Z = approximateContactPoint[2];
705 
706  //if the minimum is positive, there is no contact: return false
707  return (computeShape(contactPoint) < 0);
708 
709 }
Mdouble X
the vector components
Definition: Vector.h:65
virtual bool isSphericalParticle() const
Definition: BaseParticle.h:642
SuperQuadricParticle()
Basic Particle constructor, creates a superquadric with axes (1,1,1) and exponents (2...
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
SmallVector< 4 > getContactPoint(const SuperQuadricParticle *p, BaseInteraction *C) const
Compute the contact point between this and the given superquadric particle.
Mdouble computeShape(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the shape-functiion at the given (lab-fixed) position.
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:345
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...
Mdouble Y
Definition: Vector.h:65
Definition: Vector.h:49
Mdouble Z
Definition: Vector.h:65
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

Definition at line 392 of file SuperQuadricParticle.cc.

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

Referenced by getInteractionWithSuperQuad().

394 {
395  Mdouble alphaI = 0;
396  LabFixedCoordinates xEdge = contactPoint + alphaI * normal;
397  Mdouble dampingCoefficient = 1;
398  while (std::abs(computeShape(xEdge)) > 1e-10)
399  {
400  SmallVector<3> gradientShape = computeShapeGradientLabFixed(xEdge);
401  LabFixedCoordinates gradient;
402  gradient.X = gradientShape(0);
403  gradient.Y = gradientShape(1);
404  gradient.Z = gradientShape(2);
405  Mdouble newValue = alphaI - dampingCoefficient * computeShape(xEdge) / Vec3D::dot(gradient, normal);
406  const LabFixedCoordinates xEdgeNew = contactPoint + newValue * normal;
407 
408  if (std::abs(computeShape(xEdgeNew)) < std::abs(computeShape(xEdge)))
409  {
410  alphaI = newValue;
411  xEdge = xEdgeNew;
412  dampingCoefficient = 1;
413  }
414  else
415  {
416  dampingCoefficient /= 2;
417  if (dampingCoefficient < 1e-10)
418  {
419  logger(ERROR, "overlap detection does not converge");
420  }
421  }
422  }
423  return alphaI;
424 }
Mdouble X
the vector components
Definition: Vector.h:65
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:76
Mdouble computeShape(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the shape-functiion at the given (lab-fixed) position.
Mdouble Y
Definition: Vector.h:65
Definition: Vector.h:49
Mdouble Z
Definition: Vector.h:65
SmallVector< 3 > computeShapeGradientLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the gradient of the shape-function at the given (lab-fixed) position.
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.

Definition at line 122 of file SuperQuadricParticle.cc.

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

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 }
void read(std::istream &is) override
Particle read function, which accepts an std::istream as input.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
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.

Definition at line 164 of file SuperQuadricParticle.cc.

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

165 {
166  setAxes({a1,a2,a3});
167 }
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...
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.

Definition at line 154 of file SuperQuadricParticle.cc.

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

155 {
156  axes_ = axes;
158  if (getSpecies() != nullptr)
159  {
160  getSpecies()->computeMass(this);
161  }
162 }
void setBoundingRadius()
Get the radius of the sphere that fits precisely around the particle.
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
void computeMass(BaseParticle *p) const
Compute Particle mass function, which required a reference to the Species vector. It computes the Par...
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
void SuperQuadricParticle::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.

Definition at line 133 of file SuperQuadricParticle.cc.

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

Referenced by getContactPointPlanB().

135 {
136  eps1_ = eps1;
137  eps2_ = eps2;
138  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
139  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
140 
141  setAxes(a1,a2,a3);
142 }
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Mdouble eps1_
Blockiness parameters.
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...
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.

Definition at line 144 of file SuperQuadricParticle.cc.

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

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 }
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Mdouble eps1_
Blockiness parameters.
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...
void SuperQuadricParticle::setBoundingRadius ( )
private

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

Todo:
Currently only implemented for ellipsoids

Definition at line 244 of file SuperQuadricParticle.cc.

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().

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 }
Mdouble X
the vector components
Definition: Vector.h:65
double Mdouble
Definition: GeneralDefine.h:34
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species) ...
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
bool isEqual(Mdouble v1, Mdouble v2, Mdouble absError)
Compares the difference of two Mdouble with an absolute error, useful in UnitTests.
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
Mdouble Y
Definition: Vector.h:65
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
Mdouble Z
Definition: Vector.h:65
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.

Definition at line 169 of file SuperQuadricParticle.cc.

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

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 }
void setBoundingRadius()
Get the radius of the sphere that fits precisely around the particle.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
void computeMass(BaseParticle *p) const
Compute Particle mass function, which required a reference to the Species vector. It computes the Par...
Mdouble eps1_
Blockiness parameters.
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.

Definition at line 216 of file SuperQuadricParticle.cc.

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.

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 }
Mdouble X
the vector components
Definition: Vector.h:65
virtual void setInertia()
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Mdouble eps1_
Blockiness parameters.
Mdouble Y
Definition: Vector.h:65
Mdouble getDensity() const
Allows density_ to be accessed.
Mdouble Z
Definition: Vector.h:65
Implementation of a 3D symmetric matrix.
Mdouble XX
The six distinctive matrix elements.
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.

Definition at line 239 of file SuperQuadricParticle.cc.

References ERROR, and logger.

240 {
241  logger(ERROR,"This function should not be used");
242 }
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
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.

Definition at line 100 of file SuperQuadricParticle.cc.

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

Referenced by getContactPointPlanB().

101 {
103  os << " axes " << axes_
104  << " exp1 " << eps1_
105  << " exp2 " << eps2_;
106 }
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
void write(std::ostream &os) const override
Particle print function, which accepts an std::ostream as input.
Mdouble eps1_
Blockiness parameters.
void SuperQuadricParticle::writeDebugMessageMiddleOfLoop ( const SuperQuadricParticle p1,
const SuperQuadricParticle p2,
SmallVector< 4 > &  contactPointPlanB,
const unsigned int &  counter 
) const

Definition at line 869 of file SuperQuadricParticle.cc.

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

Referenced by getContactPointPlanB().

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

Definition at line 953 of file SuperQuadricParticle.cc.

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

Referenced by getContactPointPlanB().

954 {
955  logger(VERBOSE, "---------------------");
956  logger(VERBOSE, "STEP 1");
957  logger(VERBOSE, "---------------------\n");
958 
959  logger(VERBOSE, "Position of particle 1: % ", getPosition());
960  logger(VERBOSE, "Position of particle 2 (pQuad): % ", pQuad->getPosition());
961  logger(VERBOSE, " ");
962 
963  logger(VERBOSE, "Radius of particle 1: % ", getInteractionRadius(pQuad));
964  logger(VERBOSE, "Radius of particle 2 (pQuad): % \n", pQuad->getInteractionRadius(this));
965 
966  logger(VERBOSE, "Orientation of particle 1: % ", getOrientation());
967  logger(VERBOSE, "Orientation of particle 2 (pQuad): % ", pQuad->getOrientation());
968  logger(VERBOSE, " ");
969 
970  logger(VERBOSE, "Particle 1 axes: %", getAxes());
971  logger(VERBOSE, "Particle 2 axes (pQuad): % \n", pQuad->getAxes());
972 
973  logger(VERBOSE, "Analytical solution for two equivalent spheres in contact: % \n", contactPointPlanB);
974 }
Mdouble getInteractionRadius(const BaseParticle *particle) const
returns the radius plus half the interactionDistance of the mixed species
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Vec3D getAxes() const override
Get the axes-lengths of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 ...
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
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

Definition at line 915 of file SuperQuadricParticle.cc.

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

Referenced by getContactPointPlanB().

918 {
919  logger(VERBOSE, "---------------------");
920  logger(VERBOSE, "STEP 2");
921  logger(VERBOSE, "---------------------");
922  logger(VERBOSE, " ");
923 
924  logger(VERBOSE, "Particle 1 x-scale: %", pQuad->getAxes().X);
925  logger(VERBOSE, "Particle 1 y-scale: %", pQuad->getAxes().Y);
926  logger(VERBOSE, "Particle 1 z-scale: %", pQuad->getAxes().Z);
927  logger(VERBOSE, "Particle 1 n1: %", 2. / pQuad->getExponentEps1());
928  logger(VERBOSE, "Particle 1 n2: %", 2. / pQuad->getExponentEps2());
929  logger(VERBOSE, " ");
930 
931  logger(VERBOSE, "Particle 1 x-scale increment: %", dAxesThis.X);
932  logger(VERBOSE, "Particle 1 y-scale increment: %", dAxesThis.Y);
933  logger(VERBOSE, "Particle 1 z-scale increment: %", dAxesThis.Z);
934  logger(VERBOSE, "Particle 1 n1 increment: %", dn11);
935  logger(VERBOSE, "Particle 1 n2 increment: %", dn12);
936  logger(VERBOSE, " ");
937 
938  logger(VERBOSE, "Particle 2 x-scale: %", getAxes().X);
939  logger(VERBOSE, "Particle 2 y-scale: %", getAxes().Y);
940  logger(VERBOSE, "Particle 2 z-scale: %", getAxes().Z);
941  logger(VERBOSE, "Particle 2 n1: %", 2. / getExponentEps1());
942  logger(VERBOSE, "Particle 2 n2: %", 2. / getExponentEps2());
943  logger(VERBOSE, " ");
944 
945  logger(VERBOSE, "Particle 2 x-scale increment: %", dAxesOther.X);
946  logger(VERBOSE, "Particle 2 y-scale increment: %", dAxesOther.Y);
947  logger(VERBOSE, "Particle 2 z-scale increment: %", dAxesOther.Z);
948  logger(VERBOSE, "Particle 2 n1 increment: %", dn21);
949  logger(VERBOSE, "Particle 2 n2 increment: %", dn22);
950  logger(VERBOSE, " ");
951 }
Mdouble X
the vector components
Definition: Vector.h:65
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Vec3D getAxes() const override
Get the axes-lengths of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 ...
Mdouble getExponentEps1() const override
Get the first exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter ...
Mdouble Y
Definition: Vector.h:65
Mdouble getExponentEps2() const override
Get the second exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter...
Mdouble Z
Definition: Vector.h:65
void SuperQuadricParticle::writeDebugMessageStep3 ( const Vec3D axesThis,
const Mdouble n11,
const Mdouble n12,
const Vec3D axesOther,
const Mdouble n21,
const Mdouble n22 
) const

Definition at line 891 of file SuperQuadricParticle.cc.

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

Referenced by getContactPointPlanB().

893 {
894  logger(VERBOSE, "-----------------------------------");
895  logger(VERBOSE, "STEP 3: First increment");
896  logger(VERBOSE, "-----------------------------------");
897  logger(VERBOSE, " ");
898 
899  logger(VERBOSE, "Particle 1 x-scale after increment: %", axesThis.X);
900  logger(VERBOSE, "Particle 1 y-scale after increment: %", axesThis.Y);
901  logger(VERBOSE, "Particle 1 z-scale after increment: %", axesThis.Z);
902  logger(VERBOSE, "Particle 1 n1 after increment: %", n11);
903  logger(VERBOSE, "Particle 1 n2 after increment: %", n12);
904  logger(VERBOSE, " ");
905 
906  logger(VERBOSE, "Particle 2 x-scale after increment: %", axesOther.X);
907  logger(VERBOSE, "Particle 2 y-scale after increment: %", axesOther.Y);
908  logger(VERBOSE, "Particle 2 z-scale after increment: %", axesOther.Z);
909  logger(VERBOSE, "Particle 2 n1 after increment: %", n21);
910  logger(VERBOSE, "Particle 2 n2 after increment: %", n22);
911  logger(VERBOSE, " ");
912 }
Mdouble X
the vector components
Definition: Vector.h:65
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Mdouble Y
Definition: Vector.h:65
Mdouble Z
Definition: Vector.h:65

Member Data Documentation

Vec3D SuperQuadricParticle::axes_
private
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.

Definition at line 292 of file SuperQuadricParticle.h.

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


The documentation for this class was generated from the following files: