ParticleHandler.cc

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.


#include <limits>
#include <Math/Helpers.h>
#include "ParticleHandler.h"
#include "DPMBase.h"
#include "SpeciesHandler.h"
#include "Species/ParticleSpecies.h"


ParticleHandler::ParticleHandler()
{
    largestParticle_ = nullptr;
    smallestParticle_ = nullptr;
    logger(DEBUG, "ParticleHandler::ParticleHandler() finished");
}

ParticleHandler::ParticleHandler(const ParticleHandler& PH)
    : BaseHandler<BaseParticle>() 
{
    clear();
    setDPMBase(PH.getDPMBase());
    largestParticle_ = nullptr;
    smallestParticle_ = nullptr;
    copyContentsFromOtherHandler(PH);
    if (objects_.size() != 0)
    {
        computeLargestParticle();
        computeSmallestParticle();
    }
    logger(DEBUG, "ParticleHandler::ParticleHandler(const ParticleHandler &PH) finished");
}

ParticleHandler ParticleHandler::operator =(const ParticleHandler& rhs)
{
    if (this != &rhs)
    {
        clear();
        largestParticle_ = nullptr;
        smallestParticle_ = nullptr;
        copyContentsFromOtherHandler(rhs);
        if (objects_.size() != 0)
        {
            computeLargestParticle();
            computeSmallestParticle();
        }
    }
    logger(DEBUG, "ParticleHandler::operator = (const ParticleHandler& rhs) finished");
    return *this;
}

ParticleHandler::~ParticleHandler()
{
    //First reset the pointers, such that they are not checked twice when removing particles
    largestParticle_ = nullptr;
    smallestParticle_ = nullptr;
    logger(DEBUG, "ParticleHandler::~ParticleHandler() finished");
}

void ParticleHandler::addObject(BaseParticle* P)
{
    //\todo I diasbled the MERCURY_DEPRECATED attribute until we actually fixed the deprecation issue. This will scare users away.
    // if (P->getSpecies() == nullptr)
    // {
    //     logger(WARN, "WARNING: The particle with ID % that is added in "
    //             "ParticleHandler::addObject does not have a species yet. "
    //             "Please make sure that you have "
    //             "set the species somewhere in the driver code.", P->getId());
    // }
    //Puts the particle in the Particle list
    BaseHandler<BaseParticle>::addObject(P);
    if (getDPMBase() != nullptr)
    {
        //This places the particle in this grid
        getDPMBase()->hGridInsertParticle(P);
        //This computes where the particle currently is in the grid
        getDPMBase()->hGridUpdateParticle(P);
    }
    //set the particleHandler pointer
    P->setHandler(this);
    //compute mass of the particle
    P->getSpecies()->computeMass(P) ;
    //Check if this particle has new extrema
    checkExtrema(P);
}

void ParticleHandler::removeObject(unsigned const int id)
{
#ifdef CONTACT_LIST_HGRID
    getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getObject(id));
#endif
    getDPMBase()->hGridRemoveParticle(getObject(id));
    BaseHandler<BaseParticle>::removeObject(id);
}

void ParticleHandler::removeLastObject()
{
#ifdef CONTACT_LIST_HGRID
    getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getLastObject());
#endif
    getDPMBase()->hGridRemoveParticle(getLastObject());
    BaseHandler<BaseParticle>::removeLastObject();
}

void ParticleHandler::computeSmallestParticle()
{
    if (getNumberOfObjects() == 0)
    {
        logger(WARN, "No particles, so cannot compute the smallest particle.");
        return;
    }
    Mdouble min = std::numeric_limits<Mdouble>::max();
    smallestParticle_ = nullptr;
    for (BaseParticle* const particle : objects_)
    {
        if (particle->getInteractionRadius() < min)
        {
            min = particle->getInteractionRadius();
            smallestParticle_ = particle;
        }
    }
}

void ParticleHandler::computeLargestParticle() 
{
    if (getNumberOfObjects() == 0)
    {
        logger(WARN, "No particles, so cannot compute the largest particle.");
        return;
    }
    Mdouble max = -std::numeric_limits<Mdouble>::max();
    largestParticle_ = nullptr;
    for (BaseParticle* const particle : objects_)
    {
        if (particle->getInteractionRadius() > max)
        {
            max = particle->getInteractionRadius();
            largestParticle_ = particle;
        }
    }
}

BaseParticle* ParticleHandler::getSmallestParticle() const
{
    if (smallestParticle_ == nullptr)
    {
        logger(WARN, "No particles to return in getSmallestParticle()");
    }
    return smallestParticle_;
}

BaseParticle* ParticleHandler::getLargestParticle() const
{
    if (largestParticle_ == nullptr)
    {
        logger(WARN, "No particles to set get_LargestParticle()");
    }
    return largestParticle_;
}

BaseParticle* ParticleHandler::getFastestParticle() const
{
    if (getNumberOfObjects() == 0)
    {
        logger(WARN, "No particles to set getFastestParticle()" );
        return nullptr;
    }
    BaseParticle* p = nullptr;
    Mdouble maxSpeed = -std::numeric_limits<Mdouble>::max();
    for (BaseParticle* const pLoop : objects_)
    {
        if ((pLoop->getVelocity().getLength()) > maxSpeed)
        {
            maxSpeed = pLoop->getVelocity().getLength();
            p = pLoop;
        }
    }
    return p;
}

BaseParticle* ParticleHandler::getLowestPositionComponentParticle(const int i) const
{
    if (getNumberOfObjects() == 0)
    {
        logger(WARN, "No getLowestPositionComponentParticle(const int i) since there are no particles.");
        return nullptr;
    }
    BaseParticle* p = nullptr;
    Mdouble min = std::numeric_limits<Mdouble>::max();
    for (BaseParticle* const pLoop : objects_)
    {
        if (pLoop->getPosition().getComponent(i) < min)
        {
            min = pLoop->getPosition().getComponent(i);
            p = pLoop;
        }
    }
    return p;
}

BaseParticle* ParticleHandler::getHighestPositionComponentParticle(const int i) const
{
    if (getNumberOfObjects() == 0)
    {
        logger(WARN, "No getHighestPositionComponentParticle(const int i) since there are no particles.");
        return nullptr;
    }
    BaseParticle* p = nullptr;
    Mdouble max = -std::numeric_limits<Mdouble>::max();
    for (BaseParticle* const pLoop : objects_)
    {
        if (pLoop->getPosition().getComponent(i) > max)
        {
            max = pLoop->getPosition().getComponent(i);
            p = pLoop;
        }
    }
    
    return p;
}

BaseParticle* ParticleHandler::getLowestVelocityComponentParticle(const int i) const
{
    if (getNumberOfObjects() == 0)
    {
        logger(WARN, "No getLowestVelocityComponentParticle(const int i) since there are no particles");
        return nullptr;
    }
    BaseParticle* p = nullptr;
    Mdouble min = std::numeric_limits<Mdouble>::max();
    for (BaseParticle* const pLoop : objects_)
    {
        if (pLoop->getVelocity().getComponent(i) < min)
        {
            min = pLoop->getVelocity().getComponent(i);
            p = pLoop;
        }
    }
    return p;
}

BaseParticle* ParticleHandler::getHighestVelocityComponentParticle(const int i) const
{
    if (!getNumberOfObjects())
    {
        logger(WARN, "No getHighestVelocityComponentParticle(const int i) since there are no particles");
        return nullptr;
    }
    BaseParticle* p = nullptr;
    Mdouble max = -std::numeric_limits<Mdouble>::max();
    for (BaseParticle* const pLoop : objects_)
    {
        if (pLoop->getVelocity().getComponent(i) > max)
        {
            max = pLoop->getVelocity().getComponent(i);
            p = pLoop;
        }
    }
    return p;
}

BaseParticle* ParticleHandler::getLightestParticle() const
{
    if (getNumberOfObjects() == 0)
    {
        logger(WARN, "No particles to set getLightestParticle()");
        return nullptr;
    }
    BaseParticle* p = nullptr;
    Mdouble minMass = std::numeric_limits<Mdouble>::max();
    for (BaseParticle* const pLoop : objects_)
    {
        if (pLoop->getMass() < minMass)
        {
            minMass = pLoop->getMass();
            p = pLoop;
        }
    }
    return p;
}

void ParticleHandler::clear()
{
    smallestParticle_ = nullptr;
    largestParticle_ = nullptr;
    BaseHandler<BaseParticle>::clear();
}

void ParticleHandler::readObject(std::istream& is)
{
    std::string type;
    is >> type;
    if (type == "BaseParticle")
    {
        BaseParticle baseParticle;
        is >> baseParticle;
        copyAndAddObject(baseParticle);
        getLastObject()->setId(baseParticle.getId()); //to ensure old id
    }
    else if (type == "BP")
    {
        BaseParticle baseParticle;
        baseParticle.oldRead(is);
        copyAndAddObject(baseParticle);
        getLastObject()->setId(baseParticle.getId()); //to ensure old id
    }
    else if (isdigit(type[0]))
    {
        readOldObject(type, is);
    }
    else
    {
        logger(ERROR, "Particle type % not understood in restart file. Particle has not been read.", type);
        return;
    }
}

void ParticleHandler::readOldObject(std::string type, std::istream& is)
{
    //read in next line
    std::stringstream line(std::stringstream::in | std::stringstream::out);
    helpers::getLineFromStringStream(is, line);
    logger(VERBOSE, line.str());
    //std::cout << line.str() << std::endl;

    BaseParticle particle;
    
    //Declare all properties of the particle
    unsigned int indSpecies;
    Mdouble radius, inverseMass, inverseInertia;
    Vec3D position, velocity, orientation, angularVelocity;
    
    //Read all values
    position.X = atof(type.c_str());

    line >> position.Y >> position.Z >> velocity >> radius >> orientation >> angularVelocity >> inverseMass >> inverseInertia >> indSpecies;
    
    //Put the values in the particle    
    particle.setIndSpecies(indSpecies);
    particle.setPosition(position);
    particle.setVelocity(velocity);
    particle.setRadius(radius);
    particle.setOrientation(orientation);
    particle.setAngularVelocity(angularVelocity);
    if (inverseMass == 0.0)
        particle.fixParticle();
    else
    {
        particle.setInertia(1./inverseInertia);
    }
    
    //Put the particle in the  handler
    copyAndAddObject(particle);
}

void ParticleHandler::write(std::ostream& os) const
{
    os << "Particles " << getNumberOfObjects() << std::endl;
    for (BaseParticle* it : *this)
        os << (*it) << std::endl;
}

void ParticleHandler::checkExtrema(BaseParticle* P)
{
    if (P == largestParticle_) 
    {
        //if the properties of the largest particle changes
        computeLargestParticle();
    } 
    else if (!largestParticle_ || P->getInteractionRadius() > largestParticle_->getInteractionRadius())
    {
        largestParticle_ = P;
    }
    
    if (P == smallestParticle_) 
    {
        //if the properties of the smallest particle changes
        computeSmallestParticle();
    } 
    else if (!smallestParticle_ || P->getInteractionRadius() < smallestParticle_->getInteractionRadius())
    {
        smallestParticle_ = P;
    }
}

void ParticleHandler::checkExtremaOnDelete(BaseParticle* P)
{
    if (P == largestParticle_)
    {
        computeLargestParticle();
    }
    if (P == smallestParticle_)
    {
        computeSmallestParticle();
    }
}

void ParticleHandler::computeAllMasses(unsigned int indSpecies)
{
     for (BaseParticle* particle : objects_)
     {
        if (particle->getIndSpecies() == indSpecies)
        {
            particle->getSpecies()->computeMass(particle);
        }
     }
}

void ParticleHandler::computeAllMasses()
{
    for (BaseParticle* particle : objects_)
    {
        particle->getSpecies()->computeMass(particle);
    }
}

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