HorizontalScrew 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 <HorizontalScrew.h>

+ Inheritance diagram for HorizontalScrew:

Public Member Functions

 HorizontalScrew ()
 Default constructor: make a screw with default parameters. More...
 
 HorizontalScrew (const HorizontalScrew &other)
 Copy constructor, copies another HorizontalScrew. More...
 
 HorizontalScrew (Vec3D start, Mdouble l, Mdouble minR, Mdouble lowerR, Mdouble diffR, Mdouble n, Mdouble omega, Mdouble thickness, const ParticleSpecies *s)
 Constructor in which all parameters of the screw are set. More...
 
 ~HorizontalScrew ()
 Default destructor. More...
 
HorizontalScrewcopy () 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 HorizontalScrew 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...
 
Mdouble getLength () const
 Returns the length of the HorizontalScrew. More...
 
Vec3D getStart () const
 Returns the starting position of the HorizontalScrew. More...
 
void move_time (Mdouble dt)
 Rotate the HorizontalScrew for a period dt, so that the offset_ changes with omega_*dt. More...
 
void read (std::istream &is) override
 Reads a HorizontalScrew from an input stream, for example a restart file. More...
 
void write (std::ostream &os) const override
 Writes this HorizontalScrew to an output stream, for example a restart file. More...
 
std::string getName () const final
 Returns the name of the object, here the string "HorizontalScrew". More...
 
void writeVTK (VTKContainer &vtk) const override
 
void setBlades (const Mdouble bladeWidth, const Mdouble bladeLength, const std::vector< Mdouble > bladeMounts)
 
- 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...
 
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 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 HorizontalScrew. More...
 
Mdouble minR_
 The outer radius of the HorizontalScrew. More...
 
Mdouble lowerR_
 
Mdouble diffR_
 
Mdouble n_
 The number of revelations. More...
 
Mdouble omega_
 Rotation speed in rev/s. More...
 
Mdouble offset_
 The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation. More...
 
Mdouble thickness_
 The thickness of the HorizontalScrew. More...
 
Mdouble bladeWidth_
 The maximum radial width of a blade (in r). More...
 
Mdouble bladeLength_
 The length of a blade (in the q-coordinate, which is a linear mapping from start.Z<z<start.Z+l to 0<q<1). More...
 
std::vector< MdoublebladeMounts_
 The starting point of a blade (in the q-coordinate, which is a linear mapping from start.Z<z<start.Z+l to 0<q<1) 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 HorizontalScrew be made more clear? I don't understand them.

Constructor & Destructor Documentation

◆ HorizontalScrew() [1/3]

HorizontalScrew::HorizontalScrew ( )

Default constructor: make a screw with default parameters.

Make a HorizontalScrew 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.

17 {
18  start_.setZero();
19  l_ = 1.0;
20  minR_ = 0.0;
21  lowerR_ = 1.0;
22  diffR_ = 0.0;
23  n_ = 1.0;
24  omega_ = 1.0;
25  offset_ = 0.0;
26  thickness_ = 0.0;
27  bladeLength_=0;
28  bladeWidth_=0;
29  bladeMounts_={};
30  logger(DEBUG, "HorizontalScrew() constructor finished.");
31 }
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.
@ DEBUG
Mdouble bladeLength_
The length of a blade (in the q-coordinate, which is a linear mapping from start.Z<z<start....
Definition: HorizontalScrew.h:123
Mdouble diffR_
Definition: HorizontalScrew.h:99
Mdouble thickness_
The thickness of the HorizontalScrew.
Definition: HorizontalScrew.h:115
Vec3D start_
The centre of the lower end of the screw.
Definition: HorizontalScrew.h:89
Mdouble omega_
Rotation speed in rev/s.
Definition: HorizontalScrew.h:107
Mdouble offset_
The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation.
Definition: HorizontalScrew.h:111
Mdouble n_
The number of revelations.
Definition: HorizontalScrew.h:103
Mdouble minR_
The outer radius of the HorizontalScrew.
Definition: HorizontalScrew.h:97
Mdouble l_
The length of the HorizontalScrew.
Definition: HorizontalScrew.h:93
Mdouble lowerR_
Definition: HorizontalScrew.h:98
Mdouble bladeWidth_
The maximum radial width of a blade (in r).
Definition: HorizontalScrew.h:119
std::vector< Mdouble > bladeMounts_
The starting point of a blade (in the q-coordinate, which is a linear mapping from start....
Definition: HorizontalScrew.h:127
void setZero()
Sets all elements to zero.
Definition: Vector.cc:23

References bladeLength_, bladeMounts_, bladeWidth_, DEBUG, diffR_, l_, logger, lowerR_, minR_, n_, offset_, omega_, Vec3D::setZero(), start_, and thickness_.

Referenced by copy().

◆ HorizontalScrew() [2/3]

HorizontalScrew::HorizontalScrew ( const HorizontalScrew other)

Copy constructor, copies another HorizontalScrew.

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

References bladeLength_, bladeMounts_, bladeWidth_, DEBUG, diffR_, l_, logger, lowerR_, minR_, n_, offset_, omega_, start_, and thickness_.

◆ HorizontalScrew() [3/3]

HorizontalScrew::HorizontalScrew ( Vec3D  start,
Mdouble  l,
Mdouble  minR,
Mdouble  lowerR,
Mdouble  diffR,
Mdouble  n,
Mdouble  omega,
Mdouble  thickness,
const ParticleSpecies s 
)

Constructor in which all parameters of the screw are set.

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

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

65 {
66  start_ = start;
67  l_ = l;
68  minR_ = minR;
69  lowerR_ = lowerR;
70  diffR_ = diffR;
71  n_ = n;
72  omega_ = omega;
73  thickness_ = thickness;
74  offset_ = 0.0;
75  bladeLength_=0;
76  bladeWidth_=0;
77  bladeMounts_={};
78  setSpecies(s);
79  logger(DEBUG, "HorizontalScrew(Vec3D, Mdouble, Mdouble, Mdouble, Mdouble, Mdouble) constructor finished.");
80 }
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
Definition: BaseWall.cc:148
RealScalar s
Definition: level1_cplx_impl.h:130
Vector::Scalar omega(const Vector &t, const Vector &s, RealScalar angle)
Definition: IDRS.h:36
void start(const unsigned &i)
(Re-)start i-th timer
Definition: oomph_utilities.cc:243

References bladeLength_, bladeMounts_, bladeWidth_, DEBUG, diffR_, l_, logger, lowerR_, minR_, n, n_, offset_, Eigen::internal::omega(), omega_, s, BaseWall::setSpecies(), oomph::CumulativeTimings::start(), start_, and thickness_.

◆ ~HorizontalScrew()

HorizontalScrew::~HorizontalScrew ( )

Default destructor.

83 {
84  logger(DEBUG, "~HorizontalScrew() finished, destroyed the HorizontalScrew.");
85 }

References DEBUG, and logger.

Member Function Documentation

◆ copy()

HorizontalScrew * HorizontalScrew::copy ( ) const
finalvirtual

Copy this screw and return a pointer to the copy.

Returns
A pointer to a copy of this HorizontalScrew.

Implements BaseWall.

91 {
92  return new HorizontalScrew(*this);
93 }
HorizontalScrew()
Default constructor: make a screw with default parameters.
Definition: HorizontalScrew.cc:16

References HorizontalScrew().

◆ getDistanceAndNormal()

bool HorizontalScrew::getDistanceAndNormal ( const BaseParticle p,
Mdouble distance,
Vec3D normal_return 
) const
finalvirtual

Compute the distance from the HorizontalScrew 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.

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 HorizontalScrew. If there is a collision, this function also computes the distance between the BaseParticle and HorizontalScrew and the normal of the IntersectionOfWalls at the intersection point.

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

Implements BaseWall.

107 {
108  Mdouble RSquared = square(p.getPosition().X - start_.X) + square(p.getPosition().Y - start_.Y);
109  Mdouble Z = p.getPosition().Z - start_.Z;
110  //first do a simple check if particle is within the cylindrical hull of the screw
111  Mdouble maxR = std::max(minR_,lowerR_+diffR_*Z/l_)+bladeWidth_;//todo: fix
112  Mdouble interactionRadius = p.getWallInteractionRadius(this) + thickness_;
113  if (RSquared > square(maxR + interactionRadius)
114  || Z > l_ + interactionRadius
115  || Z < - interactionRadius)
116  {
117  return false;
118  }
119 
120  //else:
121  Mdouble R = sqrt(RSquared);
122  Mdouble A = atan2(p.getPosition().Y - start_.Y, p.getPosition().X - start_.X);
123 
124  //after subtracting the start position and transforming the particle position from (XYZ) into (RAZ)
125  //coordinates, we compute the distance to the wall at, located at (r,a,z)=(r,2*pi*(offset+N*q+k/2),q*L), 0<q<1.
126 
127  //To find the contact point we have to minimize (with respect to r and q)
128  //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
129  //Using polar coordinates (i.e. x-x0=R*cos(A), y-y0=R*sin(A) and Z=z-z0)
130  //distance^2=R^2+r^2-2*R*r*cos(A-2*pi*(offset+N*q))+(Z-q*L)^2
131 
132  //Assumption: d(distance)/dr=0 at minDistance (should there be also a q-derivative?)
133  //Differentiate with respect to r and solve for zero:
134  //0=2*r-2*R*cos(A-2*pi*(offset+N*q)
135  //r=R*cos(A-2*pi*(offset+N*q))
136 
137  //Substitue back
138  //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
139  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))+(Z-q*L)^2
140 
141  //So we have to minimize:
142  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))^2 + (Z-q*L)^2 = f(q)
143  //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])
144  //f''(q)=(4*pi^2*N^2)*R^2*cos(2*A-4*pi*(offset+N*q)) - 2*L*L
145  //For this we use the Euler algoritm
146 
147  Mdouble q; //Current guess
148  Mdouble dd; //Derivative at current guess
149  Mdouble ddd; //Second derivative at current guess
150  Mdouble q0 = Z / l_; //assume closest point q0 is at same z-location as the particle
151 
152  //Set initial q to the closest position on the screw with the same angle as A
153  //The initial guess will be in the minimum of the sin closest to q0
154  //Minima of the sin are at
155  //A-2*pi*(offset+N*q)=k*pi (k=integer)
156  //q=A/(2*pi*N)-k/(2*N)-offset/N (k=integer)
157 
158  Mdouble k = round(A / constants::pi - 2.0 * (offset_ + n_ * q0)); // k: |A-a(q0)-k*pi|=min
159  q = A / (2.0 * constants::pi * n_) - k / (2.0 * n_) - offset_ / n_; // q: a(q)=A
160 
161  //this makes the trioliet screw unique: only one turn
162  if (((int)k)%2==0)
163  return false;
164 
165  //Now apply Newton's method to find the rel height q of the contact point
166  do
167  {
168  Mdouble arg = 2.0 * A - 4.0 * constants::pi * (n_ * q + offset_);
169  dd = -2.0 * constants::pi * n_ * RSquared * sin(arg) - 2.0 * l_ * (Z - q * l_);
170  ddd = 8.0 * constants::sqr_pi * n_ * n_ * RSquared * cos(arg) + 2.0 * l_ * l_;
171  q -= dd / ddd;
172  } while (fabs(dd / ddd) > 1e-14);
173 
174  //Calculate r
175  Mdouble r = R * cos(2.0 * constants::pi * (offset_ + n_ * q) - A);
176  maxR = std::max(minR_,lowerR_+diffR_*q);//todo: fix
177  for (auto b : bladeMounts_)
178  {
179  Mdouble dq = q-b;
180  if (dq>0 && dq<bladeLength_)
181  {
182  maxR += bladeWidth_*(dq/bladeLength_);
183  }
184  }
185 
186  //Check if the location is actually on the screw:
187  //First possibility is that the radius is too large:
188  if (fabs(r) > maxR) //Left boundary of the coil
189  {
190  //either the contact point is too far from the screw ...
191  if (fabs(r) > maxR+thickness_)
192  return false;
193  //... or we have to compute the contact around the edge
194  r = mathsFunc::sign(r) * maxR;
195  //This case reduces to the coil problem
196  do
197  {
198  dd = -4.0 * R * r * constants::pi * n_ * sin(A - 2.0 * constants::pi * (n_ * q + offset_)) - 2.0 * l_ * (Z - q * l_);
199  ddd = 8.0 * R * r * constants::sqr_pi * n_ * n_ * cos(A - 2.0 * constants::pi * (n_ * q + offset_)) + 2.0 * l_ * l_;
200  q -= dd / ddd;
201  } while (fabs(dd / ddd) > 1e-14);
202  }
203  //Second possibility is that it occurred before the start of after the end
204  if (q < 0)
205  {
206  q = 0;
207  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
208  if (fabs(r) > maxR)
209  {
210  r = mathsFunc::sign(r) * maxR;
211  }
212  }
213  else if (q > 1)
214  {
215  q = 1;
216  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
217  if (fabs(r) > maxR)
218  {
219  r = mathsFunc::sign(r) * maxR;
220  }
221  }
222 
223  Mdouble distanceSquared = square(R*sin(A - 2 * constants::pi * (offset_ + n_ * q))) + square(Z - q * l_);
224  //If distance is too large there is no contact
225  if (distanceSquared >= square(interactionRadius))
226  {
227  return false;
228  }
229 
230  Vec3D ContactPoint;
231  distance = sqrt(distanceSquared) - thickness_;
232  ContactPoint.X = start_.X + r * cos(2.0 * constants::pi * (offset_ + n_ * q));
233  ContactPoint.Y = start_.Y + r * sin(2.0 * constants::pi * (offset_ + n_ * q));
234  ContactPoint.Z = start_.Z + q * l_;
235  normal_return = (ContactPoint - p.getPosition());
236  normal_return.normalise();
237  return true;
238 }
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.)
@ R
Definition: StatisticsVector.h:21
float * p
Definition: Tutorial_Map_using.cpp:9
Scalar * b
Definition: benchVecAdd.cpp:17
The matrix class, also used for vectors and row-vectors.
Definition: Eigen/Eigen/src/Core/Matrix.h:186
Definition: Kernel/Math/Vector.h:30
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 max(a, b)
Definition: datatypes.h:23
#define Z
Definition: icosphere.cpp:21
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)
r
Definition: UniformPSDSelfTest.py:20
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:117
const Mdouble sqr_pi
Definition: ExtendedMath.h:25
const Mdouble pi
Definition: ExtendedMath.h:23
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(), b, bladeLength_, bladeMounts_, bladeWidth_, cos(), diffR_, e(), boost::multiprecision::fabs(), k, l_, lowerR_, max, minR_, n_, Vec3D::normalise(), offset_, p, constants::pi, Eigen::numext::q, UniformPSDSelfTest::r, R, Eigen::bfloat16_impl::round(), mathsFunc::sign(), sin(), constants::sqr_pi, sqrt(), Eigen::square(), start_, thickness_, Vec3D::X, Vec3D::Y, Z, and Vec3D::Z.

◆ getLength()

Mdouble HorizontalScrew::getLength ( ) const

Returns the length of the HorizontalScrew.

241 {
242  return l_;
243 }

References l_.

◆ getName()

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

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

Returns
The string "HorizontalScrew".

Implements BaseObject.

311 {
312  return "HorizontalScrew";
313 }

◆ getStart()

Vec3D HorizontalScrew::getStart ( ) const

Returns the starting position of the HorizontalScrew.

246 {
247  return start_;
248 }

References start_.

◆ move_time()

void HorizontalScrew::move_time ( Mdouble  dt)

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

Parameters
[in]dtThe time for which the HorizontalScrew has to be turned.
254 {
255  offset_ += omega_ * dt;
256 }

References offset_, and omega_.

Referenced by HorizontalMixer::actionsAfterTimeStep(), and HorizontalMixer::setScrewWalls().

◆ read()

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

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

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

Reimplemented from BaseWall.

262 {
263  BaseWall::read(is);
264  std::string dummy;
265  unsigned n=0;
266  is >> dummy >> start_
267  >> dummy >> l_
268  >> dummy >> minR_
269  >> dummy >> lowerR_
270  >> dummy >> diffR_
271  >> dummy >> n_
272  >> dummy >> omega_
273  >> dummy >> thickness_
274  >> dummy >> offset_
275  >> dummy >> bladeLength_
276  >> dummy >> bladeWidth_
277  >> dummy >> n;
278  for (unsigned i=0; i<n; i++) {
279  Mdouble val;
280  is >> val;
281  bladeMounts_.push_back(val);
282  }
283 }
int i
Definition: BiCGSTAB_step_by_step.cpp:9
void read(std::istream &is) override
Function that reads a BaseWall from an input stream, usually a restart file.
Definition: BaseWall.cc:57
val
Definition: calibrate.py:119
std::string string(const unsigned &i)
Definition: oomph_definitions.cc:286

References bladeLength_, bladeMounts_, bladeWidth_, diffR_, i, l_, lowerR_, minR_, n, n_, offset_, omega_, BaseWall::read(), start_, oomph::Global_string_for_annotation::string(), thickness_, and calibrate::val.

◆ setBlades()

void HorizontalScrew::setBlades ( const Mdouble  bladeWidth,
const Mdouble  bladeLength,
const std::vector< Mdouble bladeMounts 
)
316 {
317  bladeWidth_ = bladeWidth;
318  bladeLength_ = bladeLength;
319  bladeMounts_ = bladeMounts;
320 }

References bladeLength_, bladeMounts_, and bladeWidth_.

◆ write()

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

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

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

Reimplemented from BaseWall.

289 {
290  BaseWall::write(os);
291  os << " start " << start_
292  << " length " << l_
293  << " minRadius " << minR_
294  << " lowerRadius " << lowerR_
295  << " diffRadius " << diffR_
296  << " revolutions " << n_
297  << " omega " << omega_
298  << " thickness " << thickness_
299  << " offset " << offset_
300  << " bladeLength " << bladeLength_
301  << " bladeWidth " << bladeWidth_
302  << " bladeMounts " << bladeMounts_.size();
303  for (auto n : bladeMounts_)
304  os << " " << n;
305 }
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 bladeLength_, bladeMounts_, bladeWidth_, diffR_, l_, lowerR_, minR_, n, n_, offset_, omega_, start_, thickness_, and BaseWall::write().

◆ writeVTK()

void HorizontalScrew::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.

323 {
324  unsigned nr = 10;
325  unsigned nz = 180;
326 
327  unsigned nPoints = vtk.points.size();
328  Vec3D contactPoint;
329  for (unsigned iz=0; iz<nz; iz++) {
330  double maxR = std::max(minR_,lowerR_+(double)iz/nz*diffR_);//todo: fix
331  for (auto b : bladeMounts_)
332  {
333  Mdouble dq = (double)iz/nz-b;
334  if (dq>0 && dq<bladeLength_)
335  {
336  maxR += bladeWidth_*(dq/bladeLength_);
337  }
338  }
339  for (unsigned ir=0; ir<nr; ir++) {
340  double q = (double)iz/nz;
341  double r = (double)ir/nr*maxR;
342  contactPoint.X = start_.X - r * cos(2.0 * constants::pi * (offset_ + n_ * q));
343  contactPoint.Y = start_.Y - r * sin(2.0 * constants::pi * (offset_ + n_ * q));
344  contactPoint.Z = start_.Z + q * l_;
345  vtk.points.push_back(contactPoint);
346  }
347  }
348 
349  for (unsigned iz=0; iz<nz-1; iz++) {
350  std::vector<double> cell;
351  cell.reserve(2*nr);
352  for (unsigned ir=0; ir<nr; ir++) {
353  cell.push_back(nPoints+ir+iz*nr);
354  cell.push_back(nPoints+ir+(iz+1)*nr);
355  }
356  vtk.triangleStrips.push_back(cell);
357  }
358 }
const unsigned nz
Definition: ConstraintElementsUnitTest.cpp:32
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:19
std::vector< Vec3D > points
Definition: BaseWall.h:18

References b, bladeLength_, bladeMounts_, bladeWidth_, cos(), diffR_, l_, lowerR_, max, minR_, n_, Mesh_Parameters::nz, offset_, constants::pi, VTKContainer::points, Eigen::numext::q, UniformPSDSelfTest::r, sin(), start_, VTKContainer::triangleStrips, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Member Data Documentation

◆ bladeLength_

Mdouble HorizontalScrew::bladeLength_
private

The length of a blade (in the q-coordinate, which is a linear mapping from start.Z<z<start.Z+l to 0<q<1).

Referenced by getDistanceAndNormal(), HorizontalScrew(), read(), setBlades(), write(), and writeVTK().

◆ bladeMounts_

std::vector<Mdouble> HorizontalScrew::bladeMounts_
private

The starting point of a blade (in the q-coordinate, which is a linear mapping from start.Z<z<start.Z+l to 0<q<1)

Referenced by getDistanceAndNormal(), HorizontalScrew(), read(), setBlades(), write(), and writeVTK().

◆ bladeWidth_

Mdouble HorizontalScrew::bladeWidth_
private

The maximum radial width of a blade (in r).

Referenced by getDistanceAndNormal(), HorizontalScrew(), read(), setBlades(), write(), and writeVTK().

◆ diffR_

Mdouble HorizontalScrew::diffR_
private

◆ l_

Mdouble HorizontalScrew::l_
private

The length of the HorizontalScrew.

Referenced by getDistanceAndNormal(), getLength(), HorizontalScrew(), read(), write(), and writeVTK().

◆ lowerR_

Mdouble HorizontalScrew::lowerR_
private

◆ minR_

Mdouble HorizontalScrew::minR_
private

The outer radius of the HorizontalScrew.

Referenced by getDistanceAndNormal(), HorizontalScrew(), read(), write(), and writeVTK().

◆ n_

Mdouble HorizontalScrew::n_
private

The number of revelations.

Referenced by getDistanceAndNormal(), HorizontalScrew(), read(), write(), and writeVTK().

◆ offset_

Mdouble HorizontalScrew::offset_
private

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

Referenced by getDistanceAndNormal(), HorizontalScrew(), move_time(), read(), write(), and writeVTK().

◆ omega_

Mdouble HorizontalScrew::omega_
private

Rotation speed in rev/s.

Referenced by HorizontalScrew(), move_time(), read(), and write().

◆ start_

Vec3D HorizontalScrew::start_
private

The centre of the lower end of the screw.

Referenced by getDistanceAndNormal(), getStart(), HorizontalScrew(), read(), write(), and writeVTK().

◆ thickness_

Mdouble HorizontalScrew::thickness_
private

The thickness of the HorizontalScrew.

Referenced by getDistanceAndNormal(), HorizontalScrew(), read(), and write().


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