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 <Particles/ThermalParticle.h>
#include "ParticleHandler.h"
#include "DPMBase.h"
#include "SpeciesHandler.h"
#include "Species/ParticleSpecies.h"
#include "Particles/LiquidFilmParticle.h"
#include "Particles/BaseParticle.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)
{
    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 index)
{
#ifdef CONTACT_LIST_HGRID
    getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getObject(id));
#endif
    getDPMBase()->hGridRemoveParticle(getObject(index));
    BaseHandler<BaseParticle>::removeObject(index);
}

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;
}

Mdouble ParticleHandler::getMeanRadius() const {
    Mdouble sumRadius = 0;
    for (BaseParticle* const p : objects_) {
        sumRadius += p->getRadius();
    }
    return sumRadius/getNumberOfObjects();
}
  
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();
}

BaseParticle* ParticleHandler::getNewObject(const std::string& type)
{
    if (type == "BaseParticle"||type == "BP"||isdigit(type[0])) 
    {
        return new BaseParticle;
    } 
    else if (type == "LiquidFilmParticle") 
    {
        return new LiquidFilmParticle;
    } 
    else 
    {
        logger(ERROR, "Particle type % not understood in restart file. Particle has not been read.", type);
        return nullptr;
    }
}

void ParticleHandler::readObject(std::istream& is)
{
    std::string type;
    is >> type;
    logger.log(Log::DEBUG, "ParticleHandler::readObject(is): reading type %.", type);
    if (type == "BaseParticle")
    {
        BaseParticle baseParticle;
        is >> baseParticle;
        //set the species pointer correctly
        baseParticle.setSpecies(getDPMBase()->speciesHandler.getObject(baseParticle.getIndSpecies()));
        copyAndAddObject(baseParticle);
        setId(getLastObject(),baseParticle.getId()); 
        //getLastObject()->setId(baseParticle.getId()); //to ensure old id
    }
    else if (type == "LiquidFilmParticle")
    {
        LiquidFilmParticle particle;
        is >> particle;
        //set the species pointer correctly
        particle.setSpecies(getDPMBase()->speciesHandler.getObject(particle.getIndSpecies()));
        copyAndAddObject(particle);
        setId(getLastObject(),particle.getId());
        //getLastObject()->setId(particle.getId()); //to ensure old id
    }
    else if (type == "ThermalParticle")
    {
        ThermalParticle particle;
        is >> particle;
        //set the species pointer correctly
        particle.setSpecies(getDPMBase()->speciesHandler.getObject(particle.getIndSpecies()));
        copyAndAddObject(particle);
        setId(getLastObject(),particle.getId());
        //getLastObject()->setId(particle.getId()); //to ensure old id
    }
    else if (type == "BP")
    {
        BaseParticle baseParticle;
        baseParticle.setHandler(this);//to be able to set the species
        baseParticle.oldRead(is);
        //set the species pointer correctly
        baseParticle.setSpecies(getDPMBase()->speciesHandler.getObject(baseParticle.getIndSpecies()));
        copyAndAddObject(baseParticle);
        setId(getLastObject(),baseParticle.getId()); 
        //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: You need to add this particle to ParticleHandler::readObject.", 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.setSpecies(getDPMBase()->speciesHandler.getObject(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";
}

void ParticleHandler::writeVTK() const
{
    static unsigned fileCounter = 0;
    const std::string fileName = getDPMBase()->getName() + "Particle_" +
            std::to_string(fileCounter++) + ".vtu";
    logger(INFO, "writing % (t=%, N=%)", fileName, getDPMBase()->getTime(), this->getNumberOfObjects());

    std::fstream file;
    file.open(fileName.c_str(), std::ios::out);
    if (file.fail())
    {
        logger(WARN,  "Error in writeToFile: file % could not be opened", fileName);
    }
    //
    file << "<?xml version=\"1.0\"?>\n\n";
    file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
    file << "<UnstructuredGrid>\n";
    file << "<Piece NumberOfPoints=\"" << this->getNumberOfObjects() << "\" NumberOfCells=\"" << 0 << "\">\n";
    file << "<Points>\n";
    file << "  <DataArray type=\"Float32\" Name=\"Position\" NumberOfComponents=\"3\" format=\"ascii\">\n";
    for (const auto& p: *this)
    {
        file << '\t' << p->getPosition() << '\n';        
    }   
    file << "  </DataArray>\n";
    file << "</Points>\n";
    file << "<PointData  Vectors=\"vector\">\n";
    file << "  <DataArray type=\"Float32\" Name=\"Velocity\" NumberOfComponents=\"3\" format=\"ascii\">\n";
    // Add velocity
    for (const auto& p: *this)
    {
        file << '\t' << p->getVelocity() << '\n';        
    }  
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"Float32\" Name=\"Radius\" format=\"ascii\">\n";
    // Add radius
    for (const auto& p: *this)
    {
        file << '\t' << p->getRadius() << '\n';        
    }  
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"Float32\" Name=\"SpeciesType\" format=\"ascii\">\n";
    // Add species type
    for (const auto& p: *this)
    {
        file << '\t' << p->getIndSpecies() << '\n';        
    }  
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"Float32\" Name=\"Info\" format=\"ascii\">\n";
    // Add species type
    for (const auto& p: *this)
    {
        file << '\t' << getDPMBase()->getInfo(*p) << '\n';
    }
    file << "  </DataArray>\n";
    //check if this type of Interaction has extra fields
    if (getLastObject()!= nullptr)
    {
        for (unsigned i = 0; i<getLastObject()->getNumberOfFieldsVTK(); i++)
        {
            file << "  <DataArray type=\"" << getLastObject()->getTypeVTK(i) << "\" Name=\"" << getLastObject()->getNameVTK(i) << "\" format=\"ascii\">\n";
            // Add species type
            for (const auto &p: *this)
            {
                for (auto f : p->getFieldVTK(i))
                    file << '\t' << f << '\n';
            }
            file << "  </DataArray>\n";
        }
    }
    file << "</PointData>\n";
    file << "<Cells>\n";
    file << "  <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";  
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n";   
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"UInt8\"  Name=\"types\" format=\"ascii\">\n";    
    file << "  </DataArray>\n";
    file << "</Cells>\n";
    file << "</Piece>\n";
    file << "</UnstructuredGrid>\n";
    file << "</VTKFile>\n";
    file.close();
}