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
 
std::vector< BaseWall * > getRenderedWalls () const
 
void removeRenderedWalls ()
 
void renderWall (VTKContainer &vtk)
 
void addParticlesAtWall (unsigned numElements=50)
 
void setVelocityControl (Vec3D forceGoal, Vec3D gainFactor, Vec3D baseVelocity)
 
virtual void writeWallDetailsVTK (VTKData &data) const
 
virtual void computeWear ()
 
- 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...
 
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
 
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 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.

Constructor & Destructor Documentation

◆ Screw() [1/3]

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.

21 {
22  start_.setZero();
23  l_ = 1.0;
24  maxR_ = 1.0;
25  n_ = 1.0;
26  omega_ = 1.0;
27  offset_ = 0.0;
28  thickness_ = 0.0;
30  setOrientationViaNormal({0, 0, 1});//default screw is in z-direction
31  logger(DEBUG, "Screw() constructor finished.");
32 }
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
@ DEBUG
void setOrientationViaNormal(Vec3D normal)
Sets the orientation of this BaseInteractable by defining the vector that results from the rotation o...
Definition: BaseInteractable.cc:177
Mdouble thickness_
The thickness of the Screw.
Definition: Screw.h:119
Vec3D start_
The centre of the lower end of the screw.
Definition: Screw.h:95
Mdouble l_
The length of the Screw.
Definition: Screw.h:99
Mdouble maxR_
The outer radius of the Screw.
Definition: Screw.h:103
Mdouble n_
The number of revelations.
Definition: Screw.h:107
ScrewType screwType_
Single or double helix screw.
Definition: Screw.h:123
Mdouble omega_
Rotation speed in rev/s.
Definition: Screw.h:111
Mdouble offset_
The angle that describes how much the Screw has turned, going from 0 to 1 for a rotation.
Definition: Screw.h:115
void setZero()
Sets all elements to zero.
Definition: Vector.cc:23

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

Referenced by copy().

◆ Screw() [2/3]

Screw::Screw ( const Screw other)

Copy constructor, copies another Screw.

Parameters
[in]otherThe Screw that has to be copied.
38  : BaseWall(other)
39 {
40  start_ = other.start_;
41  l_ = other.l_;
42  maxR_ = other.maxR_;
43  n_ = other.n_;
44  omega_ = other.omega_;
45  thickness_ = other.thickness_;
46  offset_ = other.offset_;
47  screwType_ = other.screwType_;
48  logger(DEBUG, "Screw(const Screw&) copy constructor finished.");
49 }
BaseWall()
Default constructor.
Definition: BaseWall.cc:15

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

◆ Screw() [3/3]

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.

62 {
63  start_ = start;
64  l_ = l;
65  maxR_ = r;
66  n_ = n;
67  omega_ = omega;
68  thickness_ = thickness;
69  offset_ = 0.0;
70  screwType_ = screwType;
71  logger(DEBUG, "Screw(Vec3D, Mdouble, Mdouble, Mdouble, Mdouble, Mdouble) constructor finished.");
72 }
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
Vector::Scalar omega(const Vector &t, const Vector &s, RealScalar angle)
Definition: IDRS.h:36
r
Definition: UniformPSDSelfTest.py:20
void start(const unsigned &i)
(Re-)start i-th timer
Definition: oomph_utilities.cc:243

References DEBUG, l_, logger, maxR_, n, n_, offset_, Eigen::internal::omega(), omega_, UniformPSDSelfTest::r, screwType_, oomph::CumulativeTimings::start(), start_, and thickness_.

◆ ~Screw()

Screw::~Screw ( )
override

Default destructor.

75 {
76  logger(DEBUG, "~Screw() finished, destroyed the Screw.");
77 }

References DEBUG, and logger.

Member Function Documentation

◆ copy()

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.

83 {
84  return new Screw(*this);
85 }
Screw()
Default constructor: make a screw with default parameters.
Definition: Screw.cc:20

References Screw().

◆ getDistanceAndNormal()

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.

88 {
89  //transform coordinates into position-orientation frame
90  Vec3D position = P.getPosition() - getPosition();
91  getOrientation().rotateBack(position);
94  if (getDistanceAndNormalLabCoordinates(position, P.getRadius() + s->getInteractionDistance(), distance,
95  normal_return))
96  {
97  getOrientation().rotate(normal_return);
98  return true;
99  }
100  else
101  {
102  return false;
103  }
104 }
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:733
const Quaternion & getOrientation() const
Returns the orientation of this BaseInteractable.
Definition: BaseInteractable.h:209
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
Definition: BaseInteractable.h:87
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Definition: BaseInteractable.h:197
BaseSpecies is the class from which all other species are derived.
Definition: BaseSpecies.h:29
WallHandler * getHandler() const
A function which returns the WallHandler that handles this BaseWall.
Definition: BaseWall.cc:113
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1433
void rotate(Vec3D &position) const
Definition: Quaternion.cc:550
void rotateBack(Vec3D &position) const
Definition: Quaternion.cc:597
bool getDistanceAndNormalLabCoordinates(Vec3D position, Mdouble wallInteractionRadius, Mdouble &distance, Vec3D &normal_return) const
Definition: Screw.cc:117
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
Definition: Kernel/Math/Vector.h:30
RealScalar s
Definition: level1_cplx_impl.h:130
double P
Uniform pressure.
Definition: TwenteMeshGluing.cpp:77

References getDistanceAndNormalLabCoordinates(), BaseHandler< T >::getDPMBase(), BaseWall::getHandler(), SpeciesHandler::getMixedObject(), BaseInteractable::getOrientation(), BaseInteractable::getPosition(), BaseInteractable::getSpecies(), Global_Physical_Variables::P, Quaternion::rotate(), Quaternion::rotateBack(), s, and DPMBase::speciesHandler.

◆ getDistanceAndNormalLabCoordinates()

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.
119 {
120  // compute the square of the radial distance (in yz-plane) between particle and screw.
121  const Mdouble RSquared = square(position.Y - start_.Y) + square(position.Z - start_.Z);
122  const Mdouble X = position.X - start_.X;
123 
124  //first do a simple check if particle is within the cylindrical hull of the screw
125  if (RSquared > square(maxR_ + wallInteractionRadius + thickness_)
126  || X > l_ + wallInteractionRadius + thickness_
127  || X < -wallInteractionRadius - thickness_)
128  {
129  return false;
130  }
131 
132  //else compute radial distance, and angle in yz-plane
133  const Mdouble R = sqrt(RSquared);
134  const Mdouble A = atan2(position.Z - start_.Z, position.Y - start_.Y);
135 
136  //after subtracting the start position and transforming the particle position from (XYZ) into (XRA)
137  //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.
138 
139  //To find the contact point we have to minimize (with respect to r and q)
140  //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
141  //Using polar coordinates (i.e. x-x0=R*cos(A), y-y0=R*sin(A) and Z=z-z0)
142  //distance^2=R^2+r^2-2*R*r*cos(A-2*pi*(offset+N*q))+(Z-q*L)^2
143 
144  //Assumption: d(distance)/dr=0 at minDistance (should there be also a q-derivative?)
145  //Differentiate with respect to r and solve for zero:
146  //0=2*r-2*R*cos(A-2*pi*(offset+N*q)
147  //r=R*cos(A-2*pi*(offset+N*q))
148 
149  //Substitue back
150  //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
151  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))+(Z-q*L)^2
152 
153  //So we have to minimize:
154  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))^2 + (Z-q*L)^2 = f(q)
155  //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])
156  //f''(q)=(4*pi^2*N^2)*R^2*cos(2*A-4*pi*(offset+N*q)) - 2*L*L
157  //For this we use the Euler algoritm
158 
159  Mdouble q; //Current guess
160  Mdouble dd; //Derivative at current guess
161  Mdouble ddd; //Second derivative at current guess
162 
163  //Set initial q to the closest position on the screw with the same angle as A
164  //The initial guess will be in the minimum of the sin closest to q0
165  //Minima of the sin are at
166  //A-2*pi*(offset+N*q)=k*pi (k=integer)
167  //q=A/(2*pi*N)-k/(2*N)-offset/N (k=integer)
168  {
169  const Mdouble q0 = X / l_; //assume closest point q0 is at same z-location as the particle
170  const Mdouble k = round(A / constants::pi - 2.0 * (offset_ + n_ * q0)); // k: |A-a(q0)-k*pi|=min
171  q = A / (2.0 * constants::pi * n_) - k / (2.0 * n_) - offset_ / n_; // q: a(q)=A
172  }
173 
174  //Now apply Newton's method
175  unsigned count = 0;
176  const unsigned maxCount = 30;
177  do
178  {
179  Mdouble arg = 2.0 * A - 4.0 * constants::pi * (n_ * q + offset_);
180  dd = -2.0 * constants::pi * n_ * RSquared * sin(arg) - 2.0 * l_ * (X - q * l_);
181  ddd = 8.0 * constants::sqr_pi * n_ * n_ * RSquared * cos(arg) + 2.0 * l_ * l_;
182  q -= dd / ddd;
183  if(++count>maxCount) {
184  logger(WARN,"Screw % contact detection did not converge (err=%); continuing with approximate contact point",getId(),fabs(dd/ddd));
185  break;
186  }
187  } while (fabs(dd / ddd) > 5e-14);
188 
189  //Calculate r
190  Mdouble r = R * cos(2.0 * constants::pi * (offset_ + n_ * q) - A);
191 
192  //If the particle touches the single screw center
193  if (screwType_==ScrewType::singleHelix && r>0) {
194  distance = R-thickness_;
195  if (distance>wallInteractionRadius) {
196  return false;
197  } else {
198  normal_return = Vec3D(0,-position.Y/R,-position.Z/R);
199  return true;
200  }
201  }
202 
203  //Check if the location is actually on the screw:
204  //First posibility is that the radius is too large:
205  Mdouble distanceSquared = 0;
206  if (fabs(r) > maxR_) //Left boundary of the coil
207  {
208  r = mathsFunc::sign(r) * maxR_;
209  distanceSquared = mathsFunc::square(R-maxR_);
210  count = 0;
211  //This case reduces to the coil problem
212  do
213  {
214  dd = -4.0 * R * r * constants::pi * n_ * sin(A - 2.0 * constants::pi * (n_ * q + offset_)) -
215  2.0 * l_ * (X - q * l_);
216  ddd = 8.0 * R * r * constants::sqr_pi * n_ * n_ * cos(A - 2.0 * constants::pi * (n_ * q + offset_)) +
217  2.0 * l_ * l_;
218  q -= dd / ddd;
219  if(++count>maxCount) {
220  logger(WARN,"Screw % contact detection at left boundary did not converge (err=%); continuing with approximate contact point",getId(),fabs(dd/ddd));
221  break;
222  }
223  } while (fabs(dd / ddd) > 1e-13);
224  }
225  //Second possibility is that it occured before the start of after the end
226  if (q < 0)
227  {
228  q = 0;
229  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
230  if (fabs(r) > maxR_)
231  {
232  r = mathsFunc::sign(r) * maxR_;
233  }
234  }
235  else if (q > 1)
236  {
237  q = 1;
238  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
239  if (fabs(r) > maxR_)
240  {
241  r = mathsFunc::sign(r) * maxR_;
242  }
243  }
244 
245  //compute squared distance between particle and contact point
246  distanceSquared += square(X - q * l_) + square(R*sin(A - 2 * constants::pi * (offset_ + n_ * q)));
247 
248  //If distance is too large there is no contact
249  if (distanceSquared >= square(wallInteractionRadius + thickness_))
250  {
251  return false;
252  }
253 
254  //Else
255  Vec3D ContactPoint;
256  distance = sqrt(distanceSquared) - thickness_;
257  ContactPoint.Y = start_.Y + r * cos(2.0 * constants::pi * (offset_ + n_ * q));
258  ContactPoint.Z = start_.Z + r * sin(2.0 * constants::pi * (offset_ + n_ * q));
259  ContactPoint.X = start_.X + q * l_;
260  normal_return = (ContactPoint - position);
261  normal_return.normalise();
262  return true;
263 }
AnnoyingScalar cos(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:136
AnnoyingScalar atan2(const AnnoyingScalar &y, const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:139
AnnoyingScalar sin(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:137
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
Array< double, 1, 3 > e(1./3., 0.5, 2.)
@ WARN
@ R
Definition: StatisticsVector.h:21
unsigned int getId() const
Returns the unique identifier of any particular object.
Definition: BaseObject.h:104
The matrix class, also used for vectors and row-vectors.
Definition: Eigen/Eigen/src/Core/Matrix.h:186
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
void normalise()
Makes this Vec3D unit length.
Definition: Vector.cc:103
#define X
Definition: icosphere.cpp:20
char char char int int * k
Definition: level2_impl.h:374
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 round(const bfloat16 &a)
Definition: BFloat16.h:646
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:2019
squared absolute sa ArrayBase::abs2 DOXCOMMA MatrixBase::cwiseAbs2 square(power 2)
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
const Mdouble sqr_pi
Definition: ExtendedMath.h:25
const Mdouble pi
Definition: ExtendedMath.h:23
T square(const T val)
squares a number
Definition: ExtendedMath.h:86
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 atan2(), cos(), e(), boost::multiprecision::fabs(), BaseObject::getId(), k, l_, logger, maxR_, n_, Vec3D::normalise(), offset_, constants::pi, Eigen::numext::q, UniformPSDSelfTest::r, R, Eigen::bfloat16_impl::round(), screwType_, mathsFunc::sign(), sin(), singleHelix, constants::sqr_pi, sqrt(), mathsFunc::square(), Eigen::square(), start_, thickness_, WARN, X, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by getDistanceAndNormal().

◆ getName()

std::string Screw::getName ( ) const
finalvirtual

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

Returns
The string "Screw".

Implements BaseObject.

340 {
341  return "Screw";
342 }

Referenced by writeVTK().

◆ move_time()

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.
269 {
270  //offset_ += omega_ * dt;
271 }

◆ oldRead()

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.

310 {
311  std::string dummy;
312  is >> dummy >> start_
313  >> dummy >> l_
314  >> dummy >> maxR_
315  >> dummy >> n_
316  >> dummy >> omega_
317  >> dummy >> offset_;
318 }
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286

References l_, maxR_, n_, offset_, omega_, start_, and oomph::Global_string_for_annotation::string().

◆ read()

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.

288 {
289  BaseWall::read(is);
290  std::string dummy;
291  unsigned screwType;
292  is >> dummy >> start_
293  >> dummy >> l_
294  >> dummy >> maxR_
295  >> dummy >> n_
296  >> dummy >> omega_
297  >> dummy >> thickness_
298  >> dummy >> offset_
299  >> dummy >> screwType;
300  screwType_ = static_cast<ScrewType>(screwType);
301 }
ScrewType
Definition: Screw.h:14
void read(std::istream &is) override
Function that reads a BaseWall from an input stream, usually a restart file.
Definition: BaseWall.cc:57

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

◆ rotate()

void Screw::rotate ( const Vec3D angularVelocityDt)
overridevirtual
Todo:
the move and rotate functions should only pass the time step, as the velocity can be accessed directly by the object
Parameters
angularVelocityDt

Reimplemented from BaseInteractable.

278 {
279  BaseInteractable::rotate(angularVelocityDt);
281 }
virtual void rotate(const Vec3D &angularVelocityDt)
Rotates this BaseInteractable.
Definition: BaseInteractable.cc:208
Mdouble getTimeStep() const
Returns the simulation time step.
Definition: DPMBase.cc:1241

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

◆ write()

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.

324 {
325  BaseWall::write(os);
326  os << " start " << start_
327  << " Length " << l_
328  << " Radius " << maxR_
329  << " Revolutions " << n_
330  << " Omega " << omega_
331  << " Thickness " << thickness_
332  << " Offset " << offset_
333  << " ScrewType " << static_cast<unsigned>(screwType_);
334 }
void write(std::ostream &os) const override
Function that writes a BaseWall to an output stream, usually a restart file.
Definition: BaseWall.cc:79

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

◆ writeVTK() [1/2]

void Screw::writeVTK ( std::string  filename) const
399 {
400  VTKContainer vtk;
401  writeVTK(vtk);
402 
403  std::stringstream file;
404  file << "# vtk DataFile Version 2.0\n"
405  << getName() << "\n"
406  "ASCII\n"
407  "DATASET UNSTRUCTURED_GRID\n"
408  "POINTS " << vtk.points.size() << " double\n";
409  for (const auto& vertex : vtk.points)
410  file << vertex << '\n';
411  file << "\nCELLS " << vtk.triangleStrips.size() << ' ' << 4 * vtk.triangleStrips.size() << "\n";
412  for (const auto& face : vtk.triangleStrips)
413  file << "3 " << face[0] << ' ' << face[1] << ' ' << face[2] << '\n';
414  file << "\nCELL_TYPES " << vtk.triangleStrips.size() << "\n";
415  for (const auto& face UNUSED : vtk.triangleStrips)
416  file << "5\n";
417  helpers::writeToFile(filename, file.str());
418 }
#define UNUSED
Definition: GeneralDefine.h:18
void writeVTK(VTKContainer &vtk) const override
Definition: Screw.cc:344
std::string getName() const final
Returns the name of the object, here the string "Screw".
Definition: Screw.cc:339
string filename
Definition: MergeRestartFiles.py:39
bool writeToFile(const std::string &filename, const std::string &filecontent)
Writes a string to a file.
Definition: FileIOHelpers.cc:29
Definition: BaseWall.h:17
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:19
std::vector< Vec3D > points
Definition: BaseWall.h:18

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

◆ writeVTK() [2/2]

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.

345 {
346  //number of points in radial direction (for one side of the screw)
347  unsigned nr = 5;
348  //number of points in axial direction
349  unsigned nz = static_cast<unsigned int>(99 * fabs(n_));
350 
351  unsigned long nPoints = vtk.points.size();
352  Vec3D contactPoint;
353  // either one or two helices
354  for (Mdouble offset = offset_; offset<=offset_+(screwType_==ScrewType::doubleHelix?0.5:0); offset+=0.5)
355  {
356  for (unsigned iz = 0; iz < nz; iz++)
357  {
358  for (unsigned ir = 0; ir < nr; ir++)
359  {
360  double q = (double) iz / nz;
361  double r = (double) ir / (nr - 1) * maxR_;
362  contactPoint.Y = start_.Y - r * cos(2.0 * constants::pi * (offset + n_ * q));
363  contactPoint.Z = start_.Z - r * sin(2.0 * constants::pi * (offset + n_ * q));
364  contactPoint.X = start_.X + q * l_ - thickness_;
365  getOrientation().rotate(contactPoint);
366  contactPoint += getPosition();
367  vtk.points.push_back(contactPoint);
368  }
369  for (unsigned ir = 0; ir < nr; ir++)
370  {
371  double q = (double) iz / nz;
372  double r = (double) (nr - 1 - ir) / (nr - 1) * maxR_;
373  contactPoint.Y = start_.Y - r * cos(2.0 * constants::pi * (offset + n_ * q));
374  contactPoint.Z = start_.Z - r * sin(2.0 * constants::pi * (offset + n_ * q));
375  contactPoint.X = start_.X + q * l_ + thickness_;
376  getOrientation().rotate(contactPoint);
377  contactPoint += getPosition();
378  vtk.points.push_back(contactPoint);
379  }
380  }
381  }
382 
383  for (unsigned iz = 0; iz < (screwType_==ScrewType::doubleHelix?2:1)*nz-1; iz++)
384  {
385  //skip step that would connect the two screw parts
386  if (iz==nz-1) ++iz;
387  std::vector<double> cell;
388  cell.reserve(2 * nr);
389  for (unsigned ir = 0; ir < 2*nr; ir++)
390  {
391  cell.push_back(nPoints + ir + 2*iz * nr);
392  cell.push_back(nPoints + ir + 2*(iz + 1) * nr);
393  }
394  vtk.triangleStrips.push_back(cell);
395  }
396 }
const unsigned nz
Definition: ConstraintElementsUnitTest.cpp:32

References cos(), doubleHelix, boost::multiprecision::fabs(), BaseInteractable::getOrientation(), BaseInteractable::getPosition(), l_, maxR_, n_, Mesh_Parameters::nz, offset_, constants::pi, VTKContainer::points, Eigen::numext::q, UniformPSDSelfTest::r, Quaternion::rotate(), screwType_, sin(), start_, thickness_, VTKContainer::triangleStrips, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by writeVTK().

Member Data Documentation

◆ l_

Mdouble Screw::l_
private

The length of the Screw.

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

◆ maxR_

Mdouble Screw::maxR_
private

The outer radius of the Screw.

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

◆ n_

Mdouble Screw::n_
private

The number of revelations.

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

◆ offset_

Mdouble Screw::offset_
private

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

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

◆ omega_

Mdouble Screw::omega_
private

Rotation speed in rev/s.

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

◆ screwType_

ScrewType Screw::screwType_
private

Single or double helix screw.

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

◆ start_

Vec3D Screw::start_
private

The centre of the lower end of the screw.

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

◆ thickness_

Mdouble Screw::thickness_
private

The thickness of the Screw.

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


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