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

Definition at line 38 of file HorizontalScrew.h.

Constructor & Destructor Documentation

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.

Definition at line 37 of file HorizontalScrew.cc.

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

Referenced by copy().

38 {
39  start_.setZero();
40  l_ = 1.0;
41  minR_ = 0.0;
42  lowerR_ = 1.0;
43  diffR_ = 0.0;
44  n_ = 1.0;
45  omega_ = 1.0;
46  offset_ = 0.0;
47  thickness_ = 0.0;
48  bladeLength_=0;
49  bladeWidth_=0;
50  bladeMounts_={};
51  logger(DEBUG, "HorizontalScrew() constructor finished.");
52 }
Mdouble bladeLength_
The length of a blade (in the q-coordinate, which is a linear mapping from start.Z
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
Mdouble bladeWidth_
The maximum radial width of a blade (in r).
void setZero()
Sets all elements to zero.
Definition: Vector.cc:43
Mdouble offset_
The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation...
Mdouble l_
The length of the HorizontalScrew.
Mdouble thickness_
The thickness of the HorizontalScrew.
std::vector< Mdouble > bladeMounts_
The starting point of a blade (in the q-coordinate, which is a linear mapping from start...
Vec3D start_
The centre of the lower end of the screw.
Mdouble omega_
Rotation speed in rev/s.
Mdouble n_
The number of revelations.
Mdouble minR_
The outer radius of the HorizontalScrew.
HorizontalScrew::HorizontalScrew ( const HorizontalScrew other)

Copy constructor, copies another HorizontalScrew.

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

Definition at line 57 of file HorizontalScrew.cc.

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

58  : BaseWall(other)
59 {
60  start_ = other.start_;
61  l_ = other.l_;
62  minR_ = other.minR_;
63  lowerR_ = other.lowerR_;
64  diffR_ = other.diffR_;
65  n_ = other.n_;
66  omega_ = other.omega_;
67  thickness_ = other.thickness_;
68  offset_ = other.offset_;
72  logger(DEBUG, "HorizontalScrew(const HorizontalScrew&) copy constructor finished.");
73 }
Mdouble bladeLength_
The length of a blade (in the q-coordinate, which is a linear mapping from start.Z
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
Mdouble bladeWidth_
The maximum radial width of a blade (in r).
Mdouble offset_
The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation...
Mdouble l_
The length of the HorizontalScrew.
Mdouble thickness_
The thickness of the HorizontalScrew.
BaseWall()
Default constructor.
Definition: BaseWall.cc:36
std::vector< Mdouble > bladeMounts_
The starting point of a blade (in the q-coordinate, which is a linear mapping from start...
Vec3D start_
The centre of the lower end of the screw.
Mdouble omega_
Rotation speed in rev/s.
Mdouble n_
The number of revelations.
Mdouble minR_
The outer radius of the HorizontalScrew.
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.

Definition at line 85 of file HorizontalScrew.cc.

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

86 {
87  start_ = start;
88  l_ = l;
89  minR_ = minR;
90  lowerR_ = lowerR;
91  diffR_ = diffR;
92  n_ = n;
93  omega_ = omega;
94  thickness_ = thickness;
95  offset_ = 0.0;
96  bladeLength_=0;
97  bladeWidth_=0;
98  bladeMounts_={};
99  setSpecies(s);
100  logger(DEBUG, "HorizontalScrew(Vec3D, Mdouble, Mdouble, Mdouble, Mdouble, Mdouble) constructor finished.");
101 }
Mdouble bladeLength_
The length of a blade (in the q-coordinate, which is a linear mapping from start.Z
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
Mdouble bladeWidth_
The maximum radial width of a blade (in r).
Mdouble offset_
The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation...
Mdouble l_
The length of the HorizontalScrew.
Mdouble thickness_
The thickness of the HorizontalScrew.
std::vector< Mdouble > bladeMounts_
The starting point of a blade (in the q-coordinate, which is a linear mapping from start...
Vec3D start_
The centre of the lower end of the screw.
Mdouble omega_
Rotation speed in rev/s.
Mdouble n_
The number of revelations.
Mdouble minR_
The outer radius of the HorizontalScrew.
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
Definition: BaseWall.cc:171
HorizontalScrew::~HorizontalScrew ( )

Default destructor.

Definition at line 103 of file HorizontalScrew.cc.

References DEBUG, and logger.

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

Member Function Documentation

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.

Definition at line 111 of file HorizontalScrew.cc.

References HorizontalScrew().

112 {
113  return new HorizontalScrew(*this);
114 }
HorizontalScrew()
Default constructor: make a screw with default parameters.
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.

Definition at line 127 of file HorizontalScrew.cc.

References A, bladeLength_, bladeMounts_, bladeWidth_, mathsFunc::cos(), diffR_, BaseInteractable::getPosition(), BaseParticle::getWallInteractionRadius(), l_, lowerR_, minR_, n_, Vec3D::normalise(), offset_, constants::pi, R, helpers::round(), mathsFunc::sign(), mathsFunc::sin(), constants::sqr_pi, mathsFunc::square(), start_, thickness_, Vec3D::X, Vec3D::Y, Z, and Vec3D::Z.

128 {
129  Mdouble RSquared = square(p.getPosition().X - start_.X) + square(p.getPosition().Y - start_.Y);
130  Mdouble Z = p.getPosition().Z - start_.Z;
131  //first do a simple check if particle is within the cylindrical hull of the screw
132  Mdouble maxR = std::max(minR_,lowerR_+diffR_*Z/l_)+bladeWidth_;//todo: fix
133  Mdouble interactionRadius = p.getWallInteractionRadius(this) + thickness_;
134  if (RSquared > square(maxR + interactionRadius)
135  || Z > l_ + interactionRadius
136  || Z < - interactionRadius)
137  {
138  return false;
139  }
140 
141  //else:
142  Mdouble R = sqrt(RSquared);
143  Mdouble A = atan2(p.getPosition().Y - start_.Y, p.getPosition().X - start_.X);
144 
145  //after subtracting the start position and transforming the particle position from (XYZ) into (RAZ)
146  //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.
147 
148  //To find the contact point we have to minimize (with respect to r and q)
149  //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
150  //Using polar coordinates (i.e. x-x0=R*cos(A), y-y0=R*sin(A) and Z=z-z0)
151  //distance^2=R^2+r^2-2*R*r*cos(A-2*pi*(offset+N*q))+(Z-q*L)^2
152 
153  //Assumption: d(distance)/dr=0 at minDistance (should there be also a q-derivative?)
154  //Differentiate with respect to r and solve for zero:
155  //0=2*r-2*R*cos(A-2*pi*(offset+N*q)
156  //r=R*cos(A-2*pi*(offset+N*q))
157 
158  //Substitue back
159  //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
160  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))+(Z-q*L)^2
161 
162  //So we have to minimize:
163  //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))^2 + (Z-q*L)^2 = f(q)
164  //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])
165  //f''(q)=(4*pi^2*N^2)*R^2*cos(2*A-4*pi*(offset+N*q)) - 2*L*L
166  //For this we use the Euler algoritm
167 
168  Mdouble q; //Current guess
169  Mdouble dd; //Derivative at current guess
170  Mdouble ddd; //Second derivative at current guess
171  Mdouble q0 = Z / l_; //assume closest point q0 is at same z-location as the particle
172 
173  //Set initial q to the closest position on the screw with the same angle as A
174  //The initial guess will be in the minimum of the sin closest to q0
175  //Minima of the sin are at
176  //A-2*pi*(offset+N*q)=k*pi (k=integer)
177  //q=A/(2*pi*N)-k/(2*N)-offset/N (k=integer)
178 
179  Mdouble k = round(A / constants::pi - 2.0 * (offset_ + n_ * q0)); // k: |A-a(q0)-k*pi|=min
180  q = A / (2.0 * constants::pi * n_) - k / (2.0 * n_) - offset_ / n_; // q: a(q)=A
181 
182  //this makes the trioliet screw unique: only one turn
183  if (((int)k)%2==0)
184  return false;
185 
186  //Now apply Newton's method to find the rel height q of the contact point
187  do
188  {
189  Mdouble arg = 2.0 * A - 4.0 * constants::pi * (n_ * q + offset_);
190  dd = -2.0 * constants::pi * n_ * RSquared * sin(arg) - 2.0 * l_ * (Z - q * l_);
191  ddd = 8.0 * constants::sqr_pi * n_ * n_ * RSquared * cos(arg) + 2.0 * l_ * l_;
192  q -= dd / ddd;
193  } while (fabs(dd / ddd) > 1e-14);
194 
195  //Calculate r
196  Mdouble r = R * cos(2.0 * constants::pi * (offset_ + n_ * q) - A);
197  maxR = std::max(minR_,lowerR_+diffR_*q);//todo: fix
198  for (auto b : bladeMounts_)
199  {
200  Mdouble dq = q-b;
201  if (dq>0 && dq<bladeLength_)
202  {
203  maxR += bladeWidth_*(dq/bladeLength_);
204  }
205  }
206 
207  //Check if the location is actually on the screw:
208  //First possibility is that the radius is too large:
209  if (fabs(r) > maxR) //Left boundary of the coil
210  {
211  //either the contact point is too far from the screw ...
212  if (fabs(r) > maxR+thickness_)
213  return false;
214  //... or we have to compute the contact around the edge
215  r = mathsFunc::sign(r) * maxR;
216  unsigned int steps = 0;
217  //This case reduces to the coil problem
218  do
219  {
220  dd = -4.0 * R * r * constants::pi * n_ * sin(A - 2.0 * constants::pi * (n_ * q + offset_)) - 2.0 * l_ * (Z - q * l_);
221  ddd = 8.0 * R * r * constants::sqr_pi * n_ * n_ * cos(A - 2.0 * constants::pi * (n_ * q + offset_)) + 2.0 * l_ * l_;
222  q -= dd / ddd;
223  steps++;
224  } while (fabs(dd / ddd) > 1e-14);
225  }
226  //Second possibility is that it occurred before the start of after the end
227  if (q < 0)
228  {
229  q = 0;
230  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
231  if (fabs(r) > maxR)
232  {
233  r = mathsFunc::sign(r) * maxR;
234  }
235  }
236  else if (q > 1)
237  {
238  q = 1;
239  r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
240  if (fabs(r) > maxR)
241  {
242  r = mathsFunc::sign(r) * maxR;
243  }
244  }
245 
246  Mdouble distanceSquared = square(R*sin(A - 2 * constants::pi * (offset_ + n_ * q))) + square(Z - q * l_);
247  //If distance is too large there is no contact
248  if (distanceSquared >= square(interactionRadius))
249  {
250  return false;
251  }
252 
253  Vec3D ContactPoint;
254  distance = sqrt(distanceSquared) - thickness_;
255  ContactPoint.X = start_.X + r * cos(2.0 * constants::pi * (offset_ + n_ * q));
256  ContactPoint.Y = start_.Y + r * sin(2.0 * constants::pi * (offset_ + n_ * q));
257  ContactPoint.Z = start_.Z + q * l_;
258  normal_return = (ContactPoint - p.getPosition());
259  normal_return.normalise();
260  return true;
261 }
const Vec3D & getPosition() const
Returns the position of this BaseInteractable.
Mdouble X
the vector components
Definition: Vector.h:65
Mdouble bladeLength_
The length of a blade (in the q-coordinate, which is a linear mapping from start.Z
double Mdouble
Definition: GeneralDefine.h:34
void normalise()
Makes this Vec3D unit length.
Definition: Vector.cc:123
Mdouble bladeWidth_
The maximum radial width of a blade (in r).
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 offset_
The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation...
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
Mdouble getWallInteractionRadius(const BaseWall *wall) const
returns the radius plus the interactionDistance
Definition: BaseParticle.h:386
Mdouble l_
The length of the HorizontalScrew.
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
const Mdouble pi
Definition: ExtendedMath.h:45
Mdouble thickness_
The thickness of the HorizontalScrew.
Mdouble Y
Definition: Vector.h:65
Mdouble round(const Mdouble value, unsigned precision)
Definition: Helpers.cc:598
std::vector< Mdouble > bladeMounts_
The starting point of a blade (in the q-coordinate, which is a linear mapping from start...
Vec3D start_
The centre of the lower end of the screw.
Mdouble n_
The number of revelations.
Definition: Vector.h:49
T square(const T val)
squares a number
Definition: ExtendedMath.h:106
Mdouble minR_
The outer radius of the HorizontalScrew.
Mdouble Z
Definition: Vector.h:65
const Mdouble sqr_pi
Definition: ExtendedMath.h:47
Mdouble HorizontalScrew::getLength ( ) const

Returns the length of the HorizontalScrew.

Definition at line 263 of file HorizontalScrew.cc.

References l_.

264 {
265  return l_;
266 }
Mdouble l_
The length of the HorizontalScrew.
std::string HorizontalScrew::getName ( ) const
finalvirtual

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

Returns
The string "HorizontalScrew".

Implements BaseObject.

Definition at line 333 of file HorizontalScrew.cc.

334 {
335  return "HorizontalScrew";
336 }
Vec3D HorizontalScrew::getStart ( ) const

Returns the starting position of the HorizontalScrew.

Definition at line 268 of file HorizontalScrew.cc.

References start_.

269 {
270  return start_;
271 }
Vec3D start_
The centre of the lower end of the screw.
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.

Definition at line 276 of file HorizontalScrew.cc.

References offset_, and omega_.

277 {
278  offset_ += omega_ * dt;
279 }
Mdouble offset_
The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation...
Mdouble omega_
Rotation speed in rev/s.
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.

Definition at line 284 of file HorizontalScrew.cc.

References bladeLength_, bladeMounts_, bladeWidth_, diffR_, constants::i, l_, lowerR_, minR_, n_, offset_, omega_, BaseWall::read(), start_, and thickness_.

285 {
286  BaseWall::read(is);
287  std::string dummy;
288  unsigned n=0;
289  is >> dummy >> start_
290  >> dummy >> l_
291  >> dummy >> minR_
292  >> dummy >> lowerR_
293  >> dummy >> diffR_
294  >> dummy >> n_
295  >> dummy >> omega_
296  >> dummy >> thickness_
297  >> dummy >> offset_
298  >> dummy >> bladeLength_
299  >> dummy >> bladeWidth_
300  >> dummy >> n;
301  for (unsigned i=0; i<n; i++) {
302  Mdouble val;
303  is >> val;
304  bladeMounts_.push_back(val);
305  }
306 }
Mdouble bladeLength_
The length of a blade (in the q-coordinate, which is a linear mapping from start.Z
double Mdouble
Definition: GeneralDefine.h:34
const std::complex< Mdouble > i
Definition: ExtendedMath.h:51
Mdouble bladeWidth_
The maximum radial width of a blade (in r).
Mdouble offset_
The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation...
Mdouble l_
The length of the HorizontalScrew.
Mdouble thickness_
The thickness of the HorizontalScrew.
void read(std::istream &is) override
Function that reads a BaseWall from an input stream, usually a restart file.
Definition: BaseWall.cc:80
std::vector< Mdouble > bladeMounts_
The starting point of a blade (in the q-coordinate, which is a linear mapping from start...
Vec3D start_
The centre of the lower end of the screw.
Mdouble omega_
Rotation speed in rev/s.
Mdouble n_
The number of revelations.
Mdouble minR_
The outer radius of the HorizontalScrew.
void HorizontalScrew::setBlades ( const Mdouble  bladeWidth,
const Mdouble  bladeLength,
const std::vector< Mdouble bladeMounts 
)

Definition at line 338 of file HorizontalScrew.cc.

References bladeLength_, bladeMounts_, and bladeWidth_.

339 {
340  bladeWidth_ = bladeWidth;
341  bladeLength_ = bladeLength;
342  bladeMounts_ = bladeMounts;
343 }
Mdouble bladeLength_
The length of a blade (in the q-coordinate, which is a linear mapping from start.Z
Mdouble bladeWidth_
The maximum radial width of a blade (in r).
std::vector< Mdouble > bladeMounts_
The starting point of a blade (in the q-coordinate, which is a linear mapping from start...
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.

Definition at line 311 of file HorizontalScrew.cc.

References bladeLength_, bladeMounts_, bladeWidth_, diffR_, l_, lowerR_, minR_, n_, offset_, omega_, start_, thickness_, and BaseWall::write().

312 {
313  BaseWall::write(os);
314  os << " start " << start_
315  << " length " << l_
316  << " minRadius " << minR_
317  << " lowerRadius " << lowerR_
318  << " diffRadius " << diffR_
319  << " revolutions " << n_
320  << " omega " << omega_
321  << " thickness " << thickness_
322  << " offset " << offset_
323  << " bladeLength " << bladeLength_
324  << " bladeWidth " << bladeWidth_
325  << " bladeMounts " << bladeMounts_.size();
326  for (auto n : bladeMounts_)
327  os << " " << n;
328 }
Mdouble bladeLength_
The length of a blade (in the q-coordinate, which is a linear mapping from start.Z
Mdouble bladeWidth_
The maximum radial width of a blade (in r).
void write(std::ostream &os) const override
Function that writes a BaseWall to an output stream, usually a restart file.
Definition: BaseWall.cc:102
Mdouble offset_
The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation...
Mdouble l_
The length of the HorizontalScrew.
Mdouble thickness_
The thickness of the HorizontalScrew.
std::vector< Mdouble > bladeMounts_
The starting point of a blade (in the q-coordinate, which is a linear mapping from start...
Vec3D start_
The centre of the lower end of the screw.
Mdouble omega_
Rotation speed in rev/s.
Mdouble n_
The number of revelations.
Mdouble minR_
The outer radius of the HorizontalScrew.
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.

Definition at line 345 of file HorizontalScrew.cc.

References bladeLength_, bladeMounts_, bladeWidth_, mathsFunc::cos(), diffR_, l_, lowerR_, minR_, n_, offset_, constants::pi, VTKContainer::points, mathsFunc::sin(), start_, VTKContainer::triangleStrips, Vec3D::X, Vec3D::Y, and Vec3D::Z.

346 {
347  unsigned nr = 10;
348  unsigned nz = 180;
349 
350  unsigned nPoints = vtk.points.size();
351  Vec3D contactPoint;
352  for (unsigned iz=0; iz<nz; iz++) {
353  double maxR = std::max(minR_,lowerR_+(double)iz/nz*diffR_);//todo: fix
354  for (auto b : bladeMounts_)
355  {
356  Mdouble dq = (double)iz/nz-b;
357  if (dq>0 && dq<bladeLength_)
358  {
359  maxR += bladeWidth_*(dq/bladeLength_);
360  }
361  }
362  for (unsigned ir=0; ir<nr; ir++) {
363  double q = (double)iz/nz;
364  double r = (double)ir/nr*maxR;
365  contactPoint.X = start_.X - r * cos(2.0 * constants::pi * (offset_ + n_ * q));
366  contactPoint.Y = start_.Y - r * sin(2.0 * constants::pi * (offset_ + n_ * q));
367  contactPoint.Z = start_.Z + q * l_;
368  vtk.points.push_back(contactPoint);
369  }
370  }
371 
372  unsigned nCells = vtk.triangleStrips.size();
373  //vtk.triangleStrips.reserve(nCells+(nz-1));
374  for (unsigned iz=0; iz<nz-1; iz++) {
375  std::vector<double> cell;
376  cell.reserve(2*nr);
377  for (unsigned ir=0; ir<nr; ir++) {
378  cell.push_back(nPoints+ir+iz*nr);
379  cell.push_back(nPoints+ir+(iz+1)*nr);
380  }
381  vtk.triangleStrips.push_back(cell);
382  }
383 }
Mdouble X
the vector components
Definition: Vector.h:65
Mdouble bladeLength_
The length of a blade (in the q-coordinate, which is a linear mapping from start.Z
double Mdouble
Definition: GeneralDefine.h:34
Mdouble bladeWidth_
The maximum radial width of a blade (in r).
Mdouble offset_
The angle that describes how much the HorizontalScrew has turned, going from 0 to 1 for a rotation...
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
Mdouble l_
The length of the HorizontalScrew.
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44
const Mdouble pi
Definition: ExtendedMath.h:45
std::vector< std::vector< double > > triangleStrips
Definition: BaseWall.h:39
Mdouble Y
Definition: Vector.h:65
std::vector< Vec3D > points
Definition: BaseWall.h:38
std::vector< Mdouble > bladeMounts_
The starting point of a blade (in the q-coordinate, which is a linear mapping from start...
Vec3D start_
The centre of the lower end of the screw.
Mdouble n_
The number of revelations.
Definition: Vector.h:49
Mdouble minR_
The outer radius of the HorizontalScrew.
Mdouble Z
Definition: Vector.h:65

Member Data Documentation

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

Definition at line 144 of file HorizontalScrew.h.

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

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)

Definition at line 148 of file HorizontalScrew.h.

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

Mdouble HorizontalScrew::bladeWidth_
private

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

Definition at line 140 of file HorizontalScrew.h.

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

Mdouble HorizontalScrew::diffR_
private

Definition at line 120 of file HorizontalScrew.h.

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

Mdouble HorizontalScrew::l_
private

The length of the HorizontalScrew.

Definition at line 114 of file HorizontalScrew.h.

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

Mdouble HorizontalScrew::lowerR_
private

Definition at line 119 of file HorizontalScrew.h.

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

Mdouble HorizontalScrew::minR_
private

The outer radius of the HorizontalScrew.

Definition at line 118 of file HorizontalScrew.h.

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

Mdouble HorizontalScrew::n_
private

The number of revelations.

Definition at line 124 of file HorizontalScrew.h.

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

Mdouble HorizontalScrew::offset_
private

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

Definition at line 132 of file HorizontalScrew.h.

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

Mdouble HorizontalScrew::omega_
private

Rotation speed in rev/s.

Definition at line 128 of file HorizontalScrew.h.

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

Vec3D HorizontalScrew::start_
private

The centre of the lower end of the screw.

Definition at line 110 of file HorizontalScrew.h.

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

Mdouble HorizontalScrew::thickness_
private

The thickness of the HorizontalScrew.

Definition at line 136 of file HorizontalScrew.h.

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


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