Domain.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 "Domain.h"
#include "DPMBase.h"
#include "Logger.h"
#include "MpiContainer.h"
#include "MpiDataClass.h"
#include "Particles/SphericalParticle.h"
#include "Walls/BaseWall.h"
#include "InteractionHandler.h"
#include "DomainHandler.h"
#include "Math/Vector.h"
#include <limits>
#include <utility>
#include <vector>
#include <set>
#include <Particles/LiquidFilmParticle.h>

Domain::Domain()
{
    constructor();
#ifdef DEBUG_CONSTRUCTOR
    std::cout<<"Domain::Domain() finished"<<std::endl;
#endif
}

Domain::Domain(std::vector<unsigned> globalMeshIndex) : globalMeshIndex_(std::move(globalMeshIndex))
{
    constructor();
#ifdef DEBUG_CONSTRUCTOR
    std::cout<<"Domain::Domain() finished"<<std::endl;
#endif
}

Domain::Domain(const Domain& b)
        : BaseObject(b)
{
    rank_ = b.rank_;
    domainHandler_ = b.domainHandler_;
    domainMin_ = b.domainMin_;
    domainMax_ = b.domainMax_;
    globalMeshIndex_ = b.globalMeshIndex_;
    middle_ = b.middle_;
    
    //A cube has 3^3=27 neighbours
    unsigned long numberOfNeighbours = 27;
    
    //Create all lists
    localIndexToGlobalIndexTable_ = std::vector<int>(numberOfNeighbours);
    localIndexToProcessorList_ = std::vector<int>(numberOfNeighbours);
    boundaryParticleList_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours, std::vector<BaseParticle*>(0));
    boundaryParticleListNeighbour_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours,
                                                                              std::vector<BaseParticle*>(0));
    newBoundaryParticleList_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours,
                                                                        std::vector<BaseParticle*>(0));
    newInteractionList_ = std::vector<std::vector<BaseInteraction*> >(numberOfNeighbours,
                                                                      std::vector<BaseInteraction*>(0));
    numberOfParticlesSend_ = std::vector<unsigned>(numberOfNeighbours);
    numberOfParticlesReceive_ = std::vector<unsigned>(numberOfNeighbours);
    numNewInteractionsSend_ = std::vector<unsigned>(numberOfNeighbours);
    numNewInteractionsReceive_ = std::vector<unsigned>(numberOfNeighbours);
    boundaryParticleDataSend_ = std::vector<std::vector<MPIParticle> >(numberOfNeighbours);
    boundaryParticleDataReceive_ = std::vector<std::vector<MPIParticle> >(numberOfNeighbours);
    updatePositionDataSend_ = std::vector<std::vector<MPIParticlePosition> >(numberOfNeighbours);
    updatePositionDataReceive_ = std::vector<std::vector<MPIParticlePosition> >(numberOfNeighbours);
    updateVelocityDataSend_ = std::vector<std::vector<MPIParticleVelocity> >(numberOfNeighbours);
    updateVelocityDataReceive_ = std::vector<std::vector<MPIParticleVelocity> >(numberOfNeighbours);
    interactionDataSend_ = std::vector<void*>(numberOfNeighbours);
    interactionDataReceive_ = std::vector<void*>(numberOfNeighbours);
    activeBoundaryList_ = std::vector<bool>(numberOfNeighbours, true);
    boundaryList_ = std::vector<int>(0);
#ifdef DEBUG_CONSTRUCTOR
    std::cout<<"Domain::Domain(const Domain &b) finished"<<std::endl;
#endif
}

Domain::~Domain()
{
#ifdef DEBUG_DESTRUCTOR
    std::cout << "Domain::~Domain() finished"<<std::endl;
#endif
}

void Domain::constructor()
{
    rank_ = constants::intMax;
    domainHandler_ = nullptr;
    domainMin_ = {-constants::inf, -constants::inf, -constants::inf};
    domainMax_ = {constants::inf, constants::inf, constants::inf};
    
    //A cube has 3^3=27 neighbours
    int numberOfNeighbours = 27;
    
    //Create all lists
    localIndexToGlobalIndexTable_ = std::vector<int>(numberOfNeighbours);
    localIndexToProcessorList_ = std::vector<int>(numberOfNeighbours);
    boundaryParticleList_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours, std::vector<BaseParticle*>(0));
    boundaryParticleListNeighbour_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours,
                                                                              std::vector<BaseParticle*>(0));
    newBoundaryParticleList_ = std::vector<std::vector<BaseParticle*> >(numberOfNeighbours,
                                                                        std::vector<BaseParticle*>(0));
    newInteractionList_ = std::vector<std::vector<BaseInteraction*> >(numberOfNeighbours,
                                                                      std::vector<BaseInteraction*>(0));
    numberOfParticlesSend_ = std::vector<unsigned>(numberOfNeighbours);
    numberOfParticlesReceive_ = std::vector<unsigned>(numberOfNeighbours);
    numNewInteractionsSend_ = std::vector<unsigned>(numberOfNeighbours);
    numNewInteractionsReceive_ = std::vector<unsigned>(numberOfNeighbours);
    boundaryParticleDataSend_ = std::vector<std::vector<MPIParticle> >(numberOfNeighbours);
    boundaryParticleDataReceive_ = std::vector<std::vector<MPIParticle> >(numberOfNeighbours);
    updatePositionDataSend_ = std::vector<std::vector<MPIParticlePosition> >(numberOfNeighbours);
    updatePositionDataReceive_ = std::vector<std::vector<MPIParticlePosition> >(numberOfNeighbours);
    updateVelocityDataSend_ = std::vector<std::vector<MPIParticleVelocity> >(numberOfNeighbours);
    updateVelocityDataReceive_ = std::vector<std::vector<MPIParticleVelocity> >(numberOfNeighbours);
    interactionDataSend_ = std::vector<void*>(numberOfNeighbours);
    interactionDataReceive_ = std::vector<void*>(numberOfNeighbours);
    activeBoundaryList_ = std::vector<bool>(numberOfNeighbours, true);
    boundaryList_ = std::vector<int>(0);
}

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

void Domain::read(std::istream& is)
{
    logger(WARN, "[Domain::read] should not be called");
    //BaseObject::read(is);
}

void Domain::write(std::ostream& os) const
{
    logger(WARN, "[Domain::write] should not be called");
    //BaseObject::write(os);
}

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

void Domain::setRange(Direction direction, Mdouble min, Mdouble max)
{
    if (min > max)
    {
        logger(ERROR, "[MercuryMPI ERROR]: min is larger than max. (%,%)", min, max);
    }
    
    double maxClosed;
    double minClosed;
    if (min == -constants::inf)
    {
        minClosed = getHandler()->getDPMBase()->getMin().getComponent(direction);
    }
    
    if (max == constants::inf)
    {
        maxClosed = getHandler()->getDPMBase()->getMax().getComponent(direction);
    }
    
    
    switch (direction)
    {
        case Direction::XAXIS :
            domainMin_[0] = min;
            domainMax_[0] = max;
            middle_.X = minClosed + (maxClosed - minClosed) / 2.0;
            break;
        case Direction::YAXIS :
            domainMin_[1] = min;
            domainMax_[1] = max;
            middle_.Y = minClosed + (maxClosed - minClosed) / 2.0;
            break;
        case Direction::ZAXIS :
            domainMin_[2] = min;
            domainMax_[2] = max;
            middle_.Z = minClosed + (maxClosed - minClosed) / 2.0;
            break;
        default :
            logger(ERROR, "Direction is not a valid direction. (%)", direction);
            break;
    }
}

void Domain::setBounds(std::vector<double> domainMin, std::vector<double> domainMax, bool computeMiddle)
{
    domainMin_ = domainMin;
    domainMax_ = domainMax;
    
    //Compute the middle of the closed domain
    if (computeMiddle)
    {
        double minClosed;
        double maxClosed;
        for (int i = 0; i < 3; i++)
        {
            minClosed = domainMin_[i];
            maxClosed = domainMax_[i];
            if (domainMin_[i] == -constants::inf)
            {
                minClosed = getHandler()->getDPMBase()->getMin().getComponent(i);
            }
            
            if (domainMax_[i] == constants::inf)
            {
                maxClosed = getHandler()->getDPMBase()->getMax().getComponent(i);
            }
            
            middle_.setComponent(i, minClosed + (maxClosed - minClosed) / 2.0);
        }
    }
}

void Domain::setRank(int rank)
{
    rank_ = rank;
}

std::vector<double> Domain::getDomainMin()
{
    return domainMin_;
}

std::vector<double> Domain::getDomainMax()
{
    return domainMax_;
}

int Domain::getRank()
{
    return rank_;
}

void Domain::setHandler(DomainHandler* domainHandler)
{
    domainHandler_ = domainHandler;
}

int Domain::getGlobalIndex()
{
    return globalIndex_;
}

std::vector<unsigned> Domain::getGlobalMeshIndex()
{
    return globalMeshIndex_;
}

void Domain::setGlobalMeshIndex(std::vector<unsigned> globalMeshIndex)
{
    globalMeshIndex_ = globalMeshIndex;
}

void Domain::disableBoundary(unsigned localIndex)
{
    activeBoundaryList_[localIndex] = false;
}

std::vector<bool> Domain::getActiveBoundaryList()
{
    return activeBoundaryList_;
}

bool Domain::containsParticle(BaseParticle* particle, Mdouble offset)
{
    for (unsigned i = 0; i < 3; i++)
    {
        if (!(((domainMin_[i] + offset) < particle->getPosition().getComponent(i))
              &&
              ((domainMax_[i] - offset) >= particle->getPosition().getComponent(i)))
                )
        {
            return false;
        }
    }
    return true;
}

bool Domain::isInGreaterDomain(BaseParticle* particle)
{
    return containsParticle(particle, -domainHandler_->getInteractionDistance());
}

bool Domain::isInInnerDomain(BaseParticle* particle)
{
    return containsParticle(particle, domainHandler_->getInteractionDistance());
}

bool Domain::isInCommunicationZone(BaseParticle* particle)
{
    //First check if the particle is actually in the domain
    //Secondly check if the particle is in the inner domain
    return containsParticle(particle) && !(isInInnerDomain(particle));
}

void Domain::createLookUpTable()
{
    int localIndex;
    int globalIndex;
    //Get the global decomposition vector, (nx,ny,nz) and compute the mesh multipliers: (1,nx,nx*ny)
    std::vector<unsigned> numberOfDomains = domainHandler_->getNumberOfDomains();
    std::vector<unsigned> globalMeshMultiplier = {1,
                                                  numberOfDomains[Direction::XAXIS],
                                                  numberOfDomains[Direction::XAXIS] *
                                                  numberOfDomains[Direction::YAXIS]};
    
    //Compute the globalIndex of this domain
    globalIndex_ = globalMeshMultiplier[Direction::XAXIS] * globalMeshIndex_[Direction::XAXIS] +
                   globalMeshMultiplier[Direction::YAXIS] * globalMeshIndex_[Direction::YAXIS] +
                   globalMeshMultiplier[Direction::ZAXIS] * globalMeshIndex_[Direction::ZAXIS];
    
    //Create the lookup table from localIndex to globalIndex of the neihbours
    for (int i = -1; i < 2; i++)
    {
        for (int j = -1; j < 2; j++)
        {
            for (int k = -1; k < 2; k++)
            {
                //Get the local index  
                localIndex = i + 3 * j + 9 * k + 13;
                //Compute the global index
                globalIndex = globalMeshMultiplier[Direction::XAXIS] * (globalMeshIndex_[Direction::XAXIS] + i) +
                              globalMeshMultiplier[Direction::YAXIS] * (globalMeshIndex_[Direction::YAXIS] + j) +
                              globalMeshMultiplier[Direction::ZAXIS] * (globalMeshIndex_[Direction::ZAXIS] + k);
                //Fill in the look-up table
                localIndexToGlobalIndexTable_[localIndex] = globalIndex;
            }
        }
    }
    //Create the lookup table from localIndex to processorRank
    rank_ = globalIndex_;
    localIndexToProcessorList_ = localIndexToGlobalIndexTable_;
}

int Domain::getLocalIndex(const int i, const int j, const int k)
{
    return i + 3 * j + 9 * k + 13;
}

int Domain::getLocalIndex(const std::vector<int> localMeshIndex)
{
    return localMeshIndex[Direction::XAXIS] + 3 * localMeshIndex[Direction::YAXIS] +
           9 * localMeshIndex[Direction::ZAXIS] + 13;
}

std::vector<int> Domain::getLocalIndexInverse(int localIndex) {
    localIndex -= 13;
    return {localIndex%3, (localIndex/3)%3, (localIndex/9)%3};
}

BaseParticle* Domain::findParticleInList(unsigned int identification, std::vector<BaseParticle*> particleList)
{
    for (BaseParticle* particle : particleList)
    {
        if (particle->getId() == identification)
        {
            return particle;
        }
    }
    return nullptr;
}

std::vector<int> Domain::findNearbyBoundaries(BaseParticle* particle, Mdouble offset)
{
    std::vector<int> boundaryIndex(3);
    Mdouble interactionDistance = domainHandler_->getInteractionDistance();
    
    //for x,y,z directions, find if the particle is close to the left or right wall or not at all
    for (int d = 0; d < 3; d++)
    {
        //Check if the particle is close to Lx
        if ((particle->getPosition().getComponent(d)) < (domainMin_[d] + interactionDistance + offset))
        {
            //Update switch
            boundaryIndex[d] = -1;
        }
            //Check if particle is close to Rx
        else if ((particle->getPosition().getComponent(d)) >= (domainMax_[d] - interactionDistance - offset))
        {
            //Update switch
            boundaryIndex[d] = 1;
        }
    }
    return boundaryIndex;
}

bool Domain::inBoundary(BaseParticle* particle, int localIndex_) {
    std::vector<int> boundaryIndex = findNearbyBoundaries(particle);
    std::vector<int> localIndex = getLocalIndexInverse(localIndex_);
    // if the particle is in the boundary zone into which we want to transmit
    for (int i = 0; i < 3; ++i)
    {
        if (localIndex[i]!=0 && boundaryIndex[i]!=localIndex[i]) {
            return false;
        }
    }
    return true;
}

void Domain::disableBoundaries()
{
    std::vector<unsigned> numberOfDomains = domainHandler_->getNumberOfDomains();
    std::vector<int> localMeshIndex(3);
    int localIndex;
    
    //disble own list
    activeBoundaryList_[13] = false;
    
    //Special case when numberOfDomains[d] = 1 -> disable all d direction boundaries
    for (unsigned d = 0; d < 3; d++) //Loop over all dimensions
    {
        //If there is only one domain in the d-direction, take action
        if (numberOfDomains[d] == 1)
        {
            //Loop over all locaIndices
            for (int i = -1; i < 2; i++)
            {
                localMeshIndex[Direction::XAXIS] = i;
                for (int j = -1; j < 2; j++)
                {
                    localMeshIndex[Direction::YAXIS] = j;
                    for (int k = -1; k < 2; k++)
                    {
                        localMeshIndex[Direction::ZAXIS] = k;
                        //Disable all boundaries that are not in the "middle" of the mesh in the d-direction
                        // i.e. a localMeshIndex[d] = -1 is a neighbour in the d-direction, however there is only one
                        // element in that direction so it can't be a real neighbour.
                        if (localMeshIndex[d] != 0)
                        {
                            //Disable the boundary
                            localIndex = getLocalIndex(localMeshIndex);
                            activeBoundaryList_[localIndex] = false;
                        }
                    }
                }
            }
        }
    }
    
    //Compute the boundaryIndex for all other cases where numberOfDomains is larger than 1
    std::vector<int> boundaryIndex(3);
    for (int d = 0; d < 3; d++)
    {
        //Check if the domain is on the simulationDomainMin boundary
        if (globalMeshIndex_[d] == 0)
        {
            //Update boundary index
            boundaryIndex[d] = -1;
        }
            //Check if the domain is on the simulationDomainMax boundary
        else if (globalMeshIndex_[d] == (numberOfDomains[d] - 1))
        {
            //Update boundary index
            boundaryIndex[d] = 1;
        }
    }
    
    //Disable the boundaries that are on the edge of the simulation domain
    for (int d = 0; d < 3; d++)
    {
        //If boundary is on the edge of the simulation domain
        if (boundaryIndex[d] != 0)
        {
            //Loop over all local indices
            for (int i = -1; i < 2; i++)
            {
                for (int j = -1; j < 2; j++)
                {
                    for (int k = -1; k < 2; k++)
                    {
                        localMeshIndex = {i, j, k};
                        //Set the localMeshIndex to the boundaryIndex in the d-direction
                        //Since we are only interested in the domains on this boundary
                        localMeshIndex[d] = boundaryIndex[d];
                        //Disable the boundary
                        localIndex = getLocalIndex(localMeshIndex);
                        activeBoundaryList_[localIndex] = false;
                    }
                }
            }
        }
    }
    
    //Create a list of the active boundaries
    localIndex = 0;
    for (bool active : activeBoundaryList_)
    {
        if (active)
        {
            boundaryList_.push_back(localIndex);
        }
        localIndex++;
    }
}

void Domain::addParticlesToLists(BaseParticle* particle, std::vector<std::vector<BaseParticle*> >& list)
{
    std::vector<int> boundaryIndex = findNearbyBoundaries(particle);
   
    //Compute and set complexity of the particle
    unsigned int complexity = boundaryIndex[0] + 3 * boundaryIndex[1] + 9 * boundaryIndex[2] + 13;
    unsigned int list_complexity = 0;
    for (int d = 0; d < 3; d++) //Loop over all directions
    {
        list_complexity += std::abs(boundaryIndex[d]);
    }
    particle->setCommunicationComplexity(complexity);
    //particle->setCommunicationComplexity(list_complexity);
    
    //Based on the complexity of the particle, add it to the approriate list
    switch (list_complexity)
    {
        //The particle is not close at all
        case 0:
            break;
            //The particle is close to one side
            // 1 side contribution
        case 1 :
            //Add the side contribution
            list[getLocalIndex(boundaryIndex)].push_back(particle);
            break;
            //The particle is close to two neighbouring directions
            //2 side and 1 rib contrubution
        case 2 :
        {
            //Add the two side contributions
            for (int d = 0; d < 3; d++)
            {
                std::vector<int> localMeshIndex = {0, 0, 0};
                localMeshIndex[d] = boundaryIndex[d];
                //Avoid adding the particle in the wrong direction by excluding localMeshIndex[d] = 0
                if (localMeshIndex[d] != 0)
                {
                    list[getLocalIndex(localMeshIndex)].push_back(particle);
                }
            }
        }
            //Add the rib contribution
            list[getLocalIndex(boundaryIndex)].push_back(particle);
            break;
            
            //The particle is close to three neighbouring directions
            // 3 side, 3 rib and 1 corner contribution
        case 3 :
        {
            //Add the three side contributions
            for (int d = 0; d < 3; d++)
            {
                std::vector<int> localMeshIndex = {0, 0, 0};
                //Reset index vector
                localMeshIndex[d] = boundaryIndex[d];
                list[getLocalIndex(localMeshIndex)].push_back(particle);
            }
            
            //Add the three rib contributions
            for (int d = 0; d < 3; d++)
            {
                std::vector<int> localMeshIndex = boundaryIndex;
                //All rib boundary indices are given by the boundaryIndex and setting one of the components to zero
                localMeshIndex[d] = 0;
                list[getLocalIndex(localMeshIndex)].push_back(particle);
            }
        }
            
            //Add the corner contribution
            list[getLocalIndex(boundaryIndex)].push_back(particle);
            break;
        
        default :
            logger(INFO, "boundaryIndex : %,%,% | list_complexity: %", boundaryIndex[0], boundaryIndex[1], boundaryIndex[2],
                   list_complexity);
            logger(ERROR, "Particle is in contact with the wrong number of boundaries");
            break;
    }
}

void Domain::findNewMPIParticles(const ParticleHandler& particleHandler)
{
    //For all particles inside the given domain, loop over all the particles to
    //see if we have to add the particles to the new particle list
    for (BaseParticle* particle : particleHandler)
    {
        findNewMPIParticle(particle);
    }
}

void Domain::findNewMPIParticle(BaseParticle* particle)
{
    //If the particle is an MPI particle
    if (!particle->isInMPIDomain() && !particle->isPeriodicGhostParticle())
    {
        addParticlesToLists(particle, newBoundaryParticleList_);
    }
}

bool Domain::isInNewBoundaryParticleList(BaseParticle* objectI,int localIndex) const
{
    for (BaseParticle *q : newBoundaryParticleList_[localIndex]) {
        if (q == objectI)
            return true;
    }
    return false;
}

void Domain::findNewMPIInteractions()
                            {
                                //For all active boundaries
                                for (int localIndex : boundaryList_)
                                {
                                    //Check all newly added particles for interactions with already known particles
                                    for (BaseParticle* newBoundaryParticle : newBoundaryParticleList_[localIndex])
                                    {
                                        //Loop over all interactions of the new particle
                                        std::vector<BaseInteraction*> interactions = newBoundaryParticle->getInteractions();
                                        for (BaseInteraction* interaction : interactions)
                                        {
                                            //Find out what the new particle is interacting with
                                            BaseParticle* particleP = dynamic_cast<BaseParticle*>(interaction->getP());
                                            BaseParticle* objectI = dynamic_cast<BaseParticle*>(interaction->getI());

                                            //If the P in the interaction structure is the new particle, find I
                                            if (newBoundaryParticle == particleP)
                                            {
                                                //Check if the new particle is interacting with a wall
                                                if (!objectI)
                                                {
                                                    //Check if the interaction is still valid
                                                    BaseWall* wall = dynamic_cast<BaseWall*>(interaction->getI());
                                                    Mdouble distance;
                                                    Vec3D normal;
                                                    wall->getDistanceAndNormal(*particleP, distance, normal);
                                                    if (distance <= particleP->getMaxInteractionRadius())
                                                    {
                                                        //Add the interaction to the list if there are still in contact (hence the if statement)
                                                        newInteractionList_[localIndex].push_back(interaction);
                        }
                    }
                    else //is I a particle
                    {
                        //check if particle is in mpi domain OR (TODO) if it is in the newBoundaryParticleList_
                        if (objectI->isInMPIDomain() || isInNewBoundaryParticleList(objectI, localIndex))
                        {
                            // iff the interaction partner of the new particle is in the communication zone, transmit the interaction
                            if (inBoundary(objectI,localIndex))
                            {
                                //Is the particle still interacting after the position update?
                                Vec3D branchVector = particleP->getPosition() - objectI->getPosition();
                                //Get the square of the distance between particle i and particle j
                                Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
                                Mdouble sumOfInteractionRadii =
                                    objectI->getSumOfInteractionRadii(particleP);
                                if (distanceSquared < (sumOfInteractionRadii * sumOfInteractionRadii))
                                {
                                    //Add the interaction to the list
                                    newInteractionList_[localIndex].push_back(interaction);
                                }
                            }
                        }
                    }
                }
                else //newBoundaryParticle is I in the interaction, P can only be a particle
                {
                    //it is "I" in the interaction structure, check if its in the mpi domain
                    if (particleP->isInMPIDomain() || isInNewBoundaryParticleList(particleP, localIndex))
                    {
                        //Is the particle still interacting after the position update?
                        Vec3D branchVector = particleP->getPosition() - objectI->getPosition();
                        //Get the square of the distance between particle i and particle j
                        Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
                        Mdouble sumOfInteractionRadii =
                                objectI->getSumOfInteractionRadii(particleP);
                        if (distanceSquared < (sumOfInteractionRadii * sumOfInteractionRadii))
                        {
                            //Add the interaction to the list
                            newInteractionList_[localIndex].push_back(interaction);
                        }
                    }
                }
            }
        }
    }
}

void Domain::collectBoundaryParticleData(int localIndex)
{
    //For all particles
    for (BaseParticle* particle : newBoundaryParticleList_[localIndex])
    {
        //Add the data to the transmission list
        boundaryParticleDataSend_[localIndex].push_back(copyDataFromParticleToMPIParticle(particle));
    }
}

void Domain::collectInteractionData(int localIndex)
{
    //Copy interactions to the data array
    unsigned int indexInteraction = 0;
    for (BaseInteraction* interaction : newInteractionList_[localIndex])
    {
        interaction->getMPIInteraction(interactionDataSend_[localIndex], indexInteraction);
        indexInteraction++;
    }
}

void Domain::processReceivedBoundaryParticleData(const unsigned index, std::vector<BaseParticle*>& newParticles)
{
    
    ParticleHandler& particleHandler = getHandler()->getDPMBase()->particleHandler;
    //for all MPIPparticles that have been received on the other side
    for (int i = 0; i < numberOfParticlesReceive_[index]; i++)
    {
        //copy data from data
        BaseParticle* p0 = MPIParticle::newParticle();
        copyDataFromMPIParticleToParticle(&(boundaryParticleDataReceive_[index][i]), p0, &particleHandler);
        
        //Flag that the particle is a ghost particle of a real particle in the neighbour domain
        p0->setMPIParticle(true);
        //Flag that the particle is list as a particle in the mpi domain
        p0->setInMPIDomain(true);
        //add particle to handler
        particleHandler.addGhostObject(p0);
        //Set previous position as current position
        BaseParticle* pGhost = particleHandler.getLastObject();
        pGhost->setPreviousPosition(particleHandler.getLastObject()->getPosition());
        //Add to the newParticle list
        newParticles.push_back(pGhost);
        logger(VERBOSE, "Adding mpi particle % at: %", particleHandler.getLastObject()->getId(),
               particleHandler.getLastObject()->getPosition());
        
        //add pointer to the particle to the list and set the status to idle
        boundaryParticleListNeighbour_[index].push_back(particleHandler.getLastObject());
    }
}

void Domain::processSentBoundaryParticles(const unsigned index)
{
    for (BaseParticle* particle : newBoundaryParticleList_[index])
    {
        boundaryParticleList_[index].push_back(particle);
        particle->setInMPIDomain(true);
    }
}

void Domain::processReceivedInteractionData(const unsigned localIndex, std::vector<BaseParticle*>& newParticles)
{
    InteractionHandler& iH = getHandler()->getDPMBase()->interactionHandler;
    for (unsigned int l = 0; l < numNewInteractionsReceive_[localIndex]; l++)
    {
        unsigned int identificationP;
        unsigned int identificationI;
        bool isWallInteraction;
        unsigned timeStamp;
        
        //Get the general information required to setup a new interaction
        iH.getInteractionDetails(interactionDataReceive_[localIndex],
                                 l, identificationP, identificationI,
                                 isWallInteraction, timeStamp);
        
        logger(VERBOSE, "interaction details: %, %, idP %, idI %, wall %, time %", interactionDataReceive_[localIndex],
               l, identificationP, identificationI,
               isWallInteraction, timeStamp);
        
        //Find the particle in the newParticle list
        BaseParticle* pGhost = nullptr;
        int idOther;
        for (BaseParticle* particle : newParticles)
        {
            if (particle->getId() == identificationP)
            {
                pGhost = particle;
                idOther = identificationI;
                break;
            }
            
            if (particle->getId() == identificationI)
            {
                pGhost = particle;
                idOther = identificationP;
                break;
            }
        }
        if (pGhost == nullptr)
        {
            logger(WARN, "In Domain::processReceivedInteractionData: pGhost (id %) is nullptr, the interaction data is not copied. Two particles possibly moved into domain simultaneously.", identificationP);
            continue;
        }
        
        //If it is a wall interaction, do stuff
        if (isWallInteraction)
        {
            BaseInteractable* I = getHandler()->getDPMBase()->wallHandler.getObjectById(identificationI);
            //Create interactions
            BaseInteraction* j = I->getInteractionWith(pGhost, timeStamp, &iH);
            if (j!= nullptr) j->setMPIInteraction(interactionDataReceive_[localIndex], l, false);
            
        }
        else
        {
            //Obtain potential interaction particles
            std::vector<BaseParticle*> interactingParticleList;
            getHandler()->getDPMBase()->hGridGetInteractingParticleList(pGhost, interactingParticleList);
            
            //Find the other interacting particle
            BaseParticle* otherParticle = nullptr;
            for (BaseParticle* p2 : interactingParticleList)
            {
                if (p2->getId() == idOther)
                {
                    otherParticle = p2;
                    break;
                }
            }
            if (otherParticle == nullptr) {
                //search for the other particle in the newParticles list
                for (BaseParticle *particle : newParticles) {
                    if (particle->getId() == idOther) {
                        otherParticle = particle;
                        break;
                    }
                }
                if (otherParticle == nullptr) {
                    logger(WARN,
                           "In Domain::processReceivedInteractionData: otherParticle (id %) is nullptr, the interaction data with pGhost (id %) is not copied. Two particles possibly moved into domain simultaneously. nt = %",
                           identificationI, identificationP, getHandler()->getDPMBase()->getNumberOfTimeSteps());
                    continue;
                }
            }
            //Add the interaction
            // flipped order of other and ghost, so it is in the same order as received
            BaseInteraction* j;
            if (otherParticle->getId() < pGhost->getId()) {
                j = pGhost->getInteractionWith(otherParticle, timeStamp, &iH);
            } else {
                j = otherParticle->getInteractionWith(pGhost, timeStamp, &iH);
            }
            if (j != nullptr) j->setMPIInteraction(interactionDataReceive_[localIndex], l, false);
        }
    }
}

void Domain::debugInformation() {
    std::stringstream s;
    std::stringstream m;
    for (auto p : getHandler()->getDPMBase()->particleHandler) {
        if (p->isMPIParticle())
            m << std::setw(4) << p->getId();
        else
            s << std::setw(4) << p->getId();
    }
    logger(INFO,"Particles %, MPI %",s.str(),m.str());
}


        void Domain::sendAndReceiveCount(MercuryMPITag tag, unsigned& countReceive, unsigned& countSend,
                unsigned localIndexNeighbour)
        {
            int globalIndexNeighbour = localIndexToGlobalIndexTable_[localIndexNeighbour];
            int processor = localIndexToProcessorList_[localIndexNeighbour];

            //Create communication tags
            int tagReceive = globalIndexNeighbour * MAX_PROC + globalIndex_ * 10 + tag;
            int tagSend = globalIndex_ * MAX_PROC + globalIndexNeighbour * 10 + tag;

            logger.assert_debug(tagSend > 0, "Send tag is wrong. Tag: %", tagSend);
            logger.assert_debug(tagReceive > 0, "Receive tag is wrong. Tag: %", tagReceive);
            logger.assert_debug(processor >= 0, "Processor is wrong. processor: %", processor);

            //Communicate the requests
            MPIContainer::Instance().receive(countReceive, processor, tagReceive);
            MPIContainer::Instance().send(countSend, processor, tagSend);
        }

        void Domain::prepareBoundaryDataTransmission()
        {
            //Find new particles that have entered the communication zone
            findNewMPIParticles(getHandler()->getDPMBase()->particleHandler);

            //Find interactions between new particles and other particles in the communication zone
            findNewMPIInteractions();

            //Compute number of particles to send
            //For all boundaries
            for (int localIndex : boundaryList_)
            {
                numberOfParticlesSend_[localIndex] = newBoundaryParticleList_[localIndex].size();
                numNewInteractionsSend_[localIndex] = newInteractionList_[localIndex].size();
            }

            //For all active boundaries
            for (int localIndex : boundaryList_)
            {
                //Send and receive the number of new boundary particles
                sendAndReceiveCount(MercuryMPITag::PARTICLE_COUNT, (numberOfParticlesReceive_[localIndex]),
                                    (numberOfParticlesSend_[localIndex]), localIndex);

                //Send and receive the new interactions
                sendAndReceiveCount(MercuryMPITag::INTERACTION_COUNT, (numNewInteractionsReceive_[localIndex]),
                                    (numNewInteractionsSend_[localIndex]), localIndex);
            }
        }

        void Domain::prepareBoundaryDataTransmission(BaseParticle* particle)
        {
            //Assign the particle to the correct lists if the particle belongs to this domain
            if (containsParticle(particle))
            {
                findNewMPIParticle(particle);
            }

            //Note: When inserting a single particle, it has no interactions, so no interactions to be copied either

            //Compute number of particles to send
            //For all boundaries
            for (int localIndex : boundaryList_)
            {
                numberOfParticlesSend_[localIndex] = newBoundaryParticleList_[localIndex].size();
                numNewInteractionsSend_[localIndex] = newInteractionList_[localIndex].size();
            }


            //For all active boundaries
            for (int localIndex : boundaryList_)
            {
                //Send and recieve the number of new boundary particles
                sendAndReceiveCount(MercuryMPITag::PARTICLE_COUNT, (numberOfParticlesReceive_[localIndex]),
                                    (numberOfParticlesSend_[localIndex]), localIndex);

                //Send and receive the new interactions
                sendAndReceiveCount(MercuryMPITag::INTERACTION_COUNT, (numNewInteractionsReceive_[localIndex]),
                                    (numNewInteractionsSend_[localIndex]), localIndex);
            }
        }

        void Domain::performBoundaryDataTransmission()
        {
            //For all active boundaries
            for (int localIndex : boundaryList_)
            {
                //make sure enough memory is reserved to receive the data
                boundaryParticleDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);

                //Collect the particle data
                collectBoundaryParticleData(localIndex);

                //Send the data
                sendAndReceiveMPIData(MercuryMPITag::PARTICLE_DATA, MercuryMPIType::PARTICLE,
                                      boundaryParticleDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
                                      boundaryParticleDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
                                      localIndex);


                //create arrays for sending and receiving interaction data
                interactionDataSend_[localIndex] = getHandler()->getDPMBase()->interactionHandler.createMPIInteractionDataArray(
                        numNewInteractionsSend_[localIndex]);
                interactionDataReceive_[localIndex] = getHandler()->getDPMBase()->interactionHandler.createMPIInteractionDataArray(
                        numNewInteractionsReceive_[localIndex]);

                //Collect the interaction data
                collectInteractionData(localIndex);

                //Send the data
                sendAndReceiveMPIData(MercuryMPITag::INTERACTION_DATA, MercuryMPIType::INTERACTION,
                                      interactionDataReceive_[localIndex], numNewInteractionsReceive_[localIndex],
                                      interactionDataSend_[localIndex], numNewInteractionsSend_[localIndex], localIndex);
            }
        }

        void Domain::finaliseBoundaryDataTransmission()
        {
            //For all boundaries
            for (int localIndex : boundaryList_)
            {
                //copy the received data into the particleHandler and neighbour particle list
                std::vector<BaseParticle*> newParticles;
                processReceivedBoundaryParticleData(localIndex, newParticles);
                //All particles in the current domain that have been send to the other domains need to be flagged as communicating particles
                processSentBoundaryParticles(localIndex);
                //copy the received interaction data into the standard dpm structure
                processReceivedInteractionData(localIndex, newParticles);

                //Delete all send/receive data
                boundaryParticleDataSend_[localIndex].clear();
                boundaryParticleDataReceive_[localIndex].clear();
                getHandler()->getDPMBase()->interactionHandler.deleteMPIInteractionDataArray(interactionDataSend_[localIndex]);
                getHandler()->getDPMBase()->interactionHandler.deleteMPIInteractionDataArray(
                        interactionDataReceive_[localIndex]);

                //Reset all counters
                numberOfParticlesSend_[localIndex] = 0;
                numberOfParticlesReceive_[localIndex] = 0;
                numNewInteractionsSend_[localIndex] = 0;
                numNewInteractionsReceive_[localIndex] = 0;

                //Reset all lists
                newBoundaryParticleList_[localIndex].clear();
                newInteractionList_[localIndex].clear();
            }
        }

        void Domain::updateParticles(std::set<BaseParticle*>& ghostParticlesToBeDeleted)
        {
            int complexityNew;
            std::vector<int> boundaryIndex;

            //For all active boundaries
            for (int localIndex : boundaryList_)
            {
                //Step 1A: Remove the particles from the boundaryParticleList_ which require re-assignment
                //or have just moved away from the region
                for (int p = 0; p < boundaryParticleList_[localIndex].size(); p++)
                {
                    BaseParticle* particle = boundaryParticleList_[localIndex][p];
                    //Check if the particle is still in the domain, but not in the communication zone
                    if (containsParticle(particle))
                    {
                        //check if the complexity has changed
                        boundaryIndex = findNearbyBoundaries(particle);
                        complexityNew = boundaryIndex[0] + 3 * boundaryIndex[1] + 9 * boundaryIndex[2] + 13;
                        if (particle->getCommunicationComplexity() != complexityNew)
                        {
                            logger(VERBOSE, "time: % | global index: % in list % | particle % | CURRENT DOMAIN - CHANGES "
                                            "COMPLEXITY", domainHandler_->getDPMBase()->getTime(), globalIndex_,
                                   localIndexToGlobalIndexTable_[localIndex], particle->getId());
                            //Flag the particle that it no longer participates in the communication layer
                            //so it can be reintroduced in the transmission step
                            particle->setMPIParticle(false);
                            particle->setInMPIDomain(false);
                            particle->setCommunicationComplexity(0);
                            boundaryParticleList_[localIndex][p] = nullptr;
                        }
                    }
                    else
                    {
                        logger(VERBOSE, "time: % | global index: % in list % | particle % | CURRENT DOMAIN - TO NEIGHBOURING "
                                        "DOMAIN", domainHandler_->getDPMBase()->getTime(), globalIndex_,
                               localIndexToGlobalIndexTable_[localIndex], particle->getId());

                        ghostParticlesToBeDeleted.insert(particle);
                        boundaryParticleList_[localIndex][p] = nullptr;
                    }
                }

                //Step 1B: Remove the particles from the boundaryParticleListNeightbour_ which require re-assignment
                //or have just moved away from the region
                for (int p = 0; p < boundaryParticleListNeighbour_[localIndex].size(); p++)
                {
                    BaseParticle* particle = boundaryParticleListNeighbour_[localIndex][p];
                    //Check if the particle moved out of the neighbour domain
                    if (!(domainHandler_->getObject(localIndexToGlobalIndexTable_[localIndex])->containsParticle(particle)))
                    {
                        //The particle has moved to this domain
                        if (containsParticle(particle))
                        {
                            logger(VERBOSE,
                                   "time: % | global index: % in list % | particle % | NEIGHBOURING DOMAIN - TO CURRENT DOMAIN",
                                   domainHandler_->getDPMBase()->getTime(), globalIndex_,
                                   localIndexToGlobalIndexTable_[localIndex], particle->getId());

                            //Flag the particle as not yet communicated
                            particle->setMPIParticle(false);
                            particle->setInMPIDomain(false);
                            particle->setCommunicationComplexity(0);
                            boundaryParticleListNeighbour_[localIndex][p] = nullptr;
                        }
                            //The particle has moved to a different domain
                        else
                        {
                            logger(VERBOSE,
                                   "time: % | global index: % in list % | particle % | NEIGBOURING DOMAIN - TO OTHER DOMAIN",
                                   domainHandler_->getDPMBase()->getTime(), globalIndex_,
                                   localIndexToGlobalIndexTable_[localIndex], particle->getId());
                            //Cruelly destroy the particle without any mercy.
                            ghostParticlesToBeDeleted.insert(particle);
                            boundaryParticleListNeighbour_[localIndex][p] = nullptr;
                        }
                    }
                    else
                    {
                        //check if the complexity has changed
                        boundaryIndex = domainHandler_->getObject(
                                localIndexToGlobalIndexTable_[localIndex])->findNearbyBoundaries(particle);
                        complexityNew = boundaryIndex[0] + 3 * boundaryIndex[1] + 9 * boundaryIndex[2] + 13;
                        if (particle->getCommunicationComplexity() != complexityNew)
                        {
                            logger(VERBOSE,
                                   "time: % | global index: % in list % | particle % | NEIGHBOURING DOMAIN - CHANGES COMPLEXITY",
                                   domainHandler_->getDPMBase()->getTime(), globalIndex_,
                                   localIndexToGlobalIndexTable_[localIndex], particle->getId());
                            //Cruelly destroy the particle without any mercy.
                            ghostParticlesToBeDeleted.insert(particle);
                            boundaryParticleListNeighbour_[localIndex][p] = nullptr;
                        }
                    }
                }
            }
        }

        void Domain::updateParticlePosition(int localIndex)
        {
            //process the updated information
            unsigned int index = 0;
            for (BaseParticle* particle : boundaryParticleListNeighbour_[localIndex])
            {
                logger.assert_debug(particle->getId() == updatePositionDataReceive_[localIndex][index].id,
                              "MPI particle lists are not in sync");

                //set position
        particle->setPreviousPosition(particle->getPosition());
        particle->setPosition(updatePositionDataReceive_[localIndex][index].position);
        particle->setOrientation(updatePositionDataReceive_[localIndex][index].orientation);
        if (std::is_base_of<MPILiquidFilmParticle,MPIParticle>())
            static_cast<LiquidFilmParticle*>(particle)->setLiquidVolume(updatePositionDataReceive_[localIndex][index].liquidVolume);

        //Update hGrid
        Vec3D displacement = particle->getPreviousPosition() - particle->getPosition();
        getHandler()->getDPMBase()->hGridUpdateMove(particle, displacement.getLengthSquared());
        
        index++;
    }
}

void Domain::updateParticleVelocity(int localIndex)
{
    //process the updated information
    unsigned int index = 0;
    for (BaseParticle* particle: boundaryParticleListNeighbour_[localIndex])
    {
        //set velocity
        particle->setVelocity(updateVelocityDataReceive_[localIndex][index].velocity);
        particle->setAngularVelocity(updateVelocityDataReceive_[localIndex][index].angularVelocity);
        index++;
    }
}

void Domain::preparePositionAndVelocityUpdate()
{
    //For all active boundaries
    for (int localIndex : boundaryList_)
    {
        numberOfParticlesSend_[localIndex] = boundaryParticleList_[localIndex].size();
        numberOfParticlesReceive_[localIndex] = boundaryParticleListNeighbour_[localIndex].size();
        
        //Increase capacity for the receiving data files
        updatePositionDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
        updateVelocityDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
        
        
        //Collect data
        for (BaseParticle* particle : boundaryParticleList_[localIndex])
        {
            updatePositionDataSend_[localIndex].push_back(copyPositionFrom(particle));
            updateVelocityDataSend_[localIndex].push_back(copyVelocityFrom(particle));
        }
        
        //Send and receive the data
        sendAndReceiveMPIData(MercuryMPITag::POSITION_DATA, MercuryMPIType::POSITION,
                              updatePositionDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
                              updatePositionDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
                              localIndex);
        sendAndReceiveMPIData(MercuryMPITag::VELOCITY_DATA, MercuryMPIType::VELOCITY,
                              updateVelocityDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
                              updateVelocityDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
                              localIndex);
    }
}

void Domain::finalisePositionAndVelocityUpdate(std::set<BaseParticle*>& ghostParticlesToBeDeleted)
{
    
    //For all active boundaries
    for (int localIndex : boundaryList_)
    {
        updateParticlePosition(localIndex);
        updateParticleVelocity(localIndex);
    }
    
    //Based on the new position, update the particle lists. 
    //Remove particles that left the communication zone, they will be re-communicated in a later step
    updateParticles(ghostParticlesToBeDeleted);
    
    //For all active boundaries clear the data lists
    for (int localIndex : boundaryList_)
    {
        updatePositionDataSend_[localIndex].clear();
        updateVelocityDataSend_[localIndex].clear();
        updatePositionDataReceive_[localIndex].clear();
        updateVelocityDataReceive_[localIndex].clear();
    }
    
}

void Domain::prepareVelocityUpdate()
{
    //For all active boundaries
    for (int localIndex : boundaryList_)
    {
        numberOfParticlesSend_[localIndex] = boundaryParticleList_[localIndex].size();
        numberOfParticlesReceive_[localIndex] = boundaryParticleListNeighbour_[localIndex].size();
        
        //Resize the vector to the correct size
        updateVelocityDataReceive_[localIndex].resize(numberOfParticlesReceive_[localIndex]);
        
        //Collect data
        for (BaseParticle* particle : boundaryParticleList_[localIndex])
        {
            updateVelocityDataSend_[localIndex].push_back(copyVelocityFrom(particle));
        }
        
        //Send the data
        sendAndReceiveMPIData(MercuryMPITag::VELOCITY_DATA, MercuryMPIType::VELOCITY,
                              updateVelocityDataReceive_[localIndex].data(), numberOfParticlesReceive_[localIndex],
                              updateVelocityDataSend_[localIndex].data(), numberOfParticlesSend_[localIndex],
                              localIndex);
    }
}

void Domain::finaliseVelocityUpdate()
{
    //For all active boundaries
    for (int localIndex : boundaryList_)
    {
        updateParticleVelocity(localIndex);
    }
    
    //For all active boundaries clear the data lists
    for (int localIndex : boundaryList_)
    {
        updateVelocityDataSend_[localIndex].clear();
        updateVelocityDataReceive_[localIndex].clear();
    }
}


DomainHandler* Domain::getHandler() const
{
    return domainHandler_;
}

void Domain::addNewParticles()
{
    //Step 1: For every MPIDomain boundary, create a list of particles that have to be transmitted
    //queue send and receive instructions for the number of particles
    prepareBoundaryDataTransmission();
    MPIContainer::Instance().sync();
    
    //Step 2: queue send and receive of data
    performBoundaryDataTransmission();
    MPIContainer::Instance().sync();
    
    //Step 3: Add the received particles to the particleHandler of the current domain
    finaliseBoundaryDataTransmission();
}

void Domain::addParticle(BaseParticle* particle)
{
    //Step1: check if the particle has to be sent to other processors
    prepareBoundaryDataTransmission(particle);
    MPIContainer::Instance().sync();
    
    //Step2: queue send and receive data. Note for an inserted particle, no interactions should be required
    performBoundaryDataTransmission();
    MPIContainer::Instance().sync();
    
    //Step3: Add the received particles to the particleHandler of the current domain
    finaliseBoundaryDataTransmission();
    MPIContainer::Instance().sync();
}

void Domain::updateStatus(std::set<BaseParticle*>& ghostParticlesToBeDeleted)
{
    //Collect new positions and velocities and send them to the other domains
    preparePositionAndVelocityUpdate();
    MPIContainer::Instance().sync();
    
    //Receive the new positions and velocities from other domains
    //and update the mpi flagged particles accordingly. removes and switched particles in the lists
    finalisePositionAndVelocityUpdate(ghostParticlesToBeDeleted);
    MPIContainer::Instance().sync();
}

void Domain::updateVelocity()
{
    //collect new velocity data and send
    prepareVelocityUpdate();
    MPIContainer::Instance().sync();
    
    //process the received data
    finaliseVelocityUpdate();
    MPIContainer::Instance().sync();
}

unsigned int Domain::getNumberOfMPIParticles()
{
    unsigned int count = 0;
    for (auto& index : boundaryParticleListNeighbour_)
    {
        count += index.size();
    }
    return count;
}

unsigned int Domain::getNumberOfTrueMPIParticles()
{
    unsigned int count = 0;
    for (auto& index : boundaryParticleListNeighbour_)
    {
        for (auto& p : index)
        {
            if (!p->isPeriodicGhostParticle())
            {
                count++;
            }
        }
    }
    return count;
}

void Domain::flushParticles(std::set<BaseParticle*>& toBeFlushedList)
{
    //For all active boundaries
    for (int localIndex : boundaryList_)
    {
        flushParticlesFromList(boundaryParticleList_[localIndex], toBeFlushedList);
        flushParticlesFromList(boundaryParticleListNeighbour_[localIndex], toBeFlushedList);
    }
}

void Domain::flushParticlesFromList(std::vector<BaseParticle*>& list, std::set<BaseParticle*>& toBeFlushedList)
{
    //Firstly: turn all particles that need to be flushed into nullptrs
    for (auto& p : list)
    {
        if (p != nullptr)
        {
            BaseParticle* particle1 = p;
            for (BaseParticle* particle2 : toBeFlushedList)
            {
                //If the particle was found in the list, make a nullptr
                if (particle1 == particle2)
                {
                    logger(VERBOSE, "Removing particle from mpi domain at: %", particle1->getPosition());
                    p = nullptr;
                }
            }
        }
    }
}

Vec3D Domain::getMiddle() const
{
    return middle_;
}

void Domain::cleanCommunicationLists()
{
    //For all active boundaries
    for (int i : boundaryList_)
    {
        cleanCommunicationList(boundaryParticleList_[i]);
        cleanCommunicationList(boundaryParticleListNeighbour_[i]);
    }
}

void Domain::cleanCommunicationList(std::vector<BaseParticle*>& list)
{
    for (int i = 0; i < list.size(); i++)
    {
        if (list[i] == nullptr)
        {
            list[i] = list.back();
            list.pop_back();
            i--;
        }
    }
}