This function defines a 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 ()

	Screw (Vec3D Start, double L, double R, double N, double omega, double thickness)

virtual Screw *	copy () const

bool	get_distance_and_normal (BaseParticle &P, Mdouble &distance, Vec3D &normal_return)

void	move_time (Mdouble dt)
	Allows the wall to be moved to a new position (also orthogonal to the normal), and setting the velocity. More...

void	read (std::istream &is)
	reads wall More...

void	print (std::ostream &os) const
	outputs wall More...

Vec3D	get_Velocity () const
	Todo{Implement this function correctly}. More...

Public Member Functions inherited from BaseWall
	BaseWall ()

	BaseWall (const BaseWall &p)

virtual	~BaseWall ()

virtual void	clear ()

virtual void	move (Mdouble position_ UNUSED)
	Allows the wall to be moved to a new position. More...

virtual void	move (Vec3D velocity_ UNUSED, Mdouble dt UNUSED)
	Allows the wall to be moved to a new position (also orthogonal to the normal), and setting the velocity. More...

virtual bool	get_distance_and_normal (BaseParticle &P UNUSED, Mdouble &distance UNUSED, Vec3D &normal_return UNUSED)=0

void	set_Velocity (Vec3D new_)
	access function for velocity More...

virtual void	read (std::istream &is UNUSED)=0

virtual void	print (std::ostream &os UNUSED) const =0

virtual void	moveInHandler (int newPos)

int	getIndSpecies () const

void	setIndSpecies (const int _new)

const Vec3D &	get_Force () const

const Vec3D &	get_Torque () const

void	set_Force (Vec3D _new)

void	set_Torque (Vec3D _new)

void	set_Index (int index)

void	set_Id (int id)

void	add_Force (Vec3D _new)

void	add_Torque (Vec3D _new)

Private Attributes
Vec3D	Start

double	L

double	MaxR

double	N

double	omega

double	offset

double	thickness

Additional Inherited Members
Public Attributes inherited from BaseWall
int	indSpecies

Protected Attributes inherited from BaseWall
Vec3D	velocity
	velocity of the wall (used to calculate the relative velocity in the force calculation) More...

Detailed Description

This function defines a 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

Definition at line 41 of file Screw.h.

Constructor & Destructor Documentation

Screw::Screw ( )

inline

Definition at line 45 of file Screw.h.

References L, MaxR, N, offset, omega, Vec3D::set_zero(), Start, and thickness.

Referenced by copy().

             : BaseWall()
     {
         Start.set_zero();
         L=1;
         MaxR=1;
         N=1;
         omega=1;
         offset=0;
         thickness=0;
         #ifdef CONSTUCTOR_OUTPUT
             std::cerr << "Screw() finished" << std::endl;
         #endif              
     }   

Screw::Screw	(	Vec3D	Start,
		double	L,
		double	R,
		double	N,
		double	omega,
		double	thickness
	)

inline

Definition at line 59 of file Screw.h.

References L, MaxR, N, offset, omega, R, Start, and thickness.

                                                                                     : BaseWall()
     {
         this->Start=Start;
         this->L=L;
         this->MaxR=R;
         this->N=N;
         this->omega=omega;
         this->thickness=thickness;
         this->offset=0.0;
         #ifdef CONSTUCTOR_OUTPUT
             std::cerr << "Screw() finished" << std::endl;
         #endif              
     }

Member Function Documentation

virtual Screw* Screw::copy ( ) const

inlinevirtual

Implements BaseWall.

Definition at line 73 of file Screw.h.

References Screw().

     {
         return new Screw(*this);
     }   

bool Screw::get_distance_and_normal	(	BaseParticle &	P,
		Mdouble &	distance,
		Vec3D &	normal_return
	)

inline

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 coordinated (i.e. x-x0=R*cos(alpha), y-y0=R*sin(alpha) and dz=z-z0) Distance^2=R^2+r^2-2*R*r*cos(alpha-2*Pi*(offset+N*q))+(dz-q*L)^2

Differentiation with respect to r and solve for zero: 0=2*r-2*R*cos(alpha-2*Pi*(offset+N*q) r=R*cos(alpha-2*Pi*(offset+N*q))

Substitue back Distance^2=R^2+R^2*cos^2(alpha-2*Pi*(offset+N*q))-2*R^2*cos(alpha-2*Pi*(offset+N*q))*cos(alpha-2*Pi*(offset+N*q))+(dz-q*L)^2 Distance^2=R^2*sin^2(alpha-2*Pi*(offset+N*q))+(dz-q*L)^2

So we have to minimize: Distance^2=R^2*sin^2(alpha-2*Pi*(offset+N*q))^2+(dz-q*L)^2 For this we use the Euler algoritm

The initial guess will be in the minimum of the sin closest to the minimum of the parabolic part Minima of the sin are at alpha-2*Pi*(offset+N*q)=k*Pi (k=integer) q=alpha/(2*Pi*N)-k/(2*N)-offset/N (k=integer)

Definition at line 78 of file Screw.h.

References BaseParticle::get_Position(), BaseParticle::get_WallInteractionRadius(), Vec3D::GetLength(), L, MaxR, N, offset, constants::pi, R, mathsFunc::sign(), constants::sqr_pi, Start, thickness, Vec3D::X, Vec3D::Y, and Vec3D::Z.

     {
         double Rsqr=pow(P.get_Position().X-Start.X,2)+pow(P.get_Position().Y-Start.Y,2);
         if(Rsqr>pow(MaxR+P.get_WallInteractionRadius()+thickness,2)||P.get_Position().Z>L+Start.Z+P.get_WallInteractionRadius()+thickness||P.get_Position().Z<Start.Z-P.get_WallInteractionRadius()-thickness)
         {
             //std::cout<<"Particle is out of first bound checking"<<std::endl;
             //std::cout<<"Rule 1: "<<Rsqr<<"<"<<pow(r+P.get_Radius()+thickness,2)<<std::endl;
             //std::cout<<"Rule 2: "<<Start.Z-P.get_Radius()-thickness<<"<"<<P.get_Position().Z<<"<"<<L+Start.Z+P.get_Radius()+thickness<<std::endl;
             return false;
         }
         double R=sqrt(Rsqr);
         double alpha=atan2(P.get_Position().Y-Start.Y,P.get_Position().X-Start.X);
         double dz=P.get_Position().Z-Start.Z;
         
         
         
         
         
         double q;       //Current guess
         double dd;      //Derivative at current guess
         double ddd;     //Second derivative at current guess
         double q0=dz/L; //Minimum of the parabolic part
         
         
         double k=round(alpha/constants::pi-2.0*(offset+N*q0));
         q=alpha/(2.0*constants::pi*N)-k/(2.0*N)-offset/N;
         
         //Now apply Newton's method
         do
         {
             dd =-2.0*Rsqr*constants::pi*N*sin(2.0*alpha-4.0*constants::pi*(N*q+offset))-2.0*L*(dz-q*L);
             ddd= 8.0*Rsqr*constants::sqr_pi*N*N*cos(2.0*alpha-4.0*constants::pi*(N*q+offset))+2.0*L*L;
             q-=dd/ddd;
         } while(fabs(dd/ddd)>1e-14);
         
         //Calculate r
         double r=R*cos(2.0*constants::pi*(offset+N*q)-alpha);
         
         //Check if the location is actually on the screw:
         //First posibility is that the radius is to large:
         if(fabs(r)>MaxR) //Left boundary of the coil
         {
             r=mathsFunc::sign(r)*MaxR;
             int steps=0;
             //This case reduces to the coil problem
             do
             {
                 dd =-4.0*R*r*constants::pi    *N  *sin(alpha-2.0*constants::pi*(N*q+offset))-2.0*L*(dz-q*L);
                 ddd= 8.0*R*r*constants::sqr_pi*N*N*cos(alpha-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(alpha-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(alpha-2.0*constants::pi*(offset+q*N));
             if(fabs(r)>MaxR)
             {
                 r=mathsFunc::sign(r)*MaxR;
             }
         }
                     
         distance=R*R*pow(sin(alpha-2*constants::pi*(offset+N*q)),2)+pow(dz-q*L,2);
         //If distance is to large there is no contact
         if (distance>=(P.get_WallInteractionRadius()+thickness)*(P.get_WallInteractionRadius()+thickness))
         {
             return false;
         }
 
         Vec3D ContactPoint;
         distance=sqrt(distance)-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.get_Position();
         normal_return/=normal_return.GetLength();
         return true;
     }

Vec3D Screw::get_Velocity ( ) const

inlinevirtual

Todo{Implement this function correctly}.

Implements BaseWall.

Definition at line 205 of file Screw.h.

205 {return Vec3D(0.0,0.0,0.0);}

Vec3D

Implementation of a 3D vector (by Vitaliy).

Definition: Vector.h:40

void Screw::move_time ( Mdouble dt )

inlinevirtual

Allows the wall to be moved to a new position (also orthogonal to the normal), and setting the velocity.

Reimplemented from BaseWall.

Definition at line 186 of file Screw.h.

References offset, and omega.

     {
         offset+=omega*dt;
     }

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

inline

outputs wall

Definition at line 199 of file Screw.h.

References L, MaxR, N, offset, omega, Start, and thickness.

     { 
         os << "Screw start " << Start << " Length "<<L<<" Radius "<<MaxR<<" Revolutions "<<N<<" Omega "<<omega<<" Thickness "<<thickness<<" Offset "<<offset;
     }   

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

inline

reads wall

Definition at line 192 of file Screw.h.

References L, MaxR, N, offset, omega, and Start.

     {
         std::string dummy;
         is >> dummy >> Start >> dummy >> L >> dummy >> MaxR >> dummy >> N >> dummy >> omega >> dummy >> offset;
     }