Screw.h

Go to the documentation of this file.

//Copyright (c) 2013-2014, The MercuryDPM Developers Team. All rights reserved.
//For the list of developers, see <http://www.MercuryDPM.org/Team>.
//
//Redistribution and use in source and binary forms, with or without
//modification, are permitted provided that the following conditions are met:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name MercuryDPM nor the
//    names of its contributors may be used to endorse or promote products
//    derived from this software without specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
//ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
//WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
//DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
//LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
//ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
//(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
//SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#ifndef SCREW_H
#define SCREW_H

#include "BaseWall.h"
#include "ExtendedMath.h"





class Screw : public BaseWall
{
    public:
    
    Screw() : 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(Vec3D Start, double L, double R, double N, double omega,double 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              
    }
    
    virtual Screw* copy() const
    {
        return new Screw(*this);
    }   
 
    bool get_distance_and_normal(BaseParticle &P, Mdouble &distance, Vec3D &normal_return)
    {
        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;
    }

    void move_time(Mdouble dt)
    {
        offset+=omega*dt;
    }
    
    void read(std::istream& is)
    {
        std::string dummy;
        is >> dummy >> Start >> dummy >> L >> dummy >> MaxR >> dummy >> N >> dummy >> omega >> dummy >> offset;
    }

    void print(std::ostream& os) const
    { 
        os << "Screw start " << Start << " Length "<<L<<" Radius "<<MaxR<<" Revolutions "<<N<<" Omega "<<omega<<" Thickness "<<thickness<<" Offset "<<offset;
    }   
    
    Vec3D get_Velocity() const {return Vec3D(0.0,0.0,0.0);}
    
    private:
        Vec3D Start;
        double L,MaxR,N,omega,offset,thickness;
};

#endif