MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Screw Class Reference

This function defines an Archimedes' screw in the z-direction from a (constant) starting point, a (constant) length L, a (constant) radius r, a (constant) number or revelations N and a (constant) rotation speed (rev/s) More...

#include <Screw.h>

+ Inheritance diagram for Screw:

Public Member Functions

 Screw ()
 Default constructor: make a screw with default parameters. More...
 
 Screw (const Screw &other)
 Copy constructor, copies another Screw. More...
 
 Screw (Vec3D start, Mdouble l, Mdouble r, Mdouble n, Mdouble omega, Mdouble thickness, ScrewType screwType=ScrewType::doubleHelix)
 Constructor in which all parameters of the screw are set. More...
 
 ~Screw () override
 Default destructor. More...
 
Screwcopy () const final
 Copy this screw and return a pointer to the copy. More...
 
bool getDistanceAndNormal (const BaseParticle &P, Mdouble &distance, Vec3D &normal_return) const final
 Compute the distance from the Screw for a given BaseParticle and return if there is a collision. If there is a collision, also return the normal vector of the interaction point. More...
 
bool getDistanceAndNormalLabCoordinates (Vec3D position, Mdouble wallInteractionRadius, Mdouble &distance, Vec3D &normal_return) const
 
void move_time (Mdouble dt)
 Rotate the Screw for a period dt, so that the offset_ changes with omega_*dt. More...
 
void rotate (const Vec3D &angularVelocityDt) override
 
void read (std::istream &is) override
 Reads a Screw from an input stream, for example a restart file. More...
 
void oldRead (std::istream &is)
 Reads a Screw in the old style from an input stream, for example a restart file old style. More...
 
void write (std::ostream &os) const override
 Writes this Screw to an output stream, for example a restart file. More...
 
std::string getName () const final
 Returns the name of the object, here the string "Screw". More...
 
void writeVTK (VTKContainer &vtk) const override
 
void writeVTK (std::string filename) const
 
- Public Member Functions inherited from BaseWall
 BaseWall ()
 Default constructor. More...
 
 BaseWall (const BaseWall &w)
 Copy constructor. More...
 
 ~BaseWall () override
 Default destructor. More...
 
virtual bool getDistanceNormalOverlap (const BaseParticle &P, Mdouble &distance, Vec3D &normal_return, Mdouble &overlap) const
 
virtual bool getDistanceNormalOverlapSuperquadric (const SuperQuadricParticle &p, Mdouble &distance, Vec3D &normal_return, Mdouble &overlap) const
 
virtual Vec3D getFurthestPointSuperQuadric (const Vec3D &normalBodyFixed, const Vec3D &axes, Mdouble eps1, Mdouble eps2) const
 
virtual void setHandler (WallHandler *handler)
 A function which sets the WallHandler for this BaseWall. More...
 
WallHandlergetHandler () const
 A function which returns the WallHandler that handles this BaseWall. More...
 
void setIndSpecies (unsigned int indSpecies) override
 Define the species of this wall using the index of the species in the SpeciesHandler in this DPMBase. More...
 
void setSpecies (const ParticleSpecies *species)
 Defines the species of the current wall. More...
 
bool isFixed () const override
 
void setForceControl (Vec3D forceGoal, Vec3D gainFactor, Vec3D baseVelocity={0, 0, 0})
 Slowly adjusts the force on a wall towards a specified goal, by adjusting (prescribing) the velocity of the wall. More...
 
virtual bool isLocal (Vec3D &min, Vec3D &max) const
 
bool getLinePlaneIntersect (Vec3D &intersect, const Vec3D &p0, const Vec3D &p1, const Vec3D &n, const Vec3D &p)
 
bool isInsideWallVTK (const Vec3D &point, const Vec3D &normal, const Vec3D &position) const
 
void projectOntoWallVTK (Vec3D &point0, const Vec3D &point1, const Vec3D &normal, const Vec3D &position) const
 
void intersectVTK (std::vector< Vec3D > &points, Vec3D normal, Vec3D position) const
 
virtual BaseInteractiongetInteractionWithSuperQuad (SuperQuadricParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler)
 
void getVTK (std::vector< Vec3D > &points, std::vector< std::vector< double >> &triangleStrips)
 
const Vec3D getAxis () const
 
BaseInteractiongetInteractionWith (BaseParticle *p, unsigned timeStamp, InteractionHandler *interactionHandler) override
 Returns the interaction between this wall and a given particle, nullptr if there is no interaction. More...
 
virtual void actionsOnRestart ()
 No implementation but can be overidden in its derived classes. More...
 
virtual void actionsAfterParticleGhostUpdate ()
 No implementation but can be overidden in its derived classes. More...
 
virtual void handleParticleAddition (unsigned int id, BaseParticle *p)
 Handles the addition of particles to the particleHandler. More...
 
virtual void handleParticleRemoval (unsigned int id)
 Handles the addition of particles to the particleHandler. More...
 
virtual void checkInteractions (InteractionHandler *interactionHandler, unsigned int timeStamp)
 Check if all interactions are valid. More...
 
bool getVTKVisibility () const
 
void setVTKVisibility (bool vtkVisibility)
 
void addRenderedWall (BaseWall *w)
 
BaseWallgetRenderedWall (size_t i) const
 
void removeRenderedWalls ()
 
void renderWall (VTKContainer &vtk)
 
void addParticlesAtWall (unsigned numElements=50)
 
void setVelocityControl (Vec3D forceGoal, Vec3D gainFactor, Vec3D baseVelocity)
 
- Public Member Functions inherited from BaseInteractable
 BaseInteractable ()
 Default BaseInteractable constructor. More...
 
 BaseInteractable (const BaseInteractable &p)
 Copy constructor. More...
 
 ~BaseInteractable () override
 Destructor, it simply destructs the BaseInteractable and all the objects it contains. More...
 
unsigned int getIndSpecies () const
 Returns the index of the species associated with the interactable object. More...
 
const ParticleSpeciesgetSpecies () const
 Returns a pointer to the species of this BaseInteractable. More...
 
void setSpecies (const ParticleSpecies *species)
 Sets the species of this BaseInteractable. More...
 
const Vec3DgetForce () const
 Returns the force on this BaseInteractable. More...
 
const Vec3DgetTorque () const
 Returns the torque on this BaseInteractable. More...
 
void setForce (const Vec3D &force)
 Sets the force on this BaseInteractable. More...
 
void setTorque (const Vec3D &torque)
 Sets the torque on this BaseInteractable. More...
 
void addForce (const Vec3D &addForce)
 Adds an amount to the force on this BaseInteractable. More...
 
void addTorque (const Vec3D &addTorque)
 Adds an amount to the torque on this BaseInteractable. More...
 
void resetForceTorque (int numberOfOMPthreads)
 
void sumForceTorqueOMP ()
 
const Vec3DgetPosition () const
 Returns the position of this BaseInteractable. More...
 
const QuaterniongetOrientation () const
 Returns the orientation of this BaseInteractable. More...
 
void setPosition (const Vec3D &position)
 Sets the position of this BaseInteractable. More...
 
void setOrientationViaNormal (Vec3D normal)
 Sets the orientation of this BaseInteractable by defining the vector that results from the rotation of the (1,0,0) vector. More...
 
void setOrientationViaEuler (Vec3D eulerAngle)
 Sets the orientation of this BaseInteractable by defining the euler angles. More...
 
void setOrientation (const Quaternion &orientation)
 Sets the orientation of this BaseInteractable. More...
 
virtual void move (const Vec3D &move)
 Moves this BaseInteractable by adding an amount to the position. More...
 
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 Mdouble getInvMass () const
 
virtual Mdouble getCurvature (const Vec3D &labFixedCoordinates) 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 Attributes

Vec3D start_
 The centre of the lower end of the screw. More...
 
Mdouble l_
 The length of the Screw. More...
 
Mdouble maxR_
 The outer radius of the Screw. More...
 
Mdouble n_
 The number of revelations. More...
 
Mdouble omega_
 Rotation speed in rev/s. More...
 
Mdouble offset_
 The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation. More...
 
Mdouble thickness_
 The thickness of the Screw. More...
 
ScrewType screwType_
 Single or double helix screw. More...
 

Additional Inherited Members

- Static Public Member Functions inherited from BaseWall
static void addToVTK (const std::vector< Vec3D > &points, VTKContainer &vtk)
 Takes the points provided and adds a triangle strip connecting these points to the vtk container. More...
 

Detailed Description

This function defines an Archimedes' screw in the z-direction from a (constant) starting point, a (constant) length L, a (constant) radius r, a (constant) number or revelations N and a (constant) rotation speed (rev/s)

q is a new coordinate going from 0 to 1 and t is the time, x=xs+r*cos(2*pi*(offset+N*q)), y=ys+r*sin(2*pi*(offset+N*q)), z=zs+q*L

Todo:
IFCD: Can these details about class Screw be made more clear? I don't understand them.

Definition at line 43 of file Screw.h.

Constructor & Destructor Documentation

Screw::Screw ( )

Default constructor: make a screw with default parameters.

Make a Screw which is centered in the origin, has a length of 1, one revelation, a radius of 1, turns with 1 revelation per second, is infinitely thin and starts at its normal initial point.

Definition at line 41 of file Screw.cc.

References DEBUG, doubleHelix, l_, logger, maxR_, n_, offset_, omega_, screwType_, BaseInteractable::setOrientationViaNormal(), Vec3D::setZero(), start_, and thickness_.

Referenced by copy().

42 {
43  start_.setZero();
44  l_ = 1.0;
45  maxR_ = 1.0;
46  n_ = 1.0;
47  omega_ = 1.0;
48  offset_ = 0.0;
49  thickness_ = 0.0;
51  setOrientationViaNormal({0, 0, 1});//default screw is in z-direction
52  logger(DEBUG, "Screw() constructor finished.");
53 }
Mdouble maxR_
The outer radius of the Screw.
Definition: Screw.h:124
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
void setZero()
Sets all elements to zero.
Definition: Vector.cc:43
ScrewType screwType_
Single or double helix screw.
Definition: Screw.h:144
Vec3D start_
The centre of the lower end of the screw.
Definition: Screw.h:116
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:136
Mdouble l_
The length of the Screw.
Definition: Screw.h:120
Mdouble n_
The number of revelations.
Definition: Screw.h:128
void setOrientationViaNormal(Vec3D normal)
Sets the orientation of this BaseInteractable by defining the vector that results from the rotation o...
Mdouble thickness_
The thickness of the Screw.
Definition: Screw.h:140
Mdouble omega_
Rotation speed in rev/s.
Definition: Screw.h:132
Screw::Screw ( const Screw other)

Copy constructor, copies another Screw.

Parameters
[in]otherThe Screw that has to be copied.

Definition at line 58 of file Screw.cc.

References DEBUG, l_, logger, maxR_, n_, offset_, omega_, screwType_, start_, and thickness_.

59  : BaseWall(other)
60 {
61  start_ = other.start_;
62  l_ = other.l_;
63  maxR_ = other.maxR_;
64  n_ = other.n_;
65  omega_ = other.omega_;
66  thickness_ = other.thickness_;
67  offset_ = other.offset_;
68  screwType_ = other.screwType_;
69  logger(DEBUG, "Screw(const Screw&) copy constructor finished.");
70 }
Mdouble maxR_
The outer radius of the Screw.
Definition: Screw.h:124
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
ScrewType screwType_
Single or double helix screw.
Definition: Screw.h:144
Vec3D start_
The centre of the lower end of the screw.
Definition: Screw.h:116
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:136
Mdouble l_
The length of the Screw.
Definition: Screw.h:120
BaseWall()
Default constructor.
Definition: BaseWall.cc:36
Mdouble n_
The number of revelations.
Definition: Screw.h:128
Mdouble thickness_
The thickness of the Screw.
Definition: Screw.h:140
Mdouble omega_
Rotation speed in rev/s.
Definition: Screw.h:132
Screw::Screw ( Vec3D  start,
Mdouble  l,
Mdouble  r,
Mdouble  n,
Mdouble  omega,
Mdouble  thickness,
ScrewType  screwType = ScrewType::doubleHelix 
)

Constructor in which all parameters of the screw are set.

Parameters
[in]startA Vec3D which denotes the centre of the lower end of the Screw.
[in]lThe length of the Screw, must be positive.
[in]rThe radius of the Screw, must be positive.
[in]nThe number of revelations of the Screw, must be positive.
[in]omegaThe rotation speed of the Screw in rev/s.
[in]thicknessThe thickness of the Screw, must be non-negative.

Make a Screw by assigning all input parameters to the data-members of this class, and setting the offset_ to 0.

Definition at line 82 of file Screw.cc.

References DEBUG, l_, logger, maxR_, n_, offset_, omega_, screwType_, start_, and thickness_.

83 {
84  start_ = start;
85  l_ = l;
86  maxR_ = r;
87  n_ = n;
88  omega_ = omega;
89  thickness_ = thickness;
90  offset_ = 0.0;
91  screwType_ = screwType;
92  logger(DEBUG, "Screw(Vec3D, Mdouble, Mdouble, Mdouble, Mdouble, Mdouble) constructor finished.");
93 }
Mdouble maxR_
The outer radius of the Screw.
Definition: Screw.h:124
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
ScrewType screwType_
Single or double helix screw.
Definition: Screw.h:144
Vec3D start_
The centre of the lower end of the screw.
Definition: Screw.h:116
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:136
Mdouble l_
The length of the Screw.
Definition: Screw.h:120
Mdouble n_
The number of revelations.
Definition: Screw.h:128
Mdouble thickness_
The thickness of the Screw.
Definition: Screw.h:140
Mdouble omega_
Rotation speed in rev/s.
Definition: Screw.h:132
Screw::~Screw ( )
override

Default destructor.

Definition at line 95 of file Screw.cc.

References DEBUG, and logger.

96 {
97  logger(DEBUG, "~Screw() finished, destroyed the Screw.");
98 }
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...

Member Function Documentation

Screw * Screw::copy ( ) const
finalvirtual

Copy this screw and return a pointer to the copy.

Returns
A pointer to a copy of this Screw.

Implements BaseWall.

Definition at line 103 of file Screw.cc.

References Screw().

104 {
105  return new Screw(*this);
106 }
Screw()
Default constructor: make a screw with default parameters.
Definition: Screw.cc:41
bool Screw::getDistanceAndNormal ( const BaseParticle P,
Mdouble distance,
Vec3D normal_return 
) const
finalvirtual

Compute the distance from the Screw for a given BaseParticle and return if there is a collision. If there is a collision, also return the normal vector of the interaction point.

Todo:
do this for all walls

Implements BaseWall.

Definition at line 108 of file Screw.cc.

References getDistanceAndNormalLabCoordinates(), BaseHandler< T >::getDPMBase(), BaseWall::getHandler(), BaseSpecies::getInteractionDistance(), SpeciesHandler::getMixedObject(), BaseInteractable::getOrientation(), BaseInteractable::getPosition(), BaseParticle::getRadius(), BaseInteractable::getSpecies(), Quaternion::rotate(), Quaternion::rotateBack(), and DPMBase::speciesHandler.

109 {
110  //transform coordinates into position-orientation frame
111  Vec3D position = P.getPosition() - getPosition();
112  getOrientation().rotateBack(position);
115  if (getDistanceAndNormalLabCoordinates(position, P.getRadius() + s->getInteractionDistance(), distance,
116  normal_return))
117  {
118  getOrientation().rotate(normal_return);
119  return true;
120  }
121  else
122  {
123  return false;
124  }
125 }
std::enable_if<!std::is_pointer< typename U::MixedSpeciesType >::value, typename U::MixedSpeciesType * >::type getMixedObject(const U *S, const U *T)
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
BaseSpecies is the class from which all other species are derived.
Definition: BaseSpecies.h:49
void rotateBack(Vec3D &position) const
Applies the inverse rotation to a position.
Definition: Quaternion.cc:592
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
void rotate(Vec3D &position) const
Applies the rotation to a position.
Definition: Quaternion.cc:563
Mdouble getInteractionDistance() const
returns the largest separation distance at which adhesive short-range forces can occur.
Definition: BaseSpecies.h:146
WallHandler * getHandler() const
A function which returns the WallHandler that handles this BaseWall.
Definition: BaseWall.cc:136
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:348
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1385
bool getDistanceAndNormalLabCoordinates(Vec3D position, Mdouble wallInteractionRadius, Mdouble &distance, Vec3D &normal_return) const
Definition: Screw.cc:138
Definition: Vector.h:49
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
bool Screw::getDistanceAndNormalLabCoordinates ( Vec3D  position,
Mdouble  wallInteractionRadius,
Mdouble distance,
Vec3D normal_return 
) const
Parameters
[in]pBaseParticle we want to calculate the distance and whether it collided of.
[out]distanceThe distance of the BaseParticle to this wall.
[out]normal_returnIf there was a collision, the normal vector to this wall will be placed here.
Returns
A boolean which says whether or not there was a collision.

This function computes whether or not there is a collision between a given BaseParticle and this Screw. If there is a collision, this function also computes the distance between the BaseParticle and Screw and the normal of the IntersectionOfWalls at the intersection point.

Todo:
Make this function readable and explain the steps in the details.

Definition at line 138 of file Screw.cc.

References A, mathsFunc::cos(), BaseObject::getId(), l_, logger, maxR_, n_, Vec3D::normalise(), offset_, constants::pi, R, helpers::round(), screwType_, mathsFunc::sign(), mathsFunc::sin(), singleHelix, constants::sqr_pi, mathsFunc::square(), start_, thickness_, WARN, X, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getDistanceAndNormal().

140 {
141  // compute the square of the radial distance (in yz-plane) between particle and screw.
142  const Mdouble RSquared = square(position.Y - start_.Y) + square(position.Z - start_.Z);
143  const Mdouble X = position.X - start_.X;
144 
145  //first do a simple check if particle is within the cylindrical hull of the screw
146  if (RSquared > square(maxR_ + wallInteractionRadius + thickness_)
147  || X > l_ + wallInteractionRadius + thickness_
148  || X < -wallInteractionRadius - thickness_)
149  {
150  return false;
151  }
152 
153  //else compute radial distance, and angle in yz-plane
154  const Mdouble R = sqrt(RSquared);
155  const Mdouble A = atan2(position.Z - start_.Z, position.Y - start_.Y);
156 
157  //after subtracting the start position and transforming the particle position from (XYZ) into (XRA)
158  //coordinates, we compute the distance to the wall at, located at (r,a,z)=(q*L,r,2*pi*(offset+N*q+k/2)), 0<q<1.
159 
160  //To find the contact point we have to minimize (with respect to r and q)
161  //distance^2=(x-x0-r*cos(2*pi*(offset+N*q)))^2+(y-y0-r*sin(2*pi*(offset+N*q)))^2+(z-z0-q*L)^2
162  //Using polar coordinates (i.e. x-x0=R*cos(A), y-y0=R*sin(A) and Z=z-z0)
163  //distance^2=R^2+r^2-2*R*r*cos(A-2*pi*(offset+N*q))+(Z-q*L)^2
164 
165  //Assumption: d(distance)/dr=0 at minDistance (should there be also a q-derivative?)
166  //Differentiate with respect to r and solve for zero:
167  //0=2*r-2*R*cos(A-2*pi*(offset+N*q)
168  //r=R*cos(A-2*pi*(offset+N*q))
169 
170  //Substitue back
171  //distance^2=R^2+R^2*cos^2(A-2*pi*(offset+N*q))-2*R^2*cos^2(A-2*pi*(offset+N*q))+(Z-q*L)^2
172  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))+(Z-q*L)^2
173 
174  //So we have to minimize:
175  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))^2 + (Z-q*L)^2 = f(q)
176  //f'(q)=(-2*pi*N)*R^2*sin(2*A-4*pi*(offset+N*q)) + 2*L*(Z-q*L) (D[Sin[x]^2,x]=Sin[2x])
177  //f''(q)=(4*pi^2*N^2)*R^2*cos(2*A-4*pi*(offset+N*q)) - 2*L*L
178  //For this we use the Euler algoritm
179 
180  Mdouble q; //Current guess
181  Mdouble dd; //Derivative at current guess
182  Mdouble ddd; //Second derivative at current guess
183 
184  //Set initial q to the closest position on the screw with the same angle as A
185  //The initial guess will be in the minimum of the sin closest to q0
186  //Minima of the sin are at
187  //A-2*pi*(offset+N*q)=k*pi (k=integer)
188  //q=A/(2*pi*N)-k/(2*N)-offset/N (k=integer)
189  {
190  const Mdouble q0 = X / l_; //assume closest point q0 is at same z-location as the particle
191  const Mdouble k = round(A / constants::pi - 2.0 * (offset_ + n_ * q0)); // k: |A-a(q0)-k*pi|=min
192  q = A / (2.0 * constants::pi * n_) - k / (2.0 * n_) - offset_ / n_; // q: a(q)=A
193  }
194 
195  //Now apply Newton's method
196  unsigned count = 0;
197  const unsigned maxCount = 30;
198  do
199  {
200  Mdouble arg = 2.0 * A - 4.0 * constants::pi * (n_ * q + offset_);
201  dd = -2.0 * constants::pi * n_ * RSquared * sin(arg) - 2.0 * l_ * (X - q * l_);
202  ddd = 8.0 * constants::sqr_pi * n_ * n_ * RSquared * cos(arg) + 2.0 * l_ * l_;
203  q -= dd / ddd;
204  if(++count>maxCount) {
205  logger(WARN,"Screw % contact detection did not converge (err=%); continuing with approximate contact point",getId(),fabs(dd/ddd));
206  break;
207  }
208  } while (fabs(dd / ddd) > 5e-14);
209 
210  //Calculate r
211  Mdouble r = R * cos(2.0 * constants::pi * (offset_ + n_ * q) - A);
212 
213  //If the particle touches the single screw center
214  if (screwType_==ScrewType::singleHelix && r>0) {
215  distance = R-thickness_;
216  if (distance>wallInteractionRadius) {
217  return false;
218  } else {
219  normal_return = Vec3D(0,-position.Y/R,-position.Z/R);
220  return true;
221  }
222  }
223 
224  //Check if the location is actually on the screw:
225  //First posibility is that the radius is too large:
226  Mdouble distanceSquared = 0;
227  if (fabs(r) > maxR_) //Left boundary of the coil
228  {
229  r = mathsFunc::sign(r) * maxR_;
230  distanceSquared = mathsFunc::square(R-maxR_);
231  count = 0;
232  //This case reduces to the coil problem
233  do
234  {
235  dd = -4.0 * R * r * constants::pi * n_ * sin(A - 2.0 * constants::pi * (n_ * q + offset_)) -
236  2.0 * l_ * (X - q * l_);
237  ddd = 8.0 * R * r * constants::sqr_pi * n_ * n_ * cos(A - 2.0 * constants::pi * (n_ * q + offset_)) +
238  2.0 * l_ * l_;
239  q -= dd / ddd;
240  if(++count>maxCount) {
241  logger(WARN,"Screw % contact detection at left boundary did not converge (err=%); continuing with approximate contact point",getId(),fabs(dd/ddd));
242  break;
243  }
244  } while (fabs(dd / ddd) > 1e-13);
245  }
246  //Second possibility is that it occured before the start of after the end
247  if (q < 0)
248  {
249  q = 0;
250  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
251  if (fabs(r) > maxR_)
252  {
253  r = mathsFunc::sign(r) * maxR_;
254  }
255  }
256  else if (q > 1)
257  {
258  q = 1;
259  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
260  if (fabs(r) > maxR_)
261  {
262  r = mathsFunc::sign(r) * maxR_;
263  }
264  }
265 
266  //compute squared distance between particle and contact point
267  distanceSquared += square(X - q * l_) + square(R*sin(A - 2 * constants::pi * (offset_ + n_ * q)));
268 
269  //If distance is too large there is no contact
270  if (distanceSquared >= square(wallInteractionRadius + thickness_))
271  {
272  return false;
273  }
274 
275  //Else
276  Vec3D ContactPoint;
277  distance = sqrt(distanceSquared) - thickness_;
278  ContactPoint.Y = start_.Y + r * cos(2.0 * constants::pi * (offset_ + n_ * q));
279  ContactPoint.Z = start_.Z + r * sin(2.0 * constants::pi * (offset_ + n_ * q));
280  ContactPoint.X = start_.X + q * l_;
281  normal_return = (ContactPoint - position);
282  normal_return.normalise();
283  return true;
284 }
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:125
Mdouble X
the vector components
Definition: Vector.h:65
Mdouble maxR_
The outer radius of the Screw.
Definition: Screw.h:124
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
double Mdouble
Definition: GeneralDefine.h:34
void normalise()
Makes this Vec3D unit length.
Definition: Vector.cc:123
int sign(T val)
This is a sign function, it returns -1 for negative numbers, 1 for positive numbers and 0 for 0...
Definition: ExtendedMath.h:97
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
ScrewType screwType_
Single or double helix screw.
Definition: Screw.h:144
Vec3D start_
The centre of the lower end of the screw.
Definition: Screw.h:116
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:136
const Mdouble pi
Definition: ExtendedMath.h:45
Mdouble l_
The length of the Screw.
Definition: Screw.h:120
Mdouble Y
Definition: Vector.h:65
Mdouble round(const Mdouble value, unsigned precision)
Definition: Helpers.cc:598
Mdouble n_
The number of revelations.
Definition: Screw.h:128
Mdouble thickness_
The thickness of the Screw.
Definition: Screw.h:140
Definition: Vector.h:49
T square(const T val)
squares a number
Definition: ExtendedMath.h:106
Mdouble Z
Definition: Vector.h:65
const Mdouble sqr_pi
Definition: ExtendedMath.h:47
std::string Screw::getName ( ) const
finalvirtual

Returns the name of the object, here the string "Screw".

Returns
The string "Screw".

Implements BaseObject.

Definition at line 360 of file Screw.cc.

Referenced by writeVTK().

361 {
362  return "Screw";
363 }
void Screw::move_time ( Mdouble  dt)

Rotate the Screw for a period dt, so that the offset_ changes with omega_*dt.

Parameters
[in]dtThe time for which the Screw has to be turned.

Definition at line 289 of file Screw.cc.

290 {
291  //offset_ += omega_ * dt;
292 }
void Screw::oldRead ( std::istream &  is)

Reads a Screw in the old style from an input stream, for example a restart file old style.

Parameters
[in,out]isInput stream from which the Screw must be read.

Read the Screw in old style, please note that the thickness is not read in this function, so it has either to be set manually or it is 0.0 from the default constructor.

Definition at line 330 of file Screw.cc.

References l_, maxR_, n_, offset_, omega_, and start_.

331 {
332  std::string dummy;
333  is >> dummy >> start_
334  >> dummy >> l_
335  >> dummy >> maxR_
336  >> dummy >> n_
337  >> dummy >> omega_
338  >> dummy >> offset_;
339 }
Mdouble maxR_
The outer radius of the Screw.
Definition: Screw.h:124
Vec3D start_
The centre of the lower end of the screw.
Definition: Screw.h:116
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:136
Mdouble l_
The length of the Screw.
Definition: Screw.h:120
Mdouble n_
The number of revelations.
Definition: Screw.h:128
Mdouble omega_
Rotation speed in rev/s.
Definition: Screw.h:132
void Screw::read ( std::istream &  is)
overridevirtual

Reads a Screw from an input stream, for example a restart file.

Parameters
[in,out]isInput stream from which the Screw must be read.

Reimplemented from BaseWall.

Definition at line 308 of file Screw.cc.

References l_, maxR_, n_, offset_, omega_, BaseWall::read(), screwType_, start_, and thickness_.

309 {
310  BaseWall::read(is);
311  std::string dummy;
312  unsigned screwType;
313  is >> dummy >> start_
314  >> dummy >> l_
315  >> dummy >> maxR_
316  >> dummy >> n_
317  >> dummy >> omega_
318  >> dummy >> thickness_
319  >> dummy >> offset_
320  >> dummy >> screwType;
321  screwType_ = static_cast<ScrewType>(screwType);
322 }
Mdouble maxR_
The outer radius of the Screw.
Definition: Screw.h:124
ScrewType screwType_
Single or double helix screw.
Definition: Screw.h:144
Vec3D start_
The centre of the lower end of the screw.
Definition: Screw.h:116
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:136
Mdouble l_
The length of the Screw.
Definition: Screw.h:120
void read(std::istream &is) override
Function that reads a BaseWall from an input stream, usually a restart file.
Definition: BaseWall.cc:80
Mdouble n_
The number of revelations.
Definition: Screw.h:128
ScrewType
used to define if the screw consists of a single or double helical thread
Definition: Screw.h:35
Mdouble thickness_
The thickness of the Screw.
Definition: Screw.h:140
Mdouble omega_
Rotation speed in rev/s.
Definition: Screw.h:132
void Screw::rotate ( const Vec3D angularVelocityDt)
overridevirtual
Todo:
the move and rotate functions should only pass the time step, as teh velocity can be accessed directly by the object
Parameters
angularVelocityDt

Reimplemented from BaseInteractable.

Definition at line 298 of file Screw.cc.

References BaseHandler< T >::getDPMBase(), BaseWall::getHandler(), DPMBase::getTimeStep(), offset_, omega_, and BaseInteractable::rotate().

299 {
300  BaseInteractable::rotate(angularVelocityDt);
302 }
WallHandler * getHandler() const
A function which returns the WallHandler that handles this BaseWall.
Definition: BaseWall.cc:136
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:136
virtual void rotate(const Vec3D &angularVelocityDt)
Rotates this BaseInteractable.
Mdouble omega_
Rotation speed in rev/s.
Definition: Screw.h:132
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
Mdouble getTimeStep() const
Returns the simulation time step.
Definition: DPMBase.cc:1241
void Screw::write ( std::ostream &  os) const
overridevirtual

Writes this Screw to an output stream, for example a restart file.

Parameters
[in,out]osOutput stream to which the Screw must be written.

Reimplemented from BaseWall.

Definition at line 344 of file Screw.cc.

References l_, maxR_, n_, offset_, omega_, screwType_, start_, thickness_, and BaseWall::write().

345 {
346  BaseWall::write(os);
347  os << " start " << start_
348  << " Length " << l_
349  << " Radius " << maxR_
350  << " Revolutions " << n_
351  << " Omega " << omega_
352  << " Thickness " << thickness_
353  << " Offset " << offset_
354  << " ScrewType " << static_cast<unsigned>(screwType_);
355 }
Mdouble maxR_
The outer radius of the Screw.
Definition: Screw.h:124
void write(std::ostream &os) const override
Function that writes a BaseWall to an output stream, usually a restart file.
Definition: BaseWall.cc:102
ScrewType screwType_
Single or double helix screw.
Definition: Screw.h:144
Vec3D start_
The centre of the lower end of the screw.
Definition: Screw.h:116
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:136
Mdouble l_
The length of the Screw.
Definition: Screw.h:120
Mdouble n_
The number of revelations.
Definition: Screw.h:128
Mdouble thickness_
The thickness of the Screw.
Definition: Screw.h:140
Mdouble omega_
Rotation speed in rev/s.
Definition: Screw.h:132
void Screw::writeVTK ( VTKContainer vtk) const
overridevirtual

adds extra information to the points and triangleStrips vectors needed to plot the wall in vtk format

Parameters
pointsCoordinates of the vertices of the triangulated surfaces (in the VTK file this is called POINTS)
triangleStripsIndices of three vertices forming one triangulated surface (in the VTK file this is called CELL)

Reimplemented from BaseWall.

Definition at line 365 of file Screw.cc.

References mathsFunc::cos(), doubleHelix, BaseInteractable::getOrientation(), BaseInteractable::getPosition(), l_, maxR_, n_, offset_, constants::pi, VTKContainer::points, Quaternion::rotate(), screwType_, mathsFunc::sin(), start_, thickness_, VTKContainer::triangleStrips, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by writeVTK().

366 {
367  //number of points in radial direction (for one side of the screw)
368  unsigned nr = 5;
369  //number of points in axial direction
370  unsigned nz = static_cast<unsigned int>(99 * fabs(n_));
371 
372  unsigned long nPoints = vtk.points.size();
373  Vec3D contactPoint;
374  // either one or two helices
375  for (Mdouble offset = offset_; offset<=offset_+(screwType_==ScrewType::doubleHelix?0.5:0); offset+=0.5)
376  {
377  for (unsigned iz = 0; iz < nz; iz++)
378  {
379  for (unsigned ir = 0; ir < nr; ir++)
380  {
381  double q = (double) iz / nz;
382  double r = (double) ir / (nr - 1) * maxR_;
383  contactPoint.Y = start_.Y - r * cos(2.0 * constants::pi * (offset + n_ * q));
384  contactPoint.Z = start_.Z - r * sin(2.0 * constants::pi * (offset + n_ * q));
385  contactPoint.X = start_.X + q * l_ - thickness_;
386  getOrientation().rotate(contactPoint);
387  contactPoint += getPosition();
388  vtk.points.push_back(contactPoint);
389  }
390  for (unsigned ir = 0; ir < nr; ir++)
391  {
392  double q = (double) iz / nz;
393  double r = (double) (nr - 1 - ir) / (nr - 1) * maxR_;
394  contactPoint.Y = start_.Y - r * cos(2.0 * constants::pi * (offset + n_ * q));
395  contactPoint.Z = start_.Z - r * sin(2.0 * constants::pi * (offset + n_ * q));
396  contactPoint.X = start_.X + q * l_ + thickness_;
397  getOrientation().rotate(contactPoint);
398  contactPoint += getPosition();
399  vtk.points.push_back(contactPoint);
400  }
401  }
402  }
403 
404  unsigned long nCells = vtk.triangleStrips.size();
405  //vtk.triangleStrips.reserve(nCells + (nz - 1)*(screwType_==ScrewType::doubleHelix?2:1));
406  for (unsigned iz = 0; iz < (screwType_==ScrewType::doubleHelix?2:1)*nz-1; iz++)
407  {
408  //skip step that would connect the two screw parts
409  if (iz==nz-1) ++iz;
410  std::vector<double> cell;
411  cell.reserve(2 * nr);
412  for (unsigned ir = 0; ir < 2*nr; ir++)
413  {
414  cell.push_back(nPoints + ir + 2*iz * nr);
415  cell.push_back(nPoints + ir + 2*(iz + 1) * nr);
416  }
417  vtk.triangleStrips.push_back(cell);
418  }
419 }
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
Mdouble maxR_
The outer radius of the Screw.
Definition: Screw.h:124
double Mdouble
Definition: GeneralDefine.h:34
void rotate(Vec3D &position) const
Applies the rotation to a position.
Definition: Quaternion.cc:563
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
ScrewType screwType_
Single or double helix screw.
Definition: Screw.h:144
Vec3D start_
The centre of the lower end of the screw.
Definition: Screw.h:116
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:136
const Mdouble pi
Definition: ExtendedMath.h:45
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:39
Mdouble l_
The length of the Screw.
Definition: Screw.h:120
Mdouble Y
Definition: Vector.h:65
std::vector< Vec3D > points
Definition: BaseWall.h:38
Mdouble n_
The number of revelations.
Definition: Screw.h:128
Mdouble thickness_
The thickness of the Screw.
Definition: Screw.h:140
Definition: Vector.h:49
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Mdouble Z
Definition: Vector.h:65
void Screw::writeVTK ( std::string  filename) const

Definition at line 421 of file Screw.cc.

References getName(), VTKContainer::points, VTKContainer::triangleStrips, helpers::writeToFile(), and writeVTK().

422 {
423  VTKContainer vtk;
424  writeVTK(vtk);
425 
426  std::stringstream file;
427  file << "# vtk DataFile Version 2.0\n"
428  << getName() << "\n"
429  "ASCII\n"
430  "DATASET UNSTRUCTURED_GRID\n"
431  "POINTS " << vtk.points.size() << " double\n";
432  for (const auto& vertex : vtk.points)
433  file << vertex << '\n';
434  file << "\nCELLS " << vtk.triangleStrips.size() << ' ' << 4 * vtk.triangleStrips.size() << "\n";
435  for (const auto& face : vtk.triangleStrips)
436  file << "3 " << face[0] << ' ' << face[1] << ' ' << face[2] << '\n';
437  file << "\nCELL_TYPES " << vtk.triangleStrips.size() << "\n";
438  for (const auto& face : vtk.triangleStrips)
439  file << "5\n";
440  helpers::writeToFile(filename, file.str());
441 }
bool writeToFile(std::string filename, std::string filecontent)
Writes a string to a file.
Definition: Helpers.cc:445
void writeVTK(VTKContainer &vtk) const override
Definition: Screw.cc:365
std::string getName() const final
Returns the name of the object, here the string "Screw".
Definition: Screw.cc:360
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:39
std::vector< Vec3D > points
Definition: BaseWall.h:38

Member Data Documentation

Mdouble Screw::l_
private

The length of the Screw.

Definition at line 120 of file Screw.h.

Referenced by getDistanceAndNormalLabCoordinates(), oldRead(), read(), Screw(), write(), and writeVTK().

Mdouble Screw::maxR_
private

The outer radius of the Screw.

Definition at line 124 of file Screw.h.

Referenced by getDistanceAndNormalLabCoordinates(), oldRead(), read(), Screw(), write(), and writeVTK().

Mdouble Screw::n_
private

The number of revelations.

Definition at line 128 of file Screw.h.

Referenced by getDistanceAndNormalLabCoordinates(), oldRead(), read(), Screw(), write(), and writeVTK().

Mdouble Screw::offset_
private

The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.

Definition at line 136 of file Screw.h.

Referenced by getDistanceAndNormalLabCoordinates(), oldRead(), read(), rotate(), Screw(), write(), and writeVTK().

Mdouble Screw::omega_
private

Rotation speed in rev/s.

Definition at line 132 of file Screw.h.

Referenced by oldRead(), read(), rotate(), Screw(), and write().

ScrewType Screw::screwType_
private

Single or double helix screw.

Definition at line 144 of file Screw.h.

Referenced by getDistanceAndNormalLabCoordinates(), read(), Screw(), write(), and writeVTK().

Vec3D Screw::start_
private

The centre of the lower end of the screw.

Definition at line 116 of file Screw.h.

Referenced by getDistanceAndNormalLabCoordinates(), oldRead(), read(), Screw(), write(), and writeVTK().

Mdouble Screw::thickness_
private

The thickness of the Screw.

Definition at line 140 of file Screw.h.

Referenced by getDistanceAndNormalLabCoordinates(), read(), Screw(), write(), and writeVTK().


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