Coil.cc

Go to the documentation of this file.

//Copyright (c) 2013-2020, 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.

#include "Coil.h"
#include "InteractionHandler.h"
#include "Particles/BaseParticle.h"

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

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

Coil::Coil(Vec3D start, Mdouble l, Mdouble r, Mdouble n, Mdouble omega, Mdouble thickness)
{
    start_ = start;
    l_ = l;
    r_ = r;
    n_ = n;
    omega_ = omega;
    thickness_ = thickness;
    offset_ = 0.0;
    logger(DEBUG, "Coil(params) constructor with parameters finished.");
}

Coil::~Coil()
{
    logger(DEBUG, "~Coil() destructor finished.");
}

void
Coil::set(Vec3D start, Mdouble length, Mdouble radius, Mdouble numberOfRevelations, Mdouble omega, Mdouble thickness)
{
    start_ = start;
    l_ = length;
    r_ = radius;
    n_ = numberOfRevelations;
    omega_ = omega;
    thickness_ = thickness;
    offset_ = 0.0;
}

Coil* Coil::copy() const
{
    return new Coil(*this);
}

bool Coil::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normal_return) const
{
    const Mdouble interactionRadius = p.getWallInteractionRadius(this) + thickness_;
    Mdouble Rsqr = pow(p.getPosition().X - start_.X, 2) + pow(p.getPosition().Y - start_.Y, 2);
    if (Rsqr > pow(r_ + interactionRadius, 2) ||
        Rsqr < pow(r_ - interactionRadius, 2) ||
        p.getPosition().Z > l_ + start_.Z + interactionRadius ||
        p.getPosition().Z < start_.Z - interactionRadius)
    {
        //std::cout<<"Particle is out of first bound checking"<<std::endl;
        //std::cout<<"Rule 1: "<<pow(r-P.getRadius()-thickness_,2)<<"<"<<Rsqr<<"<"<<pow(r+P.getRadius()+thickness_,2)<<std::endl;
        //std::cout<<"Rule 2: "<<Start.Z-P.getRadius()-thickness_<<"<"<<P.getPosition().Z<<"<"<<L+Start.Z+P.getRadius()+thickness_<<std::endl;
        return false;
    }
    Mdouble R = sqrt(Rsqr);
    Mdouble alpha = atan2(p.getPosition().Y - start_.Y, p.getPosition().X - start_.X);
    Mdouble dz = p.getPosition().Z - start_.Z;
    
    //To find the contact point we have to minimize (with respect to 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
    
    //So we have to minimize:
    //Distance^2=R^2+r^2-2*R*r*cos(alpha-2*Pi*(offset+N*q))+(dz-q*L)^2
    //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 = dz / l_; //Minimum of the parabolic part
    
    //The initial guess will be in the maximum of the cos closest to the minimum of the parabolic part
    //Minima of the cos are at
    //alpha-2*Pi*(offset+N*q)=2*k*Pi (k=integer)
    //q=alpha/(2*Pi*N)-k/N-offset/N (k=integer)
    
    Mdouble k = round(alpha / 2.0 / constants::pi - (offset_ + n_ * q0));
    q = alpha / (2 * constants::pi * n_) - k / n_ - offset_ / n_;
    
    //Now apply Newton's method
    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;
    } while (fabs(dd / ddd) > 1e-14);
    
    //Check if the location is actually on the coil, otherwise a point collision with the end of the coil calculated
    if (q < 0) //Left boundary of the coil
    {
        q = 0;
    }
    else if (q > 1) //right boundary of the coil
    {
        q = 1;
    }
    
    Mdouble distanceSquared =
            R * R + r_ * r_ - 2 * R * r_ * cos(alpha - 2 * constants::pi * (offset_ + n_ * q)) + pow(dz - q * l_, 2);
    //If distance is too large there is no contact
    if (distanceSquared >= mathsFunc::square(p.getWallInteractionRadius(this) + thickness_))
    {
        //std::cout<<"Particle is out of second bound checking, distance^2="<<distance<<" max="<<(P.getRadius()+thickness_)*(P.getRadius()+thickness_)<<std::endl;
        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 /= normal_return.getLength();
    return true;
}

void Coil::move_time(Mdouble dt)
{
    offset_ += omega_ * dt;
}

void Coil::read(std::istream& is)
{
    BaseWall::read(is);
    std::string dummy;
    is >> dummy >> start_
       >> dummy >> l_
       >> dummy >> r_
       >> dummy >> n_
       >> dummy >> omega_
       >> dummy >> thickness_
       >> dummy >> offset_;
}

void Coil::oldRead(std::istream& is)
{
    std::string dummy;
    is >> dummy >> start_
       >> dummy >> l_
       >> dummy >> r_ >> dummy >> n_
       >> dummy >> omega_
       >> dummy >> offset_;
}

void Coil::write(std::ostream& os) const
{
    BaseWall::write(os);
    os << " Start " << start_
       << " Length " << l_
       << " Radius " << r_
       << " Revolutions " << n_
       << " Omega " << omega_
       << " Thickness " << thickness_
       << " Offset " << offset_;
}

std::string Coil::getName() const
{
    return "Coil";
}