PeriodicBoundary.cc

Go to the documentation of this file.

//Copyright (c) 2013-2017, 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 provid->d 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 provid->d 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 "PeriodicBoundary.h"
#include "ParticleHandler.h"
#include "Particles/BaseParticle.h"
#include "DPMBase.h"

//TODO below is for debugging purposes
#include <unistd.h>

PeriodicBoundary::PeriodicBoundary()
        : BaseBoundary()
{
    distanceLeft_ = std::numeric_limits<double>::quiet_NaN();
    distanceRight_ = std::numeric_limits<double>::quiet_NaN();
    scaleFactor_ = std::numeric_limits<double>::quiet_NaN();

#ifdef MERCURY_USE_MPI
    MPIContainer& communicator = MPIContainer::Instance();
    if (communicator.getNumberOfProcessors() > 1)
    {
        logger(WARN,"PeriodicBoundaries are currently not implemented in parallel MercuryDPM");
    }
#endif
    logger(DEBUG, "PeriodicBoundary::PeriodicBoundary() finished");
    logger(DEBUG, "PeriodicBoundary::PeriodicBoundary() finished");
}

PeriodicBoundary::~PeriodicBoundary()
{
    logger(DEBUG, "PeriodicBoundary::~PeriodicBoundary() finished");
}

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

void PeriodicBoundary::set(Vec3D normal, Mdouble distanceLeft, Mdouble distanceRight)
{
    // factor is used to set normal to unit length
    scaleFactor_ = 1. / std::sqrt(Vec3D::dot(normal, normal));
    normal_ = normal * scaleFactor_;
    distanceLeft_ = distanceLeft * scaleFactor_;
    distanceRight_ = distanceRight * scaleFactor_;
    logger.assert_always(distanceRight_ > distanceLeft_,
                         "PeriodicBoundary::set: left distance needs to be smaller than right distance");
    shift_ = normal_ * (distanceRight_ - distanceLeft_);
}

void PeriodicBoundary::set(Vec3D normal, Vec3D positionLeft, Vec3D positionRight)
{
    set(normal, Vec3D::dot(positionLeft,normal), Vec3D::dot(positionRight,normal));
}

void PeriodicBoundary::set(Vec3D normal, Mdouble distanceLeft, Mdouble distanceRight, Vec3D shiftDirection)
{
    // factor is used to set normal to unit length
    scaleFactor_ = 1. / std::sqrt(Vec3D::dot(normal, normal));
    normal_ = normal * scaleFactor_;
    distanceLeft_ = distanceLeft * scaleFactor_;
    distanceRight_ = distanceRight * scaleFactor_;
    // factor is used to set shift vector to correct length
    scaleFactor_ = (distanceRight_ - distanceLeft_) * Vec3D::dot(shiftDirection, normal_);
    shift_ = shiftDirection * scaleFactor_;
}

Vec3D PeriodicBoundary::getNormal() const
{
    return normal_;
}

Mdouble PeriodicBoundary::getDistanceLeft() const
{
    return distanceLeft_;
}

Mdouble PeriodicBoundary::getDistanceRight() const
{
    return distanceRight_;
}

Vec3D PeriodicBoundary::getShift() const
{
    return shift_;
}

void PeriodicBoundary::moveLeft(Mdouble distanceLeft)
{
    distanceLeft_ = distanceLeft * scaleFactor_;
    shift_ = normal_ * (distanceRight_ - distanceLeft_);
}

void PeriodicBoundary::moveRight(Mdouble distanceRight)
{
    distanceRight_ = distanceRight * scaleFactor_;
    shift_ = normal_ * (distanceRight_ - distanceLeft_);
}

Mdouble PeriodicBoundary::getDistance(const BaseParticle& p) const
{
    return getDistance(p.getPosition());
}

Mdouble PeriodicBoundary::getDistance(const Vec3D& position) const
{
    Mdouble distanceFromPlaneThroughOrigin = Vec3D::dot(position, normal_);
    return std::min(distanceFromPlaneThroughOrigin - distanceLeft_, distanceRight_ - distanceFromPlaneThroughOrigin);
}

void PeriodicBoundary::shiftPosition(BaseParticle* p) const
{
    if (isClosestToLeftBoundary(*p))
    {
        p->move(shift_);
    }
    else
    {
        p->move(-shift_);
    }
}

/*
 Vec3D PeriodicBoundary::getShiftedPosition(Vec3D &position)
 {
 if (left_wall)
 {
 return position + shift_;
 }
 else
 {
 return position - shift_;
 }
 }
 */

void PeriodicBoundary::shiftPositions(Vec3D& position1, Vec3D& position2) const
{
    if (isClosestToLeftBoundary(position1))
    {
        position1 += shift_;
        position2 += shift_;
    }
    else
    {
        position1 -= shift_;
        position2 -= shift_;
    }
}

/*
 * \details Returns TRUE if last particle checked is closest to the 'left'
 * wall, and FALSE if it is closest to the 'right' wall.
 * \return      true if it is closest to the left boundary, false otherwise
 */
bool PeriodicBoundary::isClosestToLeftBoundary(const BaseParticle& p) const
{
    return isClosestToLeftBoundary(p.getPosition());
}

/*
 * \details Returns TRUE if last position checked is closest to the 'left'
 * wall, and FALSE if it is closest to the 'right' wall.
 * \return      true if it is closest to the left boundary, false otherwise
 */
bool PeriodicBoundary::isClosestToLeftBoundary(const Vec3D& p) const
{
    const Mdouble distance = Vec3D::dot(p, normal_);
    return (distanceRight_ - distance > distance - distanceLeft_);
}

/*
 void PeriodicBoundary::getCloseTogether(Vec3D &P, Vec3D &Q)
 {
 Mdouble PQdotn = Vec3D::Dot(P - Q, normal_);
 Mdouble shift_norm2 = shift.GetLength2();
 //Check if P is so far from Q that a shift would move it closer
 if (mathsFunc::square(PQdotn) > .25 * shift_norm2)
 {
 //Now determine the direction of the shift
 if (PQdotn > 0.0)
 P -= shift;
 else
 P += shift;
 }
 }
 */

void PeriodicBoundary::read(std::istream& is)
{
    BaseBoundary::read(is);
    std::string dummy;
    is >> dummy >> normal_
       >> dummy >> scaleFactor_
       >> dummy >> distanceLeft_
       >> dummy >> distanceRight_
       >> dummy >> shift_;
}

void PeriodicBoundary::oldRead(std::istream& is)
{
    std::string dummy;
    is >> dummy >> normal_
       >> dummy >> scaleFactor_
       >> dummy >> distanceLeft_
       >> dummy >> distanceRight_
       >> dummy >> shift_;
}

void PeriodicBoundary::write(std::ostream& os) const
{
    BaseBoundary::write(os);
    os << " normal " << normal_
       << " scaleFactor " << scaleFactor_
       << " distanceLeft " << distanceLeft_
       << " distanceRight " << distanceRight_
       << " shift " << shift_;
}

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

/*
 Vec3D& PeriodicBoundary::getNormal()
 {
 return normal_;
 }
 */

void PeriodicBoundary::createPeriodicParticles(BaseParticle* p, ParticleHandler& pH)
{
    //note that getDistance sets closestToLeftBoundary_ to true or false depending on which side is closest
    if (getDistance(*p) < p->getInteractionRadius() + pH.getLargestParticle()->getInteractionRadius())
    {
        //Step 1: Copy the particle to new ghost particle.
        BaseParticle* pGhost = p->copy();

        //Step 2: Copy the interactions of the ghost particle.
        pGhost->copyInteractionsForPeriodicParticles(*p);

        //Step 3: Shift the ghost to the 'reflected' location.
        shiftPosition(pGhost);

        //Step 4: If Particle is double shifted, get correct original particle
        BaseParticle* from = p;
        while (from->getPeriodicFromParticle() != nullptr)
            from = from->getPeriodicFromParticle();
        pGhost->setPeriodicFromParticle(from);

        pH.addObject(pGhost);
    }
}

bool PeriodicBoundary::checkBoundaryAfterParticleMoved(BaseParticle* p, ParticleHandler& pH UNUSED)
{
    if (getDistance(*p) < 0)
    {
        shiftPosition(p);
    }
    return false;
}