ParticleHandler.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 <limits>
#include <Math/Helpers.h>
#include <Particles/SphericalParticle.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"
#include "Particles/SuperQuadricParticle.h"
#include "MpiContainer.h"
#include "MpiDataClass.h"
#include "BoundaryHandler.h"
#include "DomainHandler.h"
#include "Domain.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_.empty())
    {
        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_.empty())
        {
            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::addExistingObject(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>::addExistingObject(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::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());
    }
#ifdef MERCURY_USE_MPI
    bool insertParticle;
    //Check if the particle P should be added to the current domain
    if (NUMBER_OF_PROCESSORS == 1)
    {
        insertParticle = true;
    }
    else
    {
        insertParticle  = getDPMBase()->mpiInsertParticleCheck(P);
    }

    //Add the particle if it is in the mpi domain or if the domain is not yet defined
    if(insertParticle)
    {
#endif
        //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);
            
            // Broadcasts the existance of a new particle
            getDPMBase()->handleParticleAddition(P->getId(), 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);
        if (!P->isSphericalParticle())
        {
            getDPMBase()->setRotation(true);
        }
#ifdef MERCURY_USE_MPI
        P->setPeriodicComplexity(std::vector<int>(0));
    }
    else
    {
        logger.assert_debug(!P->isMPIParticle(),"Can't add mpi particle as it does not exist");
        logger.assert_debug(!P->isPeriodicGhostParticle(),"Can't add mpi particle as it does not exist");
        //Somehwere a really new particle has been added, so to keep the ID's globally unique, we also update
        //the unique value on all other processors
        BaseHandler<BaseParticle>::increaseId();
    }

    //Add the particle to the ghost particle lists
    getDPMBase()->insertGhostParticle(P);
    //Check if this new particle requires an update in the mpi grid (interactionDistance).
    getDPMBase()->updateGhostGrid(P);

    //Delete the particle that was supposed to be added (but was not)
    if(!insertParticle)
    {
        delete P;
    }
#endif
}

void ParticleHandler::addObject(int fromProcessor, BaseParticle* p)
{
#ifdef MERCURY_USE_MPI
    MPIContainer& communicator = MPIContainer::Instance();

    //The processor that contains the particle that needs to be copied needs to identify the target, and communicate this
    MPIParticle pInfo;
    if (communicator.getProcessorID() == fromProcessor)
    {
        pInfo.copyDataFromParticleToMPIParticle(p);
    }

    //Broadcast from processor i
    communicator.broadcast(&pInfo,MercuryMPIType::PARTICLE,fromProcessor);
    copyDataFromMPIParticleToParticle(&pInfo, p, this);

    //All processors now know have the same information and we can proceed with the collective add
    addObject(p);
#else
    logger(WARN, "Function addObject(int fromProcessor, BaseParticle* p) should not be used in serial code");
#endif
}

void ParticleHandler::addGhostObject(int fromProcessor, int toProcessor, BaseParticle* p)
{
#ifdef MERCURY_USE_MPI
    MPIContainer& communicator = MPIContainer::Instance();
    if (fromProcessor == toProcessor)
    {
        if (communicator.getProcessorID() == fromProcessor)
        {
            addGhostObject(p);
        }

        return;
    }

    //Create communication information object
    MPIParticle pInfo;
    int tag;
    if (communicator.getProcessorID() == fromProcessor)
    {
        pInfo.copyDataFromParticleToMPIParticle(p);
        tag = fromProcessor*MAX_PROC*10 + toProcessor*10 + MercuryMPITag::PARTICLE_DATA;
        communicator.send(&pInfo, MercuryMPIType::PARTICLE, 1, toProcessor, tag);
    }

    if (communicator.getProcessorID() == toProcessor)
    {
        tag = fromProcessor*MAX_PROC*10 + toProcessor*10 + MercuryMPITag::PARTICLE_DATA;
        communicator.receive(&pInfo, MercuryMPIType::PARTICLE, 1, fromProcessor, tag);
    }

    //Sync the communication
    communicator.sync();

    //Add the data to the new particle
    if (communicator.getProcessorID() == toProcessor)
    {
        copyDataFromMPIParticleToParticle(&pInfo, p, this);
    }


    //Only toProcessor adds the particle, quietly without any check with other processors
    //IMPORTANT NOTE: When adding  a periodic particle in parallel, this is performed just before
    //finding new particles. For that reason we dont have to add a ghostparticle to the mpi communication lists
    if (communicator.getProcessorID() == toProcessor)
    {
        addGhostObject(p);
    }
#else
    logger(WARN,
           "Function addGhostObject(int fromProcessor, int toProcessor, BaseParticle* p) should not be used in serial code");
#endif
}

void ParticleHandler::addGhostObject(BaseParticle* P)
{
#ifdef MERCURY_USE_MPI
    if (P->getSpecies() == nullptr)
    {
        logger(WARN, "[ParticleHandler::adGhostObject(BaseParticle*)] "
                "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());
    }

    MPIContainer& communicator = MPIContainer::Instance();
    //Puts the particle in the Particle list
    BaseHandler<BaseParticle>::addGhostObject(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);
        
        // broadcast particle addition
        getDPMBase()->handleParticleAddition(P->getId(), 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);
#else
    logger(INFO,
           "Function ParticleHandler::mpiAddObject(BaseParticle* P) should only be called when compiling with parallel build on");
#endif
}

void ParticleHandler::removeObject(const unsigned int index)
{
#ifdef MERCURY_USE_MPI
    MPIContainer& communicator = MPIContainer::Instance();
    if (communicator.getNumberOfProcessors() > 1 )
    {
        logger(WARN, "[ParticleHandler::removeObject(const unsigned int)] Using the function removeObject in parallel could lead to out-of-sync communication, Instead use deletion boundaries");
    }
#endif
#ifdef CONTACT_LIST_HGRID
    getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getObject(id));
#endif
    getDPMBase()->hGridRemoveParticle(getObject(index));
    // broadcast the particle removal
    getDPMBase()->handleParticleRemoval(getObject(index)->getId());
    BaseHandler<BaseParticle>::removeObject(index);
}

void ParticleHandler::removeGhostObject(const unsigned int index)
{
#ifdef MERCURY_USE_MPI
#ifdef CONTACT_LIST_HGRID
    getDPMBase()->getPossibleContactList().remove_ParticlePosibleContacts(getObject(id));
#endif
    getDPMBase()->hGridRemoveParticle(getObject(index));
    // broadcast the particle removal
    getDPMBase()->handleParticleRemoval(getObject(index)->getId());
    BaseHandler<BaseParticle>::removeObject(index);
#else
    logger(ERROR, "This function should only be used interally for mpi routines");
#endif
}

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

void ParticleHandler::computeSmallestParticle()
{
    if (getSize() == 0)
    {
        logger(DEBUG, "No particles, so cannot compute the smallest particle.");
        
        smallestParticle_ = nullptr;
        return;
    }
    Mdouble min = std::numeric_limits<Mdouble>::max();
    smallestParticle_ = nullptr;
    for (BaseParticle* const particle : objects_)
    {
        //if (!(particle->isMPIParticle() || particle->isPeriodicGhostParticle()))
        {
            if (particle->getMaxInteractionRadius() < min)
            {
                min = particle->getMaxInteractionRadius();
                smallestParticle_ = particle;
            }
        }
    }
}

void ParticleHandler::computeLargestParticle()
{
    if (getSize() == 0)
    {
        logger(DEBUG, "No particles, so cannot compute the largest particle.");
        largestParticle_ = nullptr;
        return;
    }
    Mdouble max = -std::numeric_limits<Mdouble>::max();
    largestParticle_ = nullptr;
    for (BaseParticle* const particle : objects_)
    {
        //if (!(particle->isMPIParticle() || particle->isPeriodicGhostParticle()))
        {
            if (particle->getMaxInteractionRadius() > max)
            {
                max = particle->getMaxInteractionRadius();
                largestParticle_ = particle;
            }
        }
    }
}

BaseParticle* ParticleHandler::getSmallestParticleLocal() const
{
    return smallestParticle_;
}

BaseParticle* ParticleHandler::getSmallestParticle() const
{
#ifdef MERCURY_USE_MPI
    logger(WARN,"getSmallestParticle should not be used in parallel; use getSmallestInteractionRadius or "
             "ParticleSpecies::getSmallestParticleMass instead");
    return nullptr;
#else
    return getSmallestParticleLocal();
#endif
}

BaseParticle* ParticleHandler::getLargestParticleLocal() const
{
    return largestParticle_;
}

BaseParticle* ParticleHandler::getLargestParticle() const
{
#ifdef MERCURY_USE_MPIO
    logger(WARN,"getLargestParticle() should not be used in parallel; use getLargestInteractionRadius instead");
    return nullptr;
#else
    return getLargestParticleLocal();
#endif
}

Mdouble ParticleHandler::getKineticEnergyLocal() const
{
    Mdouble ene = 0;
    for (auto p : *this)
    {
        if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
        {
            ene += p->getKineticEnergy();
        }
    }
    return ene;
}

Mdouble ParticleHandler::getKineticEnergy() const
{
#ifdef MERCURY_USE_MPI
    Mdouble kineticEnergyLocal = getKineticEnergyLocal();
    Mdouble kineticEnergyGlobal = 0.0;

    //sum up over all domains
    MPIContainer& communicator = MPIContainer::Instance();
    communicator.allReduce(kineticEnergyLocal, kineticEnergyGlobal, MPI_SUM);

    return kineticEnergyGlobal;
#else
    return getKineticEnergyLocal();
#endif
}

Mdouble ParticleHandler::getRotationalEnergyLocal() const
{
    Mdouble ene = 0;
    for (auto p : *this)
    {
        if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
        {
            ene += p->getRotationalEnergy();
        }
    }
    return ene;
}

Mdouble ParticleHandler::getRotationalEnergy() const
{
#ifdef MERCURY_USE_MPI
    Mdouble rotationalEnergyLocal = getRotationalEnergyLocal();
    Mdouble rotationalEnergyGlobal = 0.0;

    //sum up over all domains
    MPIContainer& communicator = MPIContainer::Instance();
    communicator.allReduce(rotationalEnergyLocal, rotationalEnergyGlobal, MPI_SUM);

    return rotationalEnergyGlobal;
#else
    return getRotationalEnergyLocal();
#endif
}

Mdouble ParticleHandler::getMassLocal() const
{
    Mdouble m = 0;
    for (auto p : *this)
        if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
            m += p->getMass();
    return m;
}

Mdouble ParticleHandler::getMass() const
{
#ifdef MERCURY_USE_MPI
    Mdouble massLocal = getMassLocal();
    Mdouble massGlobal = 0.0;

    //Sum up over all domains
    MPIContainer& communicator = MPIContainer::Instance();
    communicator.allReduce(massLocal, massGlobal, MPI_SUM);

    return massGlobal;
#else
    return getMassLocal();
#endif
}

Vec3D ParticleHandler::getMassTimesPositionLocal() const
{
    Vec3D com = {0, 0, 0};
    for (auto p : *this)
        if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
            com += p->getMass() * p->getPosition();
    return com;
}

Vec3D ParticleHandler::getMassTimesPosition() const
{
#ifdef MERCURY_USE_MPI
    Vec3D massTimesPositionLocal = getMassTimesPositionLocal();
    Vec3D massTimesPositionGlobal = {0.0, 0.0, 0.0};

    // Sum up over all domains
    MPIContainer& communicator = MPIContainer::Instance();
    communicator.allReduce(massTimesPositionLocal.X, massTimesPositionGlobal.X, MPI_SUM);
    communicator.allReduce(massTimesPositionLocal.Y, massTimesPositionGlobal.Y, MPI_SUM);
    communicator.allReduce(massTimesPositionLocal.Z, massTimesPositionGlobal.Z, MPI_SUM);

    return massTimesPositionGlobal;
#else
    return getMassTimesPositionLocal();
#endif
}

Vec3D ParticleHandler::getCentreOfMass() const
{
    Mdouble m = getMass();
    if (m == 0)
    {
        Vec3D nanvec = {constants::NaN, constants::NaN, constants::NaN};
        return nanvec;
    }
    else
        return getMassTimesPosition() / m;
}

Vec3D ParticleHandler::getMomentum() const
{
    Vec3D momentum = {0, 0, 0};
    for (auto p : *this)
        if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
            momentum += p->getMomentum();
    return getMPISum(momentum);
}

Vec3D ParticleHandler::getAngularMomentum() const
{
    Vec3D momentum = {0, 0, 0};
    for (auto p : *this)
        if (!(p->isFixed() || p->isMPIParticle() || p->isPeriodicGhostParticle()))
            momentum += p->getAngularMomentum();
    return getMPISum(momentum);
}

BaseParticle* ParticleHandler::getFastestParticleLocal() const
{
    if (getSize() == 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->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
        {
            if ((pLoop->getVelocity().getLength()) > maxSpeed)
            {
                maxSpeed = pLoop->getVelocity().getLength();
                p = pLoop;
            }
        }
    }
    return p;
}

BaseParticle* ParticleHandler::getFastestParticle() const
{
#ifdef MERCURY_USE_MPI
    logger(ERROR,"This function should not be used in parallel");
#endif
    return getFastestParticleLocal();
}

/*
 * \returns The smallest interaction radius on a local domain
 */
Mdouble ParticleHandler::getSmallestInteractionRadiusLocal() const
{
    if (!(getSmallestParticleLocal() == nullptr))
    {
        return getSmallestParticleLocal()->getMaxInteractionRadius();
    }
    else
    {
        return 0.0;
    }
}

/*
 * \returns The smallest interaction radius
 */
Mdouble ParticleHandler::getSmallestInteractionRadius() const
{
#ifdef MERCURY_USE_MPI
    //Compute the local value
    Mdouble smallestInteractionRadiusLocal = getSmallestInteractionRadiusLocal();
    Mdouble smallestInteractionRadiusGlobal = 0.0;

    //Obtain the global value
    MPIContainer& communicator = MPIContainer::Instance();
    communicator.allReduce(smallestInteractionRadiusLocal, smallestInteractionRadiusGlobal, MPI_MIN);

    return smallestInteractionRadiusGlobal;

#else
    return getSmallestInteractionRadiusLocal();
#endif
}

/*
 * \returns the largest interaction radius on a local domain
 */
Mdouble ParticleHandler::getLargestInteractionRadiusLocal() const
{
    if (!(getLargestParticleLocal() == nullptr))
    {
        return getLargestParticle()->getMaxInteractionRadius();
    }
    else
    {
        return 0.0;
    }
}

/*
 * \returns the largest interaction radius
 */
Mdouble ParticleHandler::getLargestInteractionRadius() const
{
#ifdef MERCURY_USE_MPI
    Mdouble largestInteractionRadiusLocal = getLargestInteractionRadiusLocal();
    Mdouble largestInteractionRadiusGlobal = 0.0;

    //Obtain the global value
    MPIContainer& communicator = MPIContainer::Instance();
    communicator.allReduce(largestInteractionRadiusLocal, largestInteractionRadiusGlobal, MPI_MAX);

    return largestInteractionRadiusGlobal;
#else
    return getLargestInteractionRadiusLocal();
#endif
}

Mdouble ParticleHandler::getSumRadiusLocal() const
{
    Mdouble sumRadius = 0;
    for (BaseParticle* const p : objects_)
    {
        if (!(p->isMPIParticle() || p->isPeriodicGhostParticle()))
        {
            sumRadius += p->getRadius();
        }
    }
    return sumRadius;
}

Mdouble ParticleHandler::getMeanRadius() const
{
#ifdef MERCURY_USE_MPI
    Mdouble sumRadiusLocal = getSumRadiusLocal();
    unsigned numberOfRealParticlesLocal = getNumberOfRealObjectsLocal();

    Mdouble sumRadiusGlobal = 0.0;
    unsigned numberOfRealParticlesGlobal = 0;

    //Sum up over all domains
    MPIContainer& communicator = MPIContainer::Instance();
    communicator.allReduce(sumRadiusLocal, sumRadiusGlobal, MPI_SUM);
    communicator.allReduce(numberOfRealParticlesLocal, numberOfRealParticlesGlobal, MPI_SUM);

    return sumRadiusGlobal/numberOfRealParticlesGlobal;
#else
    return getSumRadiusLocal() / getSize();
#endif
}

BaseParticle* ParticleHandler::getLowestPositionComponentParticleLocal(const int i) const
{
#ifdef MERCURY_USE_MPI
    logger(ERROR,"getLowestPositionComponentParticle() not implemented yet in parallel");
#endif
    if (getSize() == 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->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
        {
            if (pLoop->getPosition().getComponent(i) < min)
            {
                min = pLoop->getPosition().getComponent(i);
                p = pLoop;
            }
        }
    }
    return p;
}

BaseParticle* ParticleHandler::getLowestPositionComponentParticle(const int i) const
{
#ifdef MERCURY_USE_MPI
    logger(ERROR,"This function should not be used in parallel");
#endif
    return getLowestPositionComponentParticleLocal(i);
}

BaseParticle* ParticleHandler::getHighestPositionComponentParticleLocal(const int i) const
{
    if (getSize() == 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->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
        {
            if (pLoop->getPosition().getComponent(i) > max)
            {
                max = pLoop->getPosition().getComponent(i);
                p = pLoop;
            }
        }
    }
    
    return p;
}

BaseParticle* ParticleHandler::getHighestPositionComponentParticle(const int i) const
{
#ifdef MERCURY_USE_MPI
    logger(ERROR,"This function should not be used in parallel");
#endif
    return getHighestPositionComponentParticleLocal(i);
}

BaseParticle* ParticleHandler::getLowestVelocityComponentParticleLocal(const int i) const
{
    if (getSize() == 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->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
        {
            if (pLoop->getVelocity().getComponent(i) < min)
            {
                min = pLoop->getVelocity().getComponent(i);
                p = pLoop;
            }
        }
    }
    return p;
}

BaseParticle* ParticleHandler::getLowestVelocityComponentParticle(const int i) const
{
#ifdef MERCURY_USE_MPI
    logger(ERROR,"This function should not be used in parallel");
#endif
    return getLowestVelocityComponentParticleLocal(i);
}

BaseParticle* ParticleHandler::getHighestVelocityComponentParticleLocal(const int i) const
{
    if (!getSize())
    {
        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->isMPIParticle() || pLoop->isPeriodicGhostParticle()))
        {
            if (pLoop->getVelocity().getComponent(i) > max)
            {
                max = pLoop->getVelocity().getComponent(i);
                p = pLoop;
            }
        }
    }
    return p;
}

BaseParticle* ParticleHandler::getHighestVelocityComponentParticle(const int i) const
{
#ifdef MERCURY_USE_MPI
    logger(ERROR,"This function should not be used in parallel");
#endif
    return getHighestVelocityComponentParticleLocal(i);
}

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

Mdouble ParticleHandler::getHighestPositionX() const
{
    //Define the attribute function
    std::function<Mdouble(BaseParticle*)> particleAttribute = [](BaseParticle* p) { return p->getPosition().X; };
    
    //Obtain the MAX attribute
    Mdouble positionXGlobal = getParticleAttribute<Mdouble>(particleAttribute, AttributeType::MAX);
    
    return positionXGlobal;
}

unsigned int ParticleHandler::getNumberOfFixedParticles() const
{
    return NFixedParticles_;
}

unsigned int ParticleHandler::getNumberOfUnfixedParticles() const
{
    return getNumberOfObjects() - NFixedParticles_;
}

BaseParticle* ParticleHandler::createObject(const std::string& type)
{
    if (type == "BaseParticle") {
        //for backwards compatibility
        return new SphericalParticle;
    }
    else if (type == "SphericalParticle")
    {
        return new SphericalParticle;
    }
    else if (type == "LiquidFilmParticle")
    {
        return new LiquidFilmParticle;
    }
    else if (type == "SuperQuadricParticle")
    {
        return new SuperQuadricParticle;
    }
    else if (type == "ThermalParticle")
    {
        return new ThermalParticle;
    }
    else
    {
        logger(WARN, "Particle type % not understood in restart file. Particle will not be read.", type);
        return nullptr;
    }
}

BaseParticle* ParticleHandler::readAndCreateObject(std::istream& is)
{
    std::string type;
    BaseParticle* particle = nullptr;
    if (getStorageCapacity() > 0)
    {
        is >> type;
        logger(DEBUG, "ParticleHandler::readAndCreateObject(is): type %", type);
        particle = createObject(type);
        particle->setHandler(this);
        particle->read(is);
    }
    else //for insertion boundaries
    {
        is >> type;
        logger(DEBUG, "ParticleHandler::readAndCreateObject(is): type %", type);
        particle = createObject(type);
        particle->setHandler(this);
        particle->read(is);
    }
    particle->setSpecies(getDPMBase()->speciesHandler.getObject(particle->getIndSpecies()));
    return particle;
}

void ParticleHandler::readAndAddObject(std::istream& is)
{
    BaseParticle* o = readAndCreateObject(is);
    addExistingObject(o);
}

// * \param[in] type The first value of the position.
// * \param[in] is The input stream from which the information is read.
// * \details The old objects did not have their type in the beginning of the line.
// *          Instead, the first string of the file was the position in x-direction.
// *          Since we already read the first string of the file, we need to give
// *          it to this function and convert it to the position in x-direction.
// *          The rest of the stream is then read in the usual way.
// */
//void ParticleHandler::readAndCreateOldObject(std::istream& is, const std::string& type)
//{
//    //read in next line
//    std::stringstream line;
//    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, euler, angularVelocity;
//
//    //Read all values
//    position.X = atof(type.c_str());
//
//    line >> position.Y >> position.Z >> velocity >> radius >> euler >> 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);
//    Quaternion q; q.setEuler(euler);
//    particle.setOrientation(q);
//    particle.setAngularVelocity(angularVelocity);
//    if (inverseMass == 0.0)
//        particle.fixParticle();
//    else
//    {
//        particle.setInverseInertia(MatrixSymmetric3D(1,0,0,1,0,1)*inverseInertia);
//    }
//
//    //Put the particle in the  handler
//    copyAndAddObject(particle);
//}

void ParticleHandler::write(std::ostream& os) const
{
#ifdef MERCURY_USE_MPI
    os << "Particles " << getNumberOfRealObjectsLocal() << std::endl;
    for (BaseParticle* it : *this)
    {
        if (!it->isPeriodicGhostParticle() && !it->isMPIParticle())
        {
            os << (*it) << '\n';
        }
    }
#else
    os << "Particles " << getSize() << '\n';
    for (BaseParticle* it : *this)
    {
        os << (*it) << '\n';
    }
#endif
    //os.flush();
}

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

void ParticleHandler::addedFixedParticle()
{
    NFixedParticles_++;
}

void ParticleHandler::removedFixedParticle()
{
    NFixedParticles_--;
}

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

Mdouble ParticleHandler::getVolumeLocal() const
{
    Mdouble volume = 0;
    for (auto p : *this)
    {
        if (!(p->isFixed() || p->isPeriodicGhostParticle() || p->isMPIParticle()))
            volume += p->getVolume();
    }
    return volume;
}

Mdouble ParticleHandler::getVolume() const
{
#ifdef MERCURY_USE_MPI
    Mdouble volumeLocal = getVolumeLocal();
    Mdouble volumeGlobal = 0.0;

    // sum up over all domains
    MPIContainer& communicator = MPIContainer::Instance();
    communicator.allReduce(volumeLocal, volumeGlobal, MPI_SUM);

    return volumeGlobal;
#else
    return getVolumeLocal();
#endif
}

unsigned int ParticleHandler::getNumberOfRealObjectsLocal() const
{
#ifdef MERCURY_USE_MPI
    const MPIContainer& communicator = MPIContainer::Instance();
    if (communicator.getNumberOfProcessors() > 1)
    {
        unsigned int numberOfFakeParticles = 0;
        numberOfFakeParticles += getDPMBase()->domainHandler.getCurrentDomain()->getNumberOfTrueMPIParticles();
        numberOfFakeParticles += getDPMBase()->periodicBoundaryHandler.getNumberOfPeriodicGhostParticles();
        logger.assert_debug(numberOfFakeParticles <= getSize(), "More fake particles than getSize()");
        return (getSize() - numberOfFakeParticles);
    }
    else
    {
        return getSize();
    }
#else
    return getSize();
#endif
}

unsigned int ParticleHandler::getNumberOfRealObjects() const
{
#ifdef MERCURY_USE_MPI
    MPIContainer& communicator = MPIContainer::Instance();
    unsigned int numberOfRealParticles = getNumberOfRealObjectsLocal();

    // \todo MX: use an allreduce here
    //Combine the total number of Particles into one number on processor 0
    communicator.reduce(numberOfRealParticles, MPI_SUM);

    //Broadcast new number to all the processorsi
    communicator.broadcast(numberOfRealParticles);
    return numberOfRealParticles;
#else
    return getSize();
#endif
}

/*
 * \returns the size of the ParticleHandler. In parallel code it throws an error to avoid confusion with real and fake particles
 */
unsigned int ParticleHandler::getNumberOfObjects() const
{
//#ifdef MERCURY_USE_MPI
//    MPIContainer& communicator = MPIContainer::Instance();
//    if (communicator.getNumberOfProcessors() > 1)
//    {
//        logger(WARN,"When compiling with MPI please do not use getNumberOfObjects(). Instead use: getNumberOfRealObjectsLocal(), getNumberOfRealObjects() or getSize()");
//    }
//#endif
    return getSize();
}

unsigned int ParticleHandler::getNumberOfFixedObjectsLocal() const
{
    unsigned int numberOfFixedParticles = 0;
    for (BaseParticle* particle : *this)
    {
        if (particle->isFixed())
        {
            numberOfFixedParticles++;
        }
    }
    return numberOfFixedParticles;
}

unsigned int ParticleHandler::getNumberOfFixedObjects() const
{
#ifdef MERCURY_USE_MPI
    unsigned int numberOfFixedParticlesLocal = getNumberOfFixedObjectsLocal();
    unsigned int numberOfFixedParticles = 0;

    MPIContainer& communicator = MPIContainer::Instance();
    communicator.allReduce(numberOfFixedParticlesLocal, numberOfFixedParticles, MPI_SUM);
    return numberOfFixedParticles;
#else
    return getNumberOfFixedObjectsLocal();
#endif
}

void ParticleHandler::actionsAfterTimeStep()
{
    for (auto i: *this)
    {
        i->actionsAfterTimeStep();
    }
}

double ParticleHandler::getLiquidFilmVolume() const {
    double liquidVolume = 0;
    for (auto i : objects_) {
        auto j = dynamic_cast<LiquidFilmParticle*>(i);
        if (j and !j->isMPIParticle()) liquidVolume += j->getLiquidVolume();
    }
    return getMPISum(liquidVolume);
};