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)
 Base class copy constructor. Creates a Superquadric particle from a NonSphericalParticle. More...
 
 ~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 "SuperQuadric". 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 NonSphericalParticle
 NonSphericalParticle ()=default
 
 NonSphericalParticle (const NonSphericalParticle &p)=default
 
 NonSphericalParticle (const BaseParticle &p)
 Base class copy constructor. Creates a NonSphericalParticle particle from a BaseParticle. More...
 
 ~NonSphericalParticle () override=default
 
bool isSphericalParticle () const override
 
virtual Mdouble getKineticEnergy () const override
 
virtual Mdouble getRotationalEnergy () const override
 Calculates the particle's rotational kinetic energy. More...
 
virtual Mdouble getSurfaceArea () const override
 
- Public Member Functions inherited from BaseParticle
 BaseParticle ()
 Basic Particle constructor, creates an Particle at (0,0,0) with radius, mass and inertia equal to 1. More...
 
 BaseParticle (const BaseParticle &p)
 Particle copy constructor, which accepts as input a reference to a Particle. It creates a copy of this Particle and all it's information. Usually it is better to use the copy() function for polymorphism. More...
 
 BaseParticle (const ParticleSpecies *s)
 
 ~BaseParticle () override
 Particle destructor, needs to be implemented and checked if it removes tangential spring information. More...
 
void fixParticle ()
 Fix Particle function. It fixes a Particle by setting its inverse mass and inertia and velocities to zero. More...
 
bool isFixed () const override
 Is fixed Particle function. It returns whether a Particle is fixed or not, by checking its inverse Mass. More...
 
bool isMPIParticle () const
 Indicates if this particle is a ghost in the MPI domain. More...
 
void setMPIParticle (bool flag)
 Flags the mpi particle status. More...
 
bool isInMPIDomain ()
 Indicates if the particle is in the communication zone of the mpi domain. More...
 
void setInMPIDomain (bool flag)
 Flags the status of the particle if wether it is in the communication zone or not. More...
 
bool isInPeriodicDomain () const
 Indicates if the particle is in the periodic boundary communication zone. More...
 
void setInPeriodicDomain (bool flag)
 Flags the status of the particle whether 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 setTimeStamp (unsigned timeStamp)
 
unsigned getTimeStamp () const
 
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 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 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...
 
MERCURYDPM_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 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 setHandler (ParticleHandler *handler)
 Sets the pointer to the particle's ParticleHandler. More...
 
ParticleHandlergetHandler () const
 Returns pointer to the particle's ParticleHandler. More...
 
virtual void integrateBeforeForceComputation (double time, double timeStep)
 First step of Velocity Verlet integration. More...
 
virtual 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...
 
MERCURYDPM_DEPRECATED void setIndSpecies (unsigned int indSpecies) override
 
virtual void setSpecies (const ParticleSpecies *species)
 
virtual void actionsBeforeTimeStep ()
 
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 ()
 
const HGridCellgetHGridCell () const
 
BaseParticlegetClump () const
 
bool isClump () const
 Checks if particle is a clump (container) More...
 
bool isPebble () const
 Checks if particle is a pebble (belongs to a clump) More...
 
virtual Vec3D getCenterOfMass ()
 
virtual void actionsAfterAddObject ()
 
- 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...
 
virtual 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...
 
virtual 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...
 
virtual 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...
 
virtual bool isFaceContact (const Vec3D &normal) const
 
- 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

- Public Attributes inherited from BaseParticle
Mdouble radius_
 
Mdouble invMass_
 Particle radius_. More...
 
MatrixSymmetric3D invInertia_
 Inverse Particle mass (for computation optimization) More...
 
BaseParticleclumpParticle_
 Function that updates necessary quantities of a clump particle after adding a pebble. More...
 
bool isPebble_
 pointer to a clump particle (for a pebble) More...
 
bool isClump_
 The particle is pebble. More...
 

Constructor & Destructor Documentation

◆ SuperQuadricParticle() [1/3]

SuperQuadricParticle::SuperQuadricParticle ( )

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

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

19 {
20  axes_ = Vec3D(1.0, 1.0, 1.0);
21  eps1_ = 1.0;
22  eps2_ = 1.0;
23  logger(DEBUG, "SuperQuadricParticle::SuperQuadricParticle() finished");
24 }
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
@ DEBUG
NonSphericalParticle()=default
Vec3D axes_
Lengths of principal axes (a1, a2, a3).
Definition: SuperQuadricParticle.h:281
Mdouble eps1_
Blockiness parameters.
Definition: SuperQuadricParticle.h:276
Mdouble eps2_
Definition: SuperQuadricParticle.h:276
Definition: Kernel/Math/Vector.h:30

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

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

◆ SuperQuadricParticle() [2/3]

SuperQuadricParticle::SuperQuadricParticle ( const SuperQuadricParticle p)

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

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

Parameters
[in,out]pReference to the SuperQuad this one should become a copy of.
36 {
37  axes_ = p.axes_;
38  eps1_ = p.eps1_;
39  eps2_ = p.eps2_;
40 }
float * p
Definition: Tutorial_Map_using.cpp:9

References axes_, eps1_, eps2_, and p.

◆ SuperQuadricParticle() [3/3]

SuperQuadricParticle::SuperQuadricParticle ( const BaseParticle p)

Base class copy constructor. Creates a Superquadric particle from a NonSphericalParticle.

44 {
45  Mdouble radius = p.getRadius();
47  eps1_ = 1.0;
48  eps2_ = 1.0;
49 }
radius
Definition: UniformPSDSelfTest.py:15

References axes_, eps1_, eps2_, p, and UniformPSDSelfTest::radius.

◆ ~SuperQuadricParticle()

SuperQuadricParticle::~SuperQuadricParticle ( )
override

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

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

56 {
57  if (getHandler() != nullptr)
58  {
60  }
61  logger(DEBUG, "SuperQuadricParticle::SuperQuadricParticleParticle() of particle % finished.", getId());
62 
63 }
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:104
ParticleHandler * getHandler() const
Returns pointer to the particle's ParticleHandler.
Definition: BaseParticle.cc:650
void checkExtremaOnDelete(BaseParticle *P)
Checks if the extrema of this ParticleHandler needs updating when a particle is deleted.
Definition: ParticleHandler.cc:1197

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

Member Function Documentation

◆ computeContactPoint()

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

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

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

Parameters
[in|out]contactPoint The contact point we are looking for. The input is an approximate contact point that is obtained by another method, and we update this point until we reach the defined tolerance.
[in]p1The first particle of the contact we are looking for
[in]p2The second particle of the contact we are looking for
Returns
Boolean for whether or not the contact-detection was successful.
797 {
798  // Damped Newton's method: (dampingFactor 1 is undamped)
799  Mdouble dampingFactor = 1;
800  int iteration = 1;
801  const Mdouble tolerance = 1e-5;
802  const unsigned maxIterations = 100;
803 
804  logger(VERBOSE, "func to be zero value: % \n",
805  computeResidualContactDetection(contactPoint, p1, p2));
806 
807  while (computeResidualContactDetection(contactPoint, p1, p2).length() > tolerance && iteration < maxIterations)
808  {
809  logger(VERBOSE, "Iteration: %\n", iteration);
810 
811  SmallMatrix<4, 4> jacobian = getJacobianOfContactDetectionObjective(contactPoint, p1, p2);
813  jacobian.solve(func); //note: solve is in place
814 
815  SmallVector<4> newPoint = contactPoint + dampingFactor * func;
816  const auto residueOld = computeResidualContactDetection(contactPoint, p1, p2);
817  const auto residueNew = computeResidualContactDetection(newPoint, p1, p2);
818  logger(VERBOSE, "Old value: % (%) new value: % (%)",
819  computeResidualContactDetection(contactPoint, p1, p2),
820  computeResidualContactDetection(contactPoint, p1, p2).length(),
821  computeResidualContactDetection(newPoint, p1, p2),
822  computeResidualContactDetection(newPoint, p1, p2).length());
823 
824  logger(VERBOSE, "current contact point: %, new contact point: %\n", contactPoint, newPoint);
825  logger(VERBOSE, "damping factor: %", dampingFactor);
826 
827  if (residueNew.length()
828  < residueOld.length())
829  {
830  contactPoint = newPoint;
831  dampingFactor = 1;
832  }
833  else
834  {
835  dampingFactor /= 2;
836 
837  if (dampingFactor < 1e-10)
838  {
839  return false;
840  }
841  }
842 
843  iteration++;
844  }
845  return (iteration != maxIterations);
846 }
Array< double, 1, 3 > e(1./3., 0.5, 2.)
@ VERBOSE
Vector3f p1
Definition: MatrixBase_all.cpp:2
Data type for small dense matrix.
Definition: SmallMatrix.h:48
void solve(SmallMatrix< numberOfRows, numberOfRightHandSideColumns > &B) const
solves Ax=B where A is the current matrix and B is passed in. The result is returned in B.
Definition: SmallMatrix_impl.h:309
Definition: SmallVector.h:42
SmallVector< 4 > computeResidualContactDetection(const SmallVector< 4 > &position, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
Objective function for contact detection between the two given superquadrics. See Podlozhyuk et al....
Definition: SuperQuadricParticle.cc:545
SmallMatrix< 4, 4 > getJacobianOfContactDetectionObjective(const SmallVector< 4 > &contactPoint, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
Compute and return the derivative of functionThatShouldBecomeZeroForContactDetection,...
Definition: SuperQuadricParticle.cc:602
func(actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha)
Definition: benchGeometry.cpp:21

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

Referenced by getContactPoint(), and getContactPointPlanB().

◆ computeHessianLabFixed()

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

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

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

Returns
The hessian of the shape function at the given (lab-fixed) coordinates. Note, that this is the hessian to the lab-fixed coordinates, H_X (F)(X)
Todo:
Come up with good expression for when x = y = 0 and n1 < n2
506 {
507  SmallMatrix<3, 3> hessian;
508  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
509  getOrientation().rotateBack(bodyFixedCoordinates);
510 
511  const Mdouble n1 = 2.0 / eps1_;
512  const Mdouble n2 = 2.0 / eps2_;
513 
514  const Mdouble absXoa = std::abs(bodyFixedCoordinates.X / axes_.X);
515  const Mdouble absYob = std::abs(bodyFixedCoordinates.Y / axes_.Y);
516  const Mdouble absZoc = std::abs(bodyFixedCoordinates.Z / axes_.Z);
517 
518  const Mdouble nu = std::pow(absXoa, n2) + std::pow(absYob, n2);
519  const Mdouble help1 = std::pow(nu, n1 / n2 - 1.0);
520  const Mdouble help2 = std::pow(nu, n1 / n2 - 2.0);
521 
522  hessian(0, 0) = n1 * (n2 - 1) / axes_.X / axes_.X * std::pow(absXoa, n2 - 2.0) * help1
523  + ((n1 - n2) * n1 / axes_.X / axes_.X) * std::pow(absXoa, 2. * n2 - 2.0) * help2;
524  hessian(1, 1) = n1 * (n2 - 1) / axes_.Y / axes_.Y * std::pow(absYob, n2 - 2.0) * help1
525  + ((n1 - n2) * n1 / axes_.Y / axes_.Y) * std::pow(absYob, 2. * n2 - 2.0) * help2;
526  hessian(2, 2) = n1 * (n1 - 1) / axes_.Z / axes_.Z * std::pow(absZoc, n1 - 2.0);
527  hessian(1, 0) = n1 * (n1 - n2) / axes_.X / axes_.Y * std::pow(absXoa, n2 - 1.0) * std::pow(absYob, n2 - 1.0) *
528  help2 * mathsFunc::sign(bodyFixedCoordinates.X * bodyFixedCoordinates.Y);
529  hessian(0, 1) = hessian(1, 0);
532  return A * hessian * (A.transpose());
533 }
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:47
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Definition: BaseInteractable.h:209
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Definition: BaseInteractable.h:197
The matrix class, also used for vectors and row-vectors.
Definition: Eigen/Eigen/src/Core/Matrix.h:186
void rotateBack(Vec3D &position) const
Definition: Quaternion.cc:597
void getRotationMatrix(SmallMatrix< 3, 3 > &A) const
Definition: Quaternion.cc:495
Mdouble Y
Definition: Kernel/Math/Vector.h:45
Mdouble Z
Definition: Kernel/Math/Vector.h:45
Mdouble X
the vector components
Definition: Kernel/Math/Vector.h:45
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625
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:77

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

Referenced by getCurvature(), and getJacobianOfContactDetectionObjective().

◆ computeMass()

void SuperQuadricParticle::computeMass ( const ParticleSpecies s)
overridevirtual

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

Reimplemented from BaseParticle.

962  {
963  if (isFixed()) return;
964  if (getParticleDimensions()==3) {
965  Mdouble volume = getVolume();
966 
967  Mdouble help1 = mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_);
968  Mdouble help2 = mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0);
969  Mdouble help3 = mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_ + 1);
970  Mdouble help4 = mathsFunc::beta(1.5 * eps1_, eps1_ + 1.0);
971 
972  invMass_ = 1.0 / (volume * s.getDensity());
973  invInertia_.XX = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
974  (axes_.Y * axes_.Y * help1 * help2
975  + 4.0 * axes_.Z * axes_.Z * help3 * help4));
976  invInertia_.XY = 0.0;
977  invInertia_.XZ = 0.0;
978  invInertia_.YY = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
979  (axes_.X * axes_.X * help1 * help2
980  + 4.0 * axes_.Z * axes_.Z * help3 * help4));
981  invInertia_.YZ = 0.0;
982  invInertia_.ZZ = 1.0 / (s.getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
983  (axes_.X * axes_.X + axes_.Y * axes_.Y) * help1 * help2);
984  } else {
985  logger(ERROR, "[SuperQuadricParticle::computeMass] SuperQuadricParticle cannot be two-dimensional (yet)");
986  }
987 };
@ ERROR
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:72
unsigned int getParticleDimensions() const
Returns the particle's dimensions (either 2 or 3).
Definition: BaseParticle.cc:764
MatrixSymmetric3D invInertia_
Inverse Particle mass (for computation optimization)
Definition: BaseParticle.h:652
Mdouble invMass_
Particle radius_.
Definition: BaseParticle.h:651
Mdouble ZZ
Definition: MatrixSymmetric.h:21
Mdouble YY
Definition: MatrixSymmetric.h:21
Mdouble XZ
Definition: MatrixSymmetric.h:21
Mdouble XY
Definition: MatrixSymmetric.h:21
Mdouble XX
The six distinctive matrix elements.
Definition: MatrixSymmetric.h:21
Mdouble YZ
Definition: MatrixSymmetric.h:21
Mdouble getVolume() const override
Definition: SuperQuadricParticle.cc:183
RealScalar s
Definition: level1_cplx_impl.h:130
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma)
Definition: ExtendedMath.cc:143

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

◆ computeResidualContactDetection()

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

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

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

Parameters
[in]positionCurrent guess for the contact-point, where the expression above should be evaluated
[in]p1First particle for which we are looking for a contact with.
[in]p2Second particle for which we are looking for a contact with.
Returns
Residual of the objective function.
548 {
549  LabFixedCoordinates labFixedCoordinates;
550  labFixedCoordinates.X = position[0];
551  labFixedCoordinates.Y = position[1];
552  labFixedCoordinates.Z = position[2];
553  const Mdouble lagrangeMultiplier = position[3];
554 
555  // First compute the gradient of the superquadric shape function and then rotate it to the
556  // global frame of reference
557  SmallVector<3> gradThis = p1->computeShapeGradientLabFixed(labFixedCoordinates);
558 
559  SmallVector<3> gradOther = p2->computeShapeGradientLabFixed(labFixedCoordinates);
560 
561  // Eqn. (23) in Podhlozhnyuk et al, Comp. Part. Mech. 2017.
562  // where contactPoint[3] = mu
563  SmallVector<4> toBecomeZero;
564  for (unsigned int i = 0; i < 3; ++i)
565  {
566  toBecomeZero[i] = gradThis[i] + lagrangeMultiplier * lagrangeMultiplier * gradOther[i];
567  }
568  toBecomeZero[3] = p1->computeShape(labFixedCoordinates) - p2->computeShape(labFixedCoordinates);
569 
570  return toBecomeZero;
571 }
int i
Definition: BiCGSTAB_step_by_step.cpp:9
SmallVector< 3 > computeShapeGradientLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the gradient of the shape-function at the given (lab-fixed) position.
Definition: SuperQuadricParticle.cc:472
Mdouble computeShape(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the shape-functiion at the given (lab-fixed) position.
Definition: SuperQuadricParticle.cc:452

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

Referenced by computeContactPoint().

◆ computeShape()

Mdouble SuperQuadricParticle::computeShape ( const LabFixedCoordinates labFixedCoordinates) const

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

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

Returns
The value of the shape-function at the given (lab-fixed) coordinates.
453 {
454  BodyFixedCoordinates bodyFixedCoordinates = labFixedCoordinates - this->getPosition();
455  getOrientation().rotateBack(bodyFixedCoordinates);
456 
457  const Mdouble n1 = 2. / eps1_;
458  const Mdouble n2 = 2. / eps2_;
459 
460  return std::pow(std::pow(std::abs(bodyFixedCoordinates.X / axes_.X), n2)
461  + std::pow(std::abs(bodyFixedCoordinates.Y / axes_.Y), n2), n1 / n2)
462  + std::pow(std::abs(bodyFixedCoordinates.Z / axes_.Z), n1) - 1.0;
463 }

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

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

◆ computeShapeGradientLabFixed()

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

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

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

Returns
The gradient of the shape function at the given (lab-fixed) coordinates. Note, that this is the gradient to the lab-fixed coordinates, \f nabla_X F(X)
Todo:
Come up with good expression for when x = y = 0 and n1 < n2
473 {
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  const Mdouble absXoa = std::abs(bodyFixedCoordinates.X / axes_.X);
482  const Mdouble absYob = std::abs(bodyFixedCoordinates.Y / axes_.Y);
483  const Mdouble absZoc = std::abs(bodyFixedCoordinates.Z / axes_.Z);
484 
485  const Mdouble nu = std::pow(absXoa, n2) + std::pow(absYob, n2);
486 
487  const Mdouble help1 = std::pow(nu, n1 / n2 - 1.0);
488 
489  Vec3D gradientBodyFixed = {
490  (n1 / axes_.X) * std::pow(absXoa, n2 - 1.0) * help1 * mathsFunc::sign(bodyFixedCoordinates.X),
491  (n1 / axes_.Y) * std::pow(absYob, n2 - 1.0) * help1 * mathsFunc::sign(bodyFixedCoordinates.Y),
492  (n1 / axes_.Z) * std::pow(absZoc, n1 - 1.0) * mathsFunc::sign(bodyFixedCoordinates.Z)};
493  //we get the gradient in terms of the lab-fixed coordinates by \nabla_X F(X) = A \nabla_x f(x).
494  getOrientation().rotate(gradientBodyFixed);
495  return {gradientBodyFixed.X, gradientBodyFixed.Y, gradientBodyFixed.Z};
496 }
void rotate(Vec3D &position) const
Definition: Quaternion.cc:550

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

Referenced by computeResidualContactDetection(), getCurvature(), getInteractionWithSuperQuad(), getJacobianOfContactDetectionObjective(), and overlapFromContactPoint().

◆ copy()

SuperQuadricParticle * SuperQuadricParticle::copy ( ) const
overridevirtual

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

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

Returns
pointer to the particle's copy

Implements NonSphericalParticle.

70 {
71  return new SuperQuadricParticle(*this);
72 }
SuperQuadricParticle()
Basic Particle constructor, creates a superquadric with axes (1,1,1) and exponents (2,...
Definition: SuperQuadricParticle.cc:17

References SuperQuadricParticle().

Referenced by ContactDetectionTester::setupParticles().

◆ getAxes()

Vec3D SuperQuadricParticle::getAxes ( ) const
overridevirtual

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

Todo:

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

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

Reimplemented from BaseParticle.

162 {
163  return axes_;
164 }

References axes_.

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

◆ getContactPoint()

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

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

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

Parameters
[in]pthe particle with which we want to know the contact point with.
[in]CThe contact between this and the give superquadric. Does not contain information of this time-step yet. This contact is used to see if there was an interaction between p and this superquadric in the time-step before, and if so, that contact point can serve as an initial point for contact-detection between the particles.
Returns
A SmallVector of length 4, with the contact point and the lagrange-multiplier needed for Newton's method for finding the contact point.
351 {
352  //Assuming contact between two spheres
353  SmallVector<4> approximateContactPoint = getInitialGuessForContact(p, C);
354  if (computeContactPoint(approximateContactPoint, this, p))
355  {
356  return approximateContactPoint;
357  }
358  else
359  {
360  return getContactPointPlanB(p, 4);
361  }
362 }
SmallVector< 4 > getContactPointPlanB(const SuperQuadricParticle *pOther, unsigned numberOfSteps) const
If the "normal" procedure fails to find a contact point, use an alternative approach that involves st...
Definition: SuperQuadricParticle.cc:699
bool computeContactPoint(SmallVector< 4 > &contactPoint, const SuperQuadricParticle *p1, const SuperQuadricParticle *p2) const
Perform the actual Newton-iterations to find the contact point. Note, that it is given back as a para...
Definition: SuperQuadricParticle.cc:795
SmallVector< 4 > getInitialGuessForContact(const SuperQuadricParticle *pQuad, BaseInteraction *C) const
Get an initial guess for the contact-point between this particle and the given particle.
Definition: SuperQuadricParticle.cc:414
Definition: matrices.h:74

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

Referenced by getInteractionWithSuperQuad(), and isInContactWith().

◆ getContactPointPlanB()

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

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

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

Parameters
[in]pOtherParticle we want to know if the contact-point with
Returns
The point where the shape-functions of both particles are minimised.
700 {
701  logger(VERBOSE, "Number of iterations: %", numberOfSteps);
702  // Step 1: Compute contact point for two volume equivalent spheres
703  const Mdouble interactionRadiusThis = getInteractionRadius(pOther);
704  const Mdouble interactionRadiusOther = pOther->getInteractionRadius(this);
705  const Vec3D positionThis = getPosition();
706  const Vec3D positionOther = pOther->getPosition();
707 
708  //analytical solution for contact point between bounding spheres
709  const Vec3D initialPoint = (interactionRadiusOther * positionThis + interactionRadiusThis * positionOther) /
710  (interactionRadiusThis + interactionRadiusOther);
711 
712  SmallVector<4> contactPointPlanB;
713  // Analytical solution
714  contactPointPlanB[0] = initialPoint.X;
715  contactPointPlanB[1] = initialPoint.Y;
716  contactPointPlanB[2] = initialPoint.Z;
717  contactPointPlanB[3] = 1.0; //Need to check.
718 
719  writeDebugMessageStep1(pOther, contactPointPlanB);
720 
721  //Step 2: Compute the deltas
722  const Vec3D dAxesThis = (getAxes() - interactionRadiusThis * Vec3D(1, 1, 1)) / numberOfSteps;
723  const Mdouble dn11 = (2.0 / getExponentEps1() - 2.0) / numberOfSteps;
724  const Mdouble dn12 = (2.0 / getExponentEps2() - 2.0) / numberOfSteps;
725 
726  const Vec3D dAxesOther = (pOther->getAxes() - interactionRadiusOther * Vec3D(1, 1, 1)) / numberOfSteps;
727  const Mdouble dn21 = (2.0 / pOther->getExponentEps1() - 2.0) / numberOfSteps;
728  const Mdouble dn22 = (2.0 / pOther->getExponentEps2() - 2.0) / numberOfSteps;
729 
730  writeDebugMessageStep2(pOther, dAxesThis, dn11, dn12, dAxesOther, dn21, dn22);
731 
732  // Create two superquadrics with the above parameters for incrementally computing the contact point
735 
736  p1.setPosition(getPosition());
737  p2.setPosition(pOther->getPosition());
738 
739  p1.setOrientation(getOrientation());
740  p2.setOrientation(pOther->getOrientation());
741  bool success = true;
742  const unsigned maxIterations = 20;
743  for (unsigned int counter = 1; counter <= numberOfSteps; ++counter)
744  {
745 
746  // Step 3: Increment the shape and blockiness parameters
747  const Vec3D axesThis = interactionRadiusThis * Vec3D(1, 1, 1) + counter * dAxesThis;
748  const Mdouble n11 = 2.0 + counter * dn11;
749  const Mdouble n12 = 2.0 + counter * dn12;
750 
751  Vec3D axesOther = interactionRadiusOther * Vec3D(1, 1, 1) + counter * dAxesOther;
752  const Mdouble n21 = 2.0 + counter * dn21;
753  const Mdouble n22 = 2.0 + counter * dn22;
754 
755  p1.setAxesAndExponents(axesThis, 2.0 / n11, 2.0 / n12);
756  p2.setAxesAndExponents(axesOther, 2.0 / n21, 2.0 / n22);
757 
758  writeDebugMessageStep3(axesThis, n11, n12, axesOther, n21, n22);
759 
760  writeDebugMessageMiddleOfLoop(p1, p2, contactPointPlanB, counter);
761 
762  //compute the contact point of the new particles
763  success = computeContactPoint(contactPointPlanB, &p1, &p2);
764  if (!success)
765  {
766  if (numberOfSteps > maxIterations)
767  {
768  write(std::cout);
769  pOther->write(std::cout);
770  logger(ERROR, "Plan B fails even with more than 20 intermediate steps");
771  }
772  return getContactPointPlanB(pOther, 2 * numberOfSteps);
773  }
774  }
775  logger.assert_debug(p1.getAxes().X == getAxes().X, "In getContactPointPlanB, final particle for contact-detection not "
776  "the same as original particle");
777 
778  logger(VERBOSE, "Plan B contact point: %", contactPointPlanB);
779 
780  return contactPointPlanB;
781 }
virtual void setOrientation(const Quaternion &orientation)
Sets the orientation of this BaseInteractable.
Definition: BaseInteractable.h:239
virtual void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Definition: BaseInteractable.h:218
Definition: SuperQuadricParticle.h:36
void writeDebugMessageStep1(const SuperQuadricParticle *pQuad, const SmallVector< 4 > &contactPointPlanB) const
Definition: SuperQuadricParticle.cc:932
Vec3D getAxes() const override
Get the axes-lengths of this superquadric. We use the super-ellipsoid definition stated in Chapter 2 ...
Definition: SuperQuadricParticle.cc:161
Mdouble getExponentEps2() const override
Get the second exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter...
Definition: SuperQuadricParticle.cc:171
void writeDebugMessageStep3(const Vec3D &axesThis, const Mdouble &n11, const Mdouble &n12, const Vec3D &axesOther, const Mdouble &n21, const Mdouble &n22) const
Definition: SuperQuadricParticle.cc:870
void writeDebugMessageMiddleOfLoop(const SuperQuadricParticle &p1, const SuperQuadricParticle &p2, SmallVector< 4 > &contactPointPlanB, const unsigned int &counter) const
Definition: SuperQuadricParticle.cc:848
void writeDebugMessageStep2(const SuperQuadricParticle *pQuad, const Vec3D &dAxesThis, const Mdouble &dn11, const Mdouble &dn12, const Vec3D &dAxesOther, const Mdouble &dn21, const Mdouble &dn22) const
Definition: SuperQuadricParticle.cc:894
Mdouble getInteractionRadius(const BaseParticle *particle) const
returns the radius plus half the interactionDistance of the mixed species
Definition: SuperQuadricParticle.cc:955
void write(std::ostream &os) const override
Write function: write this SuperQuadric to the given output-stream, for example a restart-file.
Definition: SuperQuadricParticle.cc:79
Mdouble getExponentEps1() const override
Get the first exponent of this superquadric. We use the super-ellipsoid definition stated in Chapter ...
Definition: SuperQuadricParticle.cc:166
void setAxesAndExponents(const Mdouble &a1, const Mdouble &a2, const Mdouble &a3, const Mdouble &eps1, const Mdouble &eps2)
Set the geometrical properties of the superquadrics, namely the axes-lengths a1, a2 and a3,...
Definition: SuperQuadricParticle.cc:112
#define X
Definition: icosphere.cpp:20

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

Referenced by getContactPoint().

◆ getCurvature()

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

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

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

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

Reimplemented from BaseParticle.

644 {
645  SmallVector<3> gradientVec = computeShapeGradientLabFixed(labFixedCoordinates);
646  SmallMatrix<3, 1> gradient = gradientVec;
647  SmallMatrix<3, 3> hessian = computeHessianLabFixed(labFixedCoordinates);
648  Mdouble helper = ((gradient.transpose() * hessian) * gradient)(0, 0);
649  Mdouble gradientLength = gradientVec.length();
650  Mdouble helper2 = gradientLength * gradientLength * (hessian(0, 0) + hessian(1, 1) + hessian(2, 2));
651  Mdouble denominator = 2 * gradientLength * gradientLength * gradientLength;
652 
654  return (helper2 - helper) / denominator;
655 }
SmallMatrix< numberOfColumns, numberOfRows > transpose() const
Definition: SmallMatrix.h:310
Mdouble length() const
Definition: SmallVector.h:233
SmallMatrix< 3, 3 > computeHessianLabFixed(const LabFixedCoordinates &labFixedCoordinates) const
Compute and get the hessian ("second derivative") of the shape-function at the given (lab-fixed) posi...
Definition: SuperQuadricParticle.cc:505

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

◆ getExponentEps1()

Mdouble SuperQuadricParticle::getExponentEps1 ( ) const
overridevirtual

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

Reimplemented from BaseParticle.

167 {
168  return eps1_;
169 }

References eps1_.

Referenced by getContactPointPlanB(), and writeDebugMessageStep2().

◆ getExponentEps2()

Mdouble SuperQuadricParticle::getExponentEps2 ( ) const
overridevirtual

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

Reimplemented from BaseParticle.

172 {
173  return eps2_;
174 }

References eps2_.

Referenced by getContactPointPlanB(), and writeDebugMessageStep2().

◆ getInitialGuessForContact()

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

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

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

Parameters
pQuadother superquadric for which we want to know if there is a contact
Ccontact between this particle and the given particle
Returns
A guess for the initial contact: contact point X,Y,Z, lagrange multiplier. Note that these are given in the LAB-FIXED coordinates.
415 {
416  SmallVector<4> approximateContactPoint;
417  if (C == nullptr || C->getOverlap() < 1e-15)
418  {
419  //this is a new interaction
420  const LabFixedCoordinates branchVector = pQuad->getPosition() - getPosition();
421  //Get the square of the distance between particle i and particle j
422  const Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
423  const Mdouble distance = sqrt(distanceSquared);
424  const LabFixedCoordinates normal = (branchVector / distance);
425  const Mdouble overlap = getSumOfInteractionRadii(pQuad) - distance;
426  const LabFixedCoordinates contactPoint = (pQuad->getPosition() -
427  (pQuad->getInteractionRadius(this) - 0.5 * overlap) * normal);
428 
429  approximateContactPoint[0] = contactPoint.X;
430  approximateContactPoint[1] = contactPoint.Y;
431  approximateContactPoint[2] = contactPoint.Z;
432  approximateContactPoint[3] = 1;
433  }
434  else
435  {
436  //this is an existing interaction
437  const LabFixedCoordinates contactPoint = C->getContactPoint();
438  const Mdouble multiplier = C->getLagrangeMultiplier();
439  approximateContactPoint[0] = contactPoint.X;
440  approximateContactPoint[1] = contactPoint.Y;
441  approximateContactPoint[2] = contactPoint.Z;
442  approximateContactPoint[3] = multiplier;
443  }
444  return approximateContactPoint;
445 }
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
Mdouble getSumOfInteractionRadii(const BaseParticle *particle) const
returns the sum of the radii plus the interactionDistance
Definition: BaseParticle.h:362
Mdouble getLengthSquared() const
Calculates the squared length of this Vec3D: .
Definition: Vector.cc:164
double multiplier(const Vector< double > &xi)
Definition: disk_oscillation.cc:63
void normal(const Vector< double > &x, Vector< double > &normal)
Definition: free_surface_rotation.cc:65

References e(), getInteractionRadius(), Vec3D::getLengthSquared(), BaseInteractable::getPosition(), BaseParticle::getSumOfInteractionRadii(), Global_Physical_Variables::multiplier(), WallFunction::normal(), sqrt(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getContactPoint().

◆ getInteractionRadius()

Mdouble SuperQuadricParticle::getInteractionRadius ( const BaseParticle particle) const

returns the radius plus half the interactionDistance of the mixed species

956 {
957  const auto mixedSpecies = getSpecies()->getHandler()->getMixedObject(getSpecies(),particle->getSpecies());
958  return getRadius() + 0.5 * mixedSpecies->getInteractionDistance();
959 }
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
Definition: BaseInteractable.h:87
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:331
SpeciesHandler * getHandler() const
Returns the pointer to the handler to which this species belongs.
Definition: BaseSpecies.cc:78
std::enable_if<!std::is_pointer< typename U::MixedSpeciesType >::value, typename U::MixedSpeciesType * >::type getMixedObject(const U *S, const U *T)
Definition: SpeciesHandler.h:52

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

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

◆ getInteractionWith()

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

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

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

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

Reimplemented from BaseParticle.

267 {
268  //get the normal (from P away from the contact)
269  const LabFixedCoordinates branchVector = p->getPosition() - getPosition();
270  //Get the square of the distance between particle i and particle j
271  const Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
272 
273  const Mdouble sumOfInteractionRadii = getMaxInteractionRadius()+p->getMaxInteractionRadius();
274  if (distanceSquared < (sumOfInteractionRadii * sumOfInteractionRadii))
275  {
276  if (p->isSphericalParticle())
277  {
279  pQuad->setAxes(p->getRadius(), p->getRadius(), p->getRadius());
280  BaseInteraction* contacts = getInteractionWithSuperQuad(pQuad, timeStamp, interactionHandler);
281  delete pQuad;
282  return contacts;
283  } else {
284  SuperQuadricParticle* pQuad = static_cast<SuperQuadricParticle*>(p);
285  return getInteractionWithSuperQuad(pQuad, timeStamp, interactionHandler);
286  }
287  }
288  return nullptr;
289 }
Stores information about interactions between two interactable objects; often particles but could be ...
Definition: BaseInteraction.h:39
Mdouble getMaxInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e....
Definition: BaseParticle.h:345
BaseInteraction * getInteractionWithSuperQuad(SuperQuadricParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler)
Checks if this superquadric is in interaction with the given superquadric, and if so,...
Definition: SuperQuadricParticle.cc:301
void setAxes(const Mdouble &a1, const Mdouble &a2, const Mdouble &a3)
Set the axes-lengths to a1, a2 and a3 for this superquadric. We use the super-ellipsoid definition st...
Definition: SuperQuadricParticle.cc:143

References getInteractionWithSuperQuad(), Vec3D::getLengthSquared(), BaseParticle::getMaxInteractionRadius(), BaseInteractable::getPosition(), p, setAxes(), and SuperQuadricParticle().

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

◆ getInteractionWithSuperQuad()

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

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

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

Parameters
[in]pSuperquadric particle we want to have the contact-information for
[in]timeStampTime stamp to be assigned to the interaction object (i.e., the current time)
[in]interactionHandlerBaseInteraction container from where the interaction is retrieved, and to which it is assigned (if it is a new interaction).
Returns
A vector of Interaction* with size 1 (if there is an interaction) or 0 (if there is no interaction).
Todo:
: find correct value for 'distance'
303 {
304  BaseInteraction* const C = interactionHandler->getInteraction(p, this, timeStamp);
305  const SmallVector<4> approximateContactPoint = getContactPoint(p, C);
306 
307  //set the new contact point:
308  LabFixedCoordinates contactPoint;
309  contactPoint.X = approximateContactPoint[0];
310  contactPoint.Y = approximateContactPoint[1];
311  contactPoint.Z = approximateContactPoint[2];
312 
313  //if the minimum is positive, there is no contact: return empty vector
314  if (computeShape(contactPoint) > -1e-10)
315  {
316  interactionHandler->removeObject(C->getIndex());
317  return nullptr;
318  }
319 
320  C->setContactPoint(contactPoint);
321  C->setLagrangeMultiplier(approximateContactPoint[3]);
322 
323  SmallVector<3> gradThis = computeShapeGradientLabFixed(contactPoint).getNormalised();
325  normal.X = gradThis[0];
326  normal.Y = gradThis[1];
327  normal.Z = gradThis[2];
328  C->setNormal(normal);
329 
330  const Mdouble alphaI = overlapFromContactPoint(contactPoint, normal);
331  const Mdouble alphaJ = p->overlapFromContactPoint(contactPoint, -normal);
332  C->setOverlap(alphaI + alphaJ);
334  C->setDistance((getPosition() - p->getPosition()).getLength() - C->getOverlap());
335 
336  return C;
337 }
Matrix< Scalar, Dynamic, Dynamic > C
Definition: bench_gemm.cpp:49
virtual void removeObject(unsigned const int index)
Removes an Object from the BaseHandler.
Definition: BaseHandler.h:453
BaseInteraction * getInteraction(BaseInteractable *P, BaseInteractable *I, unsigned timeStamp)
Returns the Interaction between the BaseInteractable's P and I.
Definition: InteractionHandler.cc:126
SmallVector< 4 > getContactPoint(const SuperQuadricParticle *p, BaseInteraction *C) const
Compute the contact point between this and the given superquadric particle.
Definition: SuperQuadricParticle.cc:350
Mdouble overlapFromContactPoint(const LabFixedCoordinates &contactPoint, const LabFixedCoordinates &normal) const
Compute the distance between the contact-point and surface of this superquadric particle.
Definition: SuperQuadricParticle.cc:371

References computeShape(), computeShapeGradientLabFixed(), e(), getContactPoint(), InteractionHandler::getInteraction(), BaseInteractable::getPosition(), WallFunction::normal(), overlapFromContactPoint(), p, BaseHandler< T >::removeObject(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getInteractionWith().

◆ getJacobianOfContactDetectionObjective()

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

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

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

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

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

Parameters
[in]contactPointThe 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.
605 {
606  LabFixedCoordinates labFixedCoordinates;
607  labFixedCoordinates.X = contactPoint[0];
608  labFixedCoordinates.Y = contactPoint[1];
609  labFixedCoordinates.Z = contactPoint[2];
610  const Mdouble lagrangeMultiplier = contactPoint[3];
611 
612  SmallVector<3> gradThis = p1->computeShapeGradientLabFixed(labFixedCoordinates);
613  SmallVector<3> gradOther = p2->computeShapeGradientLabFixed(labFixedCoordinates);
614 
615  SmallMatrix<3, 3> hessianThis = p1->computeHessianLabFixed(labFixedCoordinates);
616  SmallMatrix<3, 3> hessianOther = p2->computeHessianLabFixed(labFixedCoordinates);
617 
618  SmallMatrix<3, 3> upperLeftCorner = hessianThis + lagrangeMultiplier * lagrangeMultiplier * hessianOther;
619 
620  SmallMatrix<4, 4> jacobian;
621  for (unsigned int i = 0; i < 3; ++i)
622  {
623  for (unsigned int j = 0; j < 3; ++j)
624  {
625  jacobian(i, j) = upperLeftCorner(i, j);
626  }
627  jacobian(i, 3) = 2 * lagrangeMultiplier * gradOther[i];
628  }
629 
630  for (unsigned int j = 0; j < 3; ++j)
631  {
632  jacobian(3, j) = gradThis[j] - gradOther[j];
633  }
634 
635  return jacobian;
636 }
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

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

Referenced by computeContactPoint().

◆ getName()

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

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

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

Returns
The std::string "SuperQuadricParticle".

Implements NonSphericalParticle.

92 {
93  return "SuperQuadricParticle";
94 }

◆ getVolume()

Mdouble SuperQuadricParticle::getVolume ( ) const
overridevirtual

Get the volume of this superquadric.

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

Returns
The actual volume of this superquadric.

Reimplemented from BaseParticle.

184 {
185  logger.assert_debug(getHandler() != nullptr, "[SuperQuadricParticle::getVolume] no particle handler specified");
186 
187  return (2.0 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) * mathsFunc::beta(0.5 * eps1_ + 1.0, eps1_) *
188  mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_);
189 }

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

Referenced by computeMass().

◆ isInContactWith()

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

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

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

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

Reimplemented from BaseParticle.

665 {
666  SmallVector<4> approximateContactPoint;
667 
668  if (p->isSphericalParticle())
669  {
671  pQuad->setAxes(p->getRadius(), p->getRadius(), p->getRadius());
672  approximateContactPoint = getContactPoint(pQuad, nullptr);
673  delete pQuad;
674  } else {
675  const SuperQuadricParticle* pQuad = static_cast<const SuperQuadricParticle*>(p);
676  approximateContactPoint = getContactPoint(pQuad, nullptr);
677  }
678 
679  //set the new contact point:
680  LabFixedCoordinates contactPoint;
681  contactPoint.X = approximateContactPoint[0];
682  contactPoint.Y = approximateContactPoint[1];
683  contactPoint.Z = approximateContactPoint[2];
684 
685  //if the minimum is positive, there is no contact: return false
686  return (computeShape(contactPoint) < 0);
687 
688 }

References computeShape(), getContactPoint(), p, setAxes(), SuperQuadricParticle(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by ContactDetectionTester::testEllipsoidsContact().

◆ overlapFromContactPoint()

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

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

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

Parameters
[in]contactPointThe contact point between this particle with another BaseInteractable
[in]normalThe normal direction of the contact.
Returns
Distance between contact-point and surface of this particle
373 {
374  Mdouble alphaI = 0;
375  LabFixedCoordinates xEdge = contactPoint + alphaI * normal;
376  Mdouble dampingCoefficient = 1;
377  while (std::abs(computeShape(xEdge)) > 1e-10)
378  {
379  SmallVector<3> gradientShape = computeShapeGradientLabFixed(xEdge);
380  LabFixedCoordinates gradient;
381  gradient.X = gradientShape(0);
382  gradient.Y = gradientShape(1);
383  gradient.Z = gradientShape(2);
384  Mdouble newValue = alphaI - dampingCoefficient * computeShape(xEdge) / Vec3D::dot(gradient, normal);
385  const LabFixedCoordinates xEdgeNew = contactPoint + newValue * normal;
386 
387  if (std::abs(computeShape(xEdgeNew)) < std::abs(computeShape(xEdge)))
388  {
389  alphaI = newValue;
390  xEdge = xEdgeNew;
391  dampingCoefficient = 1;
392  }
393  else
394  {
395  dampingCoefficient /= 2;
396  if (dampingCoefficient < 1e-10)
397  {
398  logger(ERROR, "overlap detection does not converge");
399  }
400  }
401  }
402  return alphaI;
403 }
static Mdouble dot(const Vec3D &a, const Vec3D &b)
Calculates the dot product of two Vec3D: .
Definition: Vector.cc:56

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

Referenced by getInteractionWithSuperQuad().

◆ read()

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

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

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

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

Reimplemented from BaseParticle.

102 {
103  BaseParticle::read(is);
104  std::string dummy;
105  is >> dummy >> axes_ >> dummy >> eps1_ >> dummy >> eps2_;
106  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
107  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
108 }
void read(std::istream &is) override
Particle read function, which accepts an std::istream as input.
Definition: BaseParticle.cc:368
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286

References axes_, e(), eps1_, eps2_, logger, BaseParticle::read(), and oomph::Global_string_for_annotation::string().

◆ setAxes() [1/2]

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

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

144 {
145  setAxes({a1,a2,a3});
146 }

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

◆ setAxes() [2/2]

void SuperQuadricParticle::setAxes ( const Vec3D axes)
overridevirtual

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

Reimplemented from BaseParticle.

134 {
135  axes_ = axes;
137  if (getSpecies() != nullptr)
138  {
139  getSpecies()->computeMass(this);
140  }
141 }
void computeMass(BaseParticle *p) const
Compute Particle mass function, which required a reference to the Species vector. It computes the Par...
Definition: ParticleSpecies.cc:147
void setBoundingRadius()
Get the radius of the sphere that fits precisely around the particle.
Definition: SuperQuadricParticle.cc:223

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

◆ setAxesAndExponents() [1/2]

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

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.

114 {
115  eps1_ = eps1;
116  eps2_ = eps2;
117  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
118  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
119 
120  setAxes(a1,a2,a3);
121 }

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

Referenced by getContactPointPlanB(), GranularCollapse::setupInitialConditions(), ContactDetectionWithWallTester::setupParticleAndWall(), and ContactDetectionTester::setupParticles().

◆ setAxesAndExponents() [2/2]

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

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

124 {
125  eps1_ = eps1;
126  eps2_ = 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 
130  setAxes(axes);
131 }

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

◆ setBoundingRadius()

void SuperQuadricParticle::setBoundingRadius ( )
private

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

Todo:
Currently only implemented for ellipsoids
224 {
225 
227  {
229  return;
230  }else
231  {
232 
233  const Mdouble axesX = std::max(axes_.X, axes_.Y);
234  const Mdouble axesY = std::min(axes_.X, axes_.Y);
235 
236  Mdouble alpha;
237 
238  const Mdouble eps1 = std::min(.96, eps1_);
239  const Mdouble eps2 = std::min(.96, eps2_);
240 
241  alpha = std::pow(axesY / axesX, 2.0 / (2.0 / eps2 - 2.0));
242 
243  const Mdouble help1 = std::pow(alpha, 2.0 / eps2);
244  const Mdouble gamma = std::pow(1.0 + help1, eps2 / eps1 - 1.0);
245  const Mdouble beta = std::pow(gamma * axes_.Z * axes_.Z / (axesX * axesX), 1.0 / (2.0 / eps1 - 2.0));
246  const Mdouble xTilde = std::pow(std::pow(1 + help1, eps2 / eps1) + std::pow(beta, 2.0 / eps1),
247  -eps1 / 2.0);
249  + mathsFunc::square(alpha * axesY * xTilde)
250  + mathsFunc::square(beta * axes_.Z * xTilde)));
251 
252  }
253 }
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species)
Definition: BaseParticle.cc:548
#define min(a, b)
Definition: datatypes.h:22
#define max(a, b)
Definition: datatypes.h:23
RealScalar alpha
Definition: level1_cplx_impl.h:151
Scalar beta
Definition: level2_cplx_impl.h:36
T square(const T val)
squares a number
Definition: ExtendedMath.h:86
bool isEqual(Mdouble v1, Mdouble v2, Mdouble absError)
Compares the difference of two Mdouble with an absolute error, useful in UnitTests.
Definition: ExtendedMath.cc:230
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
Definition: ExtendedMath.cc:116
double epsilon
Definition: osc_ring_sarah_asymptotics.h:43

References alpha, axes_, beta, eps1_, eps2_, oomph::SarahBL::epsilon, mathsFunc::gamma(), mathsFunc::isEqual(), max, min, Eigen::bfloat16_impl::pow(), BaseParticle::setRadius(), sqrt(), mathsFunc::square(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by setAxes(), and setExponents().

◆ setExponents()

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

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

Reimplemented from BaseParticle.

149 {
150  eps1_ = eps1;
151  eps2_ = eps2;
153  logger.assert_always(eps1_ < 1 + 1e-10, "epsilon1 should be at most 1");
154  logger.assert_always(eps2_ < 1 + 1e-10, "epsilon2 should be at most 1");
155  if (getSpecies() != nullptr)
156  {
157  getSpecies()->computeMass(this);
158  }
159 }

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

◆ setInertia()

void SuperQuadricParticle::setInertia ( )
overridevirtual

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

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

Reimplemented from BaseParticle.

196 {
197  MatrixSymmetric3D inertia;
198 
199  inertia.XX = getSpecies()->getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
200  (axes_.Y * axes_.Y * mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_) *
201  mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0)
202  + 4.0 * axes_.Z * axes_.Z * mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_ + 1) *
203  mathsFunc::beta(1.5 * eps1_, eps1_ + 1.0));
204 
205  inertia.YY = getSpecies()->getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
206  (axes_.X * axes_.X * mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_) *
207  mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0)
208  + 4.0 * axes_.Z * axes_.Z * mathsFunc::beta(0.5 * eps2_, 0.5 * eps2_ + 1) *
209  mathsFunc::beta(1.5 * eps1_, eps1_ + 1.0));
210 
211  inertia.ZZ = getSpecies()->getDensity() * (0.5 * axes_.X * axes_.Y * axes_.Z * eps1_ * eps2_) *
212  (axes_.X * axes_.X + axes_.Y * axes_.Y) * mathsFunc::beta(1.5 * eps2_, 0.5 * eps2_) *
213  mathsFunc::beta(0.5 * eps1_, 2.0 * eps1_ + 1.0);
214 
215  BaseParticle::setInertia(inertia);
216 }
virtual void setInertia()
Definition: BaseParticle.cc:500
Implementation of a 3D symmetric matrix.
Definition: MatrixSymmetric.h:16
Mdouble getDensity() const
Allows density_ to be accessed.
Definition: ParticleSpecies.cc:98

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

Referenced by EllipsoidsBouncingOnWallDemo::setupInitialConditions(), GranularCollapse::setupInitialConditions(), ContactDetectionWithWallTester::setupParticleAndWall(), and ContactDetectionTester::setupParticles().

◆ setRadius()

void SuperQuadricParticle::setRadius ( const Mdouble  radius)
overridevirtual

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

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

Parameters
[in]radiusthe new radius

Reimplemented from BaseParticle.

219 {
220  logger(ERROR,"This function should not be used");
221 }

References ERROR, and logger.

◆ write()

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

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

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

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

Reimplemented from BaseParticle.

80 {
82  os << " axes " << axes_
83  << " exp1 " << eps1_
84  << " exp2 " << eps2_;
85 }
void write(std::ostream &os) const override
Particle print function, which accepts an std::ostream as input.
Definition: BaseParticle.cc:319

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

Referenced by getContactPointPlanB().

◆ writeDebugMessageMiddleOfLoop()

void SuperQuadricParticle::writeDebugMessageMiddleOfLoop ( const SuperQuadricParticle p1,
const SuperQuadricParticle p2,
SmallVector< 4 > &  contactPointPlanB,
const unsigned int counter 
) const
850 {
851  logger(VERBOSE, "Position of particle 1 (p1): % \nPosition of particle 2 (p2): % \n",
852  p1.getPosition(), p2.getPosition());
853  logger(VERBOSE, "Orientation of particle 1 (p1): % \nOrientation of particle 2 (p2): %\n",
854  p1.getOrientation(), p2.getOrientation());
855 
856  // Step 4: Calculate the contact point for the updated shape parameters
857 
858  logger(VERBOSE, "----------------------------------------");
859  logger(VERBOSE, "STEP 4: Compute contact point");
860  logger(VERBOSE, "----------------------------------------");
861  logger(VERBOSE, " ");
862  logger(VERBOSE, "Counter: %", counter);
863  logger(VERBOSE, " ");
864 
865  logger(VERBOSE, "Analytical solution Contact Point: % ", contactPointPlanB);
866  logger(VERBOSE, " ");
867 }

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

Referenced by getContactPointPlanB().

◆ writeDebugMessageStep1()

void SuperQuadricParticle::writeDebugMessageStep1 ( const SuperQuadricParticle pQuad,
const SmallVector< 4 > &  contactPointPlanB 
) const
933 {
934  logger(VERBOSE, "---------------------");
935  logger(VERBOSE, "STEP 1");
936  logger(VERBOSE, "---------------------\n");
937 
938  logger(VERBOSE, "Position of particle 1: % ", getPosition());
939  logger(VERBOSE, "Position of particle 2 (pQuad): % ", pQuad->getPosition());
940  logger(VERBOSE, " ");
941 
942  logger(VERBOSE, "Radius of particle 1: % ", getInteractionRadius(pQuad));
943  logger(VERBOSE, "Radius of particle 2 (pQuad): % \n", pQuad->getInteractionRadius(this));
944 
945  logger(VERBOSE, "Orientation of particle 1: % ", getOrientation());
946  logger(VERBOSE, "Orientation of particle 2 (pQuad): % ", pQuad->getOrientation());
947  logger(VERBOSE, " ");
948 
949  logger(VERBOSE, "Particle 1 axes: %", getAxes());
950  logger(VERBOSE, "Particle 2 axes (pQuad): % \n", pQuad->getAxes());
951 
952  logger(VERBOSE, "Analytical solution for two equivalent spheres in contact: % \n", contactPointPlanB);
953 }

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

Referenced by getContactPointPlanB().

◆ writeDebugMessageStep2()

void SuperQuadricParticle::writeDebugMessageStep2 ( const SuperQuadricParticle pQuad,
const Vec3D dAxesThis,
const Mdouble dn11,
const Mdouble dn12,
const Vec3D dAxesOther,
const Mdouble dn21,
const Mdouble dn22 
) const
897 {
898  logger(VERBOSE, "---------------------");
899  logger(VERBOSE, "STEP 2");
900  logger(VERBOSE, "---------------------");
901  logger(VERBOSE, " ");
902 
903  logger(VERBOSE, "Particle 1 x-scale: %", pQuad->getAxes().X);
904  logger(VERBOSE, "Particle 1 y-scale: %", pQuad->getAxes().Y);
905  logger(VERBOSE, "Particle 1 z-scale: %", pQuad->getAxes().Z);
906  logger(VERBOSE, "Particle 1 n1: %", 2. / pQuad->getExponentEps1());
907  logger(VERBOSE, "Particle 1 n2: %", 2. / pQuad->getExponentEps2());
908  logger(VERBOSE, " ");
909 
910  logger(VERBOSE, "Particle 1 x-scale increment: %", dAxesThis.X);
911  logger(VERBOSE, "Particle 1 y-scale increment: %", dAxesThis.Y);
912  logger(VERBOSE, "Particle 1 z-scale increment: %", dAxesThis.Z);
913  logger(VERBOSE, "Particle 1 n1 increment: %", dn11);
914  logger(VERBOSE, "Particle 1 n2 increment: %", dn12);
915  logger(VERBOSE, " ");
916 
917  logger(VERBOSE, "Particle 2 x-scale: %", getAxes().X);
918  logger(VERBOSE, "Particle 2 y-scale: %", getAxes().Y);
919  logger(VERBOSE, "Particle 2 z-scale: %", getAxes().Z);
920  logger(VERBOSE, "Particle 2 n1: %", 2. / getExponentEps1());
921  logger(VERBOSE, "Particle 2 n2: %", 2. / getExponentEps2());
922  logger(VERBOSE, " ");
923 
924  logger(VERBOSE, "Particle 2 x-scale increment: %", dAxesOther.X);
925  logger(VERBOSE, "Particle 2 y-scale increment: %", dAxesOther.Y);
926  logger(VERBOSE, "Particle 2 z-scale increment: %", dAxesOther.Z);
927  logger(VERBOSE, "Particle 2 n1 increment: %", dn21);
928  logger(VERBOSE, "Particle 2 n2 increment: %", dn22);
929  logger(VERBOSE, " ");
930 }
#define Z
Definition: icosphere.cpp:21
const char Y
Definition: test/EulerAngles.cpp:32

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

Referenced by getContactPointPlanB().

◆ writeDebugMessageStep3()

void SuperQuadricParticle::writeDebugMessageStep3 ( const Vec3D axesThis,
const Mdouble n11,
const Mdouble n12,
const Vec3D axesOther,
const Mdouble n21,
const Mdouble n22 
) const
872 {
873  logger(VERBOSE, "-----------------------------------");
874  logger(VERBOSE, "STEP 3: First increment");
875  logger(VERBOSE, "-----------------------------------");
876  logger(VERBOSE, " ");
877 
878  logger(VERBOSE, "Particle 1 x-scale after increment: %", axesThis.X);
879  logger(VERBOSE, "Particle 1 y-scale after increment: %", axesThis.Y);
880  logger(VERBOSE, "Particle 1 z-scale after increment: %", axesThis.Z);
881  logger(VERBOSE, "Particle 1 n1 after increment: %", n11);
882  logger(VERBOSE, "Particle 1 n2 after increment: %", n12);
883  logger(VERBOSE, " ");
884 
885  logger(VERBOSE, "Particle 2 x-scale after increment: %", axesOther.X);
886  logger(VERBOSE, "Particle 2 y-scale after increment: %", axesOther.Y);
887  logger(VERBOSE, "Particle 2 z-scale after increment: %", axesOther.Z);
888  logger(VERBOSE, "Particle 2 n1 after increment: %", n21);
889  logger(VERBOSE, "Particle 2 n2 after increment: %", n22);
890  logger(VERBOSE, " ");
891 }

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

Referenced by getContactPointPlanB().

Member Data Documentation

◆ axes_

◆ eps1_

Mdouble SuperQuadricParticle::eps1_
private

Blockiness parameters.

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

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

◆ eps2_


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