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)
	Constructor in which all parameters of the screw are set. More...

	~Screw ()
	Default destructor. More...

Screw *	copy () 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...

void	move_time (Mdouble dt)
	Rotate the Screw for a period dt, so that the offset_ changes with omega_*dt. More...

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 (std::string filename) const

void	getTriangulation (std::vector< Vec3D > &vertex, std::vector< std::array< unsigned, 3 >> &face, unsigned nr=5, unsigned nz=15) const

Public Member Functions inherited from BaseWall
	BaseWall ()
	Default constructor. It makes an empty BaseWall. More...

	BaseWall (const BaseWall &w)
	Copy constructor. More...

virtual	~BaseWall ()
	Default destructor. More...

virtual MERCURY_DEPRECATED void	clear ()
	A function that removes all data from this BaseWall, so sets handler_ to nullptr. More...

virtual void	setHandler (WallHandler *handler)
	A function which sets the WallHandler for this BaseWall. More...

WallHandler *	getHandler () const
	A function which returns the WallHandler that handles this BaseWall. More...

void	setIndSpecies (unsigned int indSpecies)
	Define the species of this wall using the index of the species in the SpeciesHandler in this DPMBase. More...

void	setSpecies (const ParticleSpecies *species)
	Define the species of this wall. More...

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, const Vec3D normal, const Vec3D position) const

virtual void	writeVTK (VTKContainer &vtk) const

virtual std::vector < BaseInteraction * >	getInteractionWith (BaseParticle p, Mdouble timeStamp, InteractionHandler interactionHandler)

Public Member Functions inherited from BaseInteractable
	BaseInteractable ()
	Default BaseInteractable constructor, it simply creates an empty BaseInteractable. More...

	BaseInteractable (const BaseInteractable &p)
	Copy constructor. It copies the BaseInteractable and all objects it contains. More...

virtual	~BaseInteractable ()
	Destructor, it simply destructs the BaseInteractable and all the objects it contains. More...

unsigned int	getIndSpecies () const
	Returns the index of the Species of this BaseInteractable. More...

const ParticleSpecies *	getSpecies () const
	Returns a pointer to the species of this BaseInteractable. More...

void	setSpecies (const ParticleSpecies *species)
	Sets the species of this BaseInteractable. More...

const Vec3D &	getForce () const
	Returns the force on this BaseInteractable. More...

const Vec3D &	getTorque () 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...

const Vec3D &	getPosition () const
	Returns the position of this BaseInteractable. More...

const Vec3D &	getOrientation () const
	Returns the orientation of this BaseInteractable. More...

void	setPosition (const Vec3D &position)
	Sets the position of this BaseInteractable. More...

void	setOrientation (const Vec3D &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...

void	rotate (const Vec3D &rotate)
	Rotates this BaseInteractable. More...

const std::list < BaseInteraction * > &	getInteractions () const
	Returns a reference to the list of interactions in this BaseInteractable. 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 Vec3D &	getVelocity () const
	Returns the velocity of this interactable. More...

virtual const Vec3D &	getAngularVelocity () 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< Vec3D(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...

Public Member Functions inherited from BaseObject
	BaseObject ()
	Default constructor. More...

	BaseObject (const BaseObject &p)
	Copy constructor, copies all the objects BaseObject contains. More...

virtual	~BaseObject ()
	virtual destructor More...

virtual void	moveInHandler (const unsigned int index)
	Except that it is virtual, it does the same thing as setIndex() does. More...

void	setIndex (const unsigned int index)
	Allows one to assign an index to an object in the handler/container. More...

void	setId (const unsigned int 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...

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

Additional Inherited Members
Static Public Member Functions inherited from BaseWall
static void	addToVTK (const std::vector< Vec3D > &points, VTKContainer &vtk)

Public Attributes inherited from BaseWall
WallHandler *	handler_

Detailed Description

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

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

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

Definition at line 38 of file Screw.h.

Constructor & Destructor Documentation

Screw::Screw ( )

Default constructor: make a screw with default parameters.

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

Definition at line 37 of file Screw.cc.

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

Referenced by copy().

 {
     start_.setZero();
     l_ = 1.0;
     maxR_ = 1.0;
     n_ = 1.0;
     omega_ = 1.0;
     offset_ = 0.0;
     thickness_ = 0.0;
     logger(DEBUG, "Screw() constructor finished.");              
 }

Screw::Screw ( const Screw & other )

Copy constructor, copies another Screw.

Parameters

[in] other The Screw that has to be copied.

Definition at line 52 of file Screw.cc.

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

     : BaseWall(other)
 {
     start_ = other.start_;
     l_ = other.l_;
     maxR_ = other.maxR_;
     n_ = other.n_;
     omega_ = other.omega_;
     thickness_ = other.thickness_;
     offset_ = other.offset_;
     logger(DEBUG, "Screw(const Screw&) copy constructor finished.");
 }

Screw::Screw	(	Vec3D	start,
		Mdouble	l,
		Mdouble	r,
		Mdouble	n,
		Mdouble	omega,
		Mdouble	thickness
	)

Constructor in which all parameters of the screw are set.

Parameters

[in]	start	A Vec3D which denotes the centre of the lower end of the Screw.
[in]	l	The length of the Screw, must be positive.
[in]	r	The radius of the Screw, must be positive.
[in]	n	The number of revelations of the Screw, must be positive.
[in]	omega	The rotation speed of the Screw in rev/s.
[in]	thickness	The thickness of the Screw, must be non-negative.

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

Definition at line 75 of file Screw.cc.

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

 {
     start_ = start;
     l_ = l;
     maxR_ = r;
     n_ = n;
     omega_ = omega;
     thickness_ = thickness;
     offset_ = 0.0;
     logger(DEBUG, "Screw(Vec3D, Mdouble, Mdouble, Mdouble, Mdouble, Mdouble) constructor finished.");           
 }

Screw::~Screw ( )

Default destructor.

Definition at line 87 of file Screw.cc.

References DEBUG, and logger.

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

Member Function Documentation

Screw * Screw::copy ( ) const

finalvirtual

Copy this screw and return a pointer to the copy.

Returns: A pointer to a copy of this Screw.

Implements BaseWall.

Definition at line 95 of file Screw.cc.

References Screw().

 {
     return new Screw(*this);
 }

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.

Parameters

[in]	p	BaseParticle we want to calculate the distance and whether it collided of.
[out]	distance	The distance of the BaseParticle to this wall.
[out]	normal_return	If 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.

Implements BaseWall.

Definition at line 111 of file Screw.cc.

References A, mathsFunc::cos(), BaseInteractable::getPosition(), BaseParticle::getWallInteractionRadius(), l_, maxR_, n_, Vec3D::normalize(), offset_, constants::pi, R, mathsFunc::sign(), mathsFunc::sin(), constants::sqr_pi, mathsFunc::square(), start_, thickness_, Vec3D::X, Vec3D::Y, Z, and Vec3D::Z.

 {
     Mdouble RSquared = square(p.getPosition().X - start_.X) + square(p.getPosition().Y - start_.Y);
     Mdouble Z = p.getPosition().Z - start_.Z;
     //first do a simple check if particle is within the cylindrical hull of the screw
     if (RSquared > square(maxR_ + p.getWallInteractionRadius() + thickness_) 
         || Z > l_ + p.getWallInteractionRadius() + thickness_ 
         || Z < - p.getWallInteractionRadius() - thickness_)
     {
         //std::cout << "failed " << p.getPosition() << std::endl;
         return false;
     }
 
     //else:
     Mdouble R = sqrt(RSquared);
     Mdouble A = atan2(p.getPosition().Y - start_.Y, p.getPosition().X - start_.X);
     
     //after subtracting the start position and transforming the particle position from (XYZ) into (RAZ) 
     //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.
 
     //To find the contact point we have to minimize (with respect to r and q)
     //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
     //Using polar coordinates (i.e. x-x0=R*cos(A), y-y0=R*sin(A) and Z=z-z0)
     //distance^2=R^2+r^2-2*R*r*cos(A-2*pi*(offset+N*q))+(Z-q*L)^2
     
     //Assumption: d(distance)/dr=0 at minDistance (should there be also a q-derivative?)
     //Differentiate with respect to r and solve for zero:
     //0=2*r-2*R*cos(A-2*pi*(offset+N*q)
     //r=R*cos(A-2*pi*(offset+N*q))
     
     //Substitue back
     //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
     //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))+(Z-q*L)^2
     
     //So we have to minimize:
     //distance^2=R^2*sin^2(A-2*pi*(offset+N*q))^2 + (Z-q*L)^2 = f(q)
     //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])
     //f''(q)=(4*pi^2*N^2)*R^2*cos(2*A-4*pi*(offset+N*q)) - 2*L*L
     //For this we use the Euler algoritm
     
     Mdouble q; //Current guess
     Mdouble dd; //Derivative at current guess
     Mdouble ddd; //Second derivative at current guess
     Mdouble q0 = Z / l_; //assume closest point q0 is at same z-location as the particle
             
     //Set initial q to the closest position on the screw with the same angle as A
     //The initial guess will be in the minimum of the sin closest to q0
     //Minima of the sin are at
     //A-2*pi*(offset+N*q)=k*pi (k=integer)
     //q=A/(2*pi*N)-k/(2*N)-offset/N (k=integer)
     
     Mdouble k = round(A / constants::pi - 2.0 * (offset_ + n_ * q0)); // k: |A-a(q0)-k*pi|=min
     q = A / (2.0 * constants::pi * n_) - k / (2.0 * n_) - offset_ / n_; // q: a(q)=A
     
     //Now apply Newton's method
     do
     {
         Mdouble arg = 2.0 * A - 4.0 * constants::pi * (n_ * q + offset_);
         dd = -2.0 * constants::pi * n_ * RSquared * sin(arg) - 2.0 * l_ * (Z - q * l_);
         ddd = 8.0 * constants::sqr_pi * n_ * n_ * RSquared * cos(arg) + 2.0 * l_ * l_;
         q -= dd / ddd;
     } while (fabs(dd / ddd) > 1e-14);
     
     //Calculate r
     Mdouble r = R * cos(2.0 * constants::pi * (offset_ + n_ * q) - A);
     
     //Check if the location is actually on the screw:
     //First posibility is that the radius is too large:
     if (fabs(r) > maxR_) //Left boundary of the coil
     {
         r = mathsFunc::sign(r) * maxR_;
         unsigned int steps = 0;
         //This case reduces to the coil problem
         do
         {
             dd = -4.0 * R * r * constants::pi * n_ * sin(A - 2.0 * constants::pi * (n_ * q + offset_)) - 2.0 * l_ * (Z - q * l_);
             ddd = 8.0 * R * r * constants::sqr_pi * n_ * n_ * cos(A - 2.0 * constants::pi * (n_ * q + offset_)) + 2.0 * l_ * l_;
             q -= dd / ddd;
             steps++;
         } while (fabs(dd / ddd) > 1e-14);
     }
     //Second possibility is that it occured before the start of after the end
     if (q < 0)
     {
         q = 0;
         r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
         if (fabs(r) > maxR_)
         {
             r = mathsFunc::sign(r) * maxR_;
         }
     }
     else if (q > 1)
     {
         q = 1;
         r = R * cos(A - 2.0 * constants::pi * (offset_ + q * n_));
         if (fabs(r) > maxR_)
         {
             r = mathsFunc::sign(r) * maxR_;
         }
     }
     
     Mdouble distanceSquared = R * R * pow(sin(A - 2 * constants::pi * (offset_ + n_ * q)), 2) + pow(Z - q * l_, 2);
     //If distance is too large there is no contact
     if (distanceSquared >= square(p.getWallInteractionRadius() + thickness_))
     {
         return false;
     }
     
     Vec3D ContactPoint;
     distance = sqrt(distanceSquared) - thickness_;
     ContactPoint.X = start_.X + r * cos(2.0 * constants::pi * (offset_ + n_ * q));
     ContactPoint.Y = start_.Y + r * sin(2.0 * constants::pi * (offset_ + n_ * q));
     ContactPoint.Z = start_.Z + q * l_;
     normal_return = (ContactPoint - p.getPosition());
     normal_return.normalize();
     return true;
 }

std::string Screw::getName ( ) const

finalvirtual

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

Returns: The string "Screw".

Implements BaseObject.

Definition at line 288 of file Screw.cc.

Referenced by writeVTK().

 {
     return "Screw";
 }

void Screw::getTriangulation	(	std::vector< Vec3D > &	vertex,
		std::vector< std::array< unsigned, 3 >> &	face,
		unsigned	nr = `5`,
		unsigned	nz = `15`
	)		const

Definition at line 316 of file Screw.cc.

References mathsFunc::cos(), l_, maxR_, n_, offset_, constants::pi, mathsFunc::sin(), start_, Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by writeVTK().

 {
     //add vertices
     vertices.reserve(nr*nz);
     Vec3D vertex;
     for (Mdouble iz=0; iz<nz; iz++) {
         for (Mdouble ir=0; ir<nr; ir++) {
             Mdouble q = iz/(nz-1); //angular direction
             Mdouble r = ir/(nr-1)*maxR_;
             vertex.X = start_.X - r * cos(2.0 * constants::pi * (offset_ + n_ * q));
             vertex.Y = start_.Y - r * sin(2.0 * constants::pi * (offset_ + n_ * q));
             vertex.Z = start_.Z + q * l_;
             vertices.push_back(vertex);
         }
     }
     //reserve space for faces
     faces.reserve(2*(nr-1)*(nz-1));
     std::array<unsigned,3> face;
     for (unsigned iz=0; iz<nz-1; iz++) {
         for (unsigned ir=0; ir<nr-1; ir++) {
             //first face 0 1 nr
             faces.push_back({iz*nr+ir,iz*nr+(ir+1),(iz+1)*nr+ir});
             //second face 1 nr+1 nr
             faces.push_back({iz*nr+(ir+1),(iz+1)*nr+(ir+1),(iz+1)*nr+ir});
         }
     }
 }

void Screw::move_time ( Mdouble dt )

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

Parameters

[in] dt The time for which the Screw has to be turned.

Definition at line 232 of file Screw.cc.

References offset_, and omega_.

 {
     offset_ += omega_ * dt;
 }

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] is Input stream from which the Screw must be read.

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

Definition at line 259 of file Screw.cc.

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

 {
     std::string dummy;
     is >> dummy >> start_
             >> dummy >> l_
             >> dummy >> maxR_
             >> dummy >> n_
             >> dummy >> omega_
             >> dummy >> offset_;
 }

void Screw::read ( std::istream & is )

overridevirtual

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

Parameters

[in,out] is Input stream from which the Screw must be read.

Reimplemented from BaseWall.

Definition at line 240 of file Screw.cc.

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

 {
     BaseWall::read(is);
     std::string dummy;
     is >> dummy >> start_
             >> dummy >> l_
             >> dummy >> maxR_
             >> dummy >> n_
             >> dummy >> omega_
             >> dummy >> thickness_
             >> dummy >> offset_;
 }

void Screw::write ( std::ostream & os ) const

overridevirtual

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

Parameters

[in,out] os Output stream to which the Screw must be written.

Reimplemented from BaseWall.

Definition at line 273 of file Screw.cc.

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

 {
     BaseWall::write(os);
     os << " start " << start_
             << " Length " << l_
             << " Radius " << maxR_
             << " Revolutions " << n_
             << " Omega " << omega_
             << " Thickness " << thickness_
             << " Offset " << offset_;
 }

void Screw::writeVTK ( std::string filename ) const

Definition at line 293 of file Screw.cc.

References getName(), getTriangulation(), and helpers::writeToFile().

 {
     std::vector<Vec3D> vertices;
     std::vector<std::array<unsigned,3>> faces;
     getTriangulation (vertices, faces);
 
     std::stringstream file;
     file << "# vtk DataFile Version 2.0\n"
      << getName() << "\n"
      "ASCII\n"
      "DATASET UNSTRUCTURED_GRID\n"
      "POINTS " << vertices.size() << " double\n";
     for (const auto& vertex : vertices)
         file << vertex << '\n';
     file << "\nCELLS " << faces.size() << ' ' << 4*faces.size() << "\n";
     for (const auto& face : faces)
         file << "3 " << face[0] << ' ' << face[1] << ' ' << face[2] << '\n';
     file << "\nCELL_TYPES " << faces.size() << "\n";
     for (const auto& face : faces)
         file << "5\n";
     helpers::writeToFile(filename,file.str());
 }