InteractionHandler.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 <iostream>
#include <iomanip>
#include <string>
#include <cstdlib>
#include "Particles/BaseParticle.h"
#include "InteractionHandler.h"
#include "SpeciesHandler.h"
#include "Species/BaseSpecies.h"
#include "Interactions/BaseInteraction.h"
#include "DPMBase.h"
#include "MpiDataClass.h"
#ifdef MERCURY_USE_OMP
#include "omp.h"
#endif

InteractionHandler::InteractionHandler()
{
    writeVTK_ = FileType::NO_FILE;
    logger(DEBUG, "InteractionHandler::InteractionHandler() finished");
}

InteractionHandler::InteractionHandler(const InteractionHandler& IH UNUSED)
        : BaseHandler<BaseInteraction>()
{
    writeVTK_ = IH.writeVTK_;
    //By default interactions are not copied.
    logger(DEBUG, "InteractionHandler::InteractionHandler(const "
                  "InteractionHandler &IH) finished, please note that no interactions"
                  " have been copied.");
}

InteractionHandler& InteractionHandler::operator=(const InteractionHandler& rhs)
{
    if (this != &rhs)
    {
        clear();
    }
    writeVTK_ = rhs.writeVTK_;
    logger(DEBUG, "InteractionHandler::operator =(const InteractionHandler& rhs) finished.");
    return *this;
}

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

void InteractionHandler::addObject(BaseInteraction* I)
{
    //set the particleHandler pointer
    I->setHandler(this);
    //Puts the interaction in the Interaction list
    BaseHandler<BaseInteraction>::addObject(I);
    I->getP()->addInteraction(I);
    I->getI()->addInteraction(I);
}

BaseInteraction*
InteractionHandler::getExistingInteraction(const BaseInteractable* const P, const BaseInteractable* const I) const
{
    //for particle-particle collision it is assumed BaseInteractable P has a lower index then I, so we only have to check for I, not P
    {
        for (unsigned j = 0; j < P->getInteractions().size(); ++j) {
            auto i = P->getInteractions()[j];
            if (i->getI() == I) {
                return i;
            }
        }
    }
    return nullptr;
}

BaseInteraction* InteractionHandler::addInteraction(BaseInteractable* P, BaseInteractable* I, unsigned timeStamp)
{
    BaseSpecies* species = getDPMBase()->speciesHandler.getMixedObject(P->getIndSpecies(), I->getIndSpecies());
    BaseInteraction* C = species->getNewInteraction(P, I, timeStamp);
    addObject(C);
    return C;
}

void InteractionHandler::resetNewObjectsOMP() {
    #ifdef MERCURY_USE_OMP
    newObjects_.clear();
    newObjects_.resize(getDPMBase()->getNumberOfOMPThreads());
    #endif
}

void InteractionHandler::addNewObjectsOMP() {
    for (auto objects : newObjects_) {
        for (auto object : objects) {
            addObject(object);
        }
    }
}

BaseInteraction*
InteractionHandler::getInteraction(BaseInteractable* const P, BaseInteractable* const I, const unsigned timeStamp)
{
    const BaseSpecies* const species = getDPMBase()->speciesHandler.getMixedObject(P->getIndSpecies(),
                                                                                   I->getIndSpecies());
    
    BaseInteraction* C = getExistingInteraction(P, I);
    if (C == nullptr) {
        C = species->getNewInteraction(P, I, timeStamp);
        #ifdef MERCURY_USE_OMP
        if (omp_get_num_threads()>1) {
            newObjects_[omp_get_thread_num()].push_back(C);
            C->setHandler(this);
        } else {
            addObject(C);
        }
        #else
        addObject(C);
        #endif
    }
    
    //set timeStamp
    C->setTimeStamp(timeStamp);
    
    //set species of collision
    C->setSpecies(species);
    
    return C;
}

BaseInteraction* InteractionHandler::createEmptyInteraction() const
{
    //NOTE: assuming the interaction type of a species is the same for all species
    BaseSpecies* species = this->getDPMBase()->speciesHandler.getObject(0);
    BaseInteraction* emptyInteraction = species->getEmptyInteraction();
    
    return emptyInteraction;
}

void InteractionHandler::deleteEmptyInteraction(BaseInteraction* interaction) const
{
    BaseSpecies* species = this->getDPMBase()->speciesHandler.getObject(0);
    species->deleteEmptyInteraction(interaction);
}

void* InteractionHandler::createMPIInteractionDataArray(unsigned int numberOfInteractions) const
{
    BaseInteraction* emptyInteraction = this->createEmptyInteraction();
    //create a vector based on the first interaction. This is to avoid templating the whole class
    void* historyData = emptyInteraction->createMPIInteractionDataArray(numberOfInteractions);
    this->deleteEmptyInteraction(emptyInteraction);
    
    return historyData;
}

void InteractionHandler::deleteMPIInteractionDataArray(void* dataArray)
{
    BaseInteraction* emptyInteraction = this->createEmptyInteraction();
    emptyInteraction->deleteMPIInteractionDataArray(dataArray);
    this->deleteEmptyInteraction(emptyInteraction);
}

void InteractionHandler::getInteractionDetails(void* interactionData, unsigned int index, unsigned int& identificationP,
                                               unsigned int& identificationI, bool& isWallInteraction,
                                               unsigned& timeStamp)
{
    BaseInteraction* emptyInteraction = this->createEmptyInteraction();
    emptyInteraction->getInteractionDetails(interactionData, index, identificationP, identificationI, isWallInteraction,
                                            timeStamp);
    this->deleteEmptyInteraction(emptyInteraction);
}

BaseInteraction*
InteractionHandler::getInteraction(BaseInteractable* P, BaseInteractable* I, unsigned timeStamp, const Vec3D& normal)
{
    BaseInteraction* c = nullptr;
    Mdouble dNormalMax = 0.9; //minimal agreement required
    for (const auto& i : P->getInteractions())
    {
        if (i->getI() == I && i->getTimeStamp() != timeStamp)
        {
            const Mdouble dNormal = Vec3D::dot(normal, i->getNormal());
            if (dNormal > dNormalMax)
            {
                c = i;
                dNormalMax = dNormal;
            }
        }
    }
    
    BaseSpecies* species = getDPMBase()->speciesHandler.getMixedObject(P->getIndSpecies(), I->getIndSpecies());
    
    if (c)
    {
        c->setTimeStamp(timeStamp);
    }
    else
    {
        c = species->getNewInteraction(P, I, timeStamp);
        c->setSpecies(species);
        addObject(c);
        //logger(INFO,"new interaction t=%: i=%, n=%",timeStamp,getNumberOfObjects(),normal);
    }
    
    return c;
}

void InteractionHandler::removeObjectKeepingPeriodics(const unsigned int id)
{
    BaseInteraction* iMain = getObject(id);
    
    BaseParticle* P = dynamic_cast<BaseParticle*>(iMain->getP());
    BaseParticle* I = dynamic_cast<BaseParticle*>(iMain->getI());
    if (P != nullptr && I != nullptr) //check that both P and I are particles (not walls)
    {
        BaseInteractable* realI = I->getPeriodicFromParticle();
        if (realI != nullptr && !P->getPeriodicFromParticle()) //check that P is a real and I is a ghost particle
        {
            if (P->getIndex() < realI->getIndex())
            {
                BaseInteraction* iOther = getExistingInteraction(P, realI);
                //You have to also check for existing interactions of the particles in reverse order since the copySwitchPointer function can revert the order of the particles
                if (iOther == nullptr)
                {
                    iOther = getExistingInteraction(realI, P);
                }
                if (iOther != nullptr) //if the interaction existed before the ghost particles were created
                {
                    //Here we decide which of the two interactions should be kept:
                    //the interaction between the two real particles (iMain), or
                    //the interaction between the real and a ghost particle (iOther).
                    //It picks the one for which a collision has happened,
                    //i.e. the one with the newer timeStamp.
                    if (iOther->getTimeStamp() <
                        iMain->getTimeStamp()) //if the interaction has been active during the last computeForce routine, make this the new (real) interaction.
                    
                    {
                        iMain->setI(realI);
                        removeObject(iOther->getIndex());
                        return;
                    }
                    else  //if the interaction has not been active during the last computeForce routine
                    {
                        BaseHandler<BaseInteraction>::removeObject(id);
                        return;
                    }
                }
                else //if the interaction has been created during the last computeForce routine, make this a new (real) interaction. 
                {
                    iMain->setI(realI);
                    return;
                }
            }
        }
    }
    //this point is reached if either P or I are walls, or P and I are both ghost particles; in these cases, the interaction gets deleted
    BaseHandler<BaseInteraction>::removeObject(id);
}

void InteractionHandler::eraseOldInteractions(unsigned currentNTime)
{
    //std::cout<<"void InteractionHandler::eraseOldInteractions(Mdouble lastTimeStep)"<<std::endl;
    //std::cout<<"Current interactions="<<getNumberOfObjects()<<std::endl;
    //Remove_if reconstructs the vector with only elements passing the check_spring_time function
    //Erase removes the end of the vector
    for (unsigned int id = 0; id < getNumberOfObjects(); id++)
    {
        if (getObject(id)->getTimeStamp() <= currentNTime)
        {
            getObject(id)->actionsOnErase();
            removeObject(id);
            --id;
        }
    }
}

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

Mdouble InteractionHandler::getMeanOverlap() const
{
    Mdouble sum = 0;
    Mdouble n = 0;
    for (BaseInteraction* const p : objects_)
    {
        if (p->getOverlap() > 0)
        {
            sum += p->getOverlap();
            n++;
        }
    }
    return sum / n;
}

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

void InteractionHandler::write(std::ostream& os) const
{
#ifdef MERCURY_USE_MPI
    //note: this function only prints the particles on processor 0
    //The rest of the particles are processed in the restart file write function
    MPIContainer& communicator = MPIContainer::Instance();
    unsigned int totalNumberOfInteractions = getNumberOfObjects();
    os << "Interactions " << totalNumberOfInteractions << '\n';
    //logger(INFO,"% interactions",getNumberOfObjects());
    for (BaseInteraction* it : *this)
    {
    os << (*it) << '\n';
    }
#else
    os << "Interactions " << getNumberOfObjects() << '\n';
    for (BaseInteraction* i : objects_)
        os << (*i) << '\n';
#endif
}

void InteractionHandler::read(std::istream& is)
{
    clear();
    unsigned int N;
    std::string dummy;
    is >> dummy;
    std::stringstream line;
    helpers::getLineFromStringStream(is, line);
    line >> N;
    if (N > 1e5)
    {
        logger(INFO, "Reading % %", N, dummy);
    }
    logger(VERBOSE, "In %::read(is): reading in % objects.", getName(), N);
    setStorageCapacity(N);
    // set map
    particleById.clear();
    for (BaseParticle* p: getDPMBase()->particleHandler)
    {
        particleById[p->getId()] = p;
    }
    wallById.clear();
    for (BaseWall* w: getDPMBase()->wallHandler)
    {
        wallById[w->getId()] = w;
    }
    for (unsigned int i = 0; i < N; i++)
    {
        readAndAddObject(is);
    }
    particleById.clear();
    wallById.clear();
}


void InteractionHandler::readAndAddObject(std::istream& is)
{
    std::string type, dummy, idType;
    unsigned int id0, id1;
    Mdouble doubleTimeStamp;
    unsigned timeStamp;
    
    is >> type;
    logger(DEBUG, "InteractionHandler::readAndAddObject(is): reading type %.", type);
    Mdouble timeStampDouble;
    is >> idType >> id0 >> id1 >> dummy >> timeStampDouble;
    timeStamp = timeStampDouble; //in order to read old restart files

    BaseInteraction* C;

#ifdef MERCURY_USE_MPI
    if (idType.compare("particleIds") == 0)
    {
        std::vector<BaseParticle*> list0 = getDPMBase()->particleHandler.getObjectsById(id0);
        std::vector<BaseParticle*> list1 = getDPMBase()->particleHandler.getObjectsById(id1);
        for (int p0 = 0; p0 < list0.size(); p0++)
        {
            for (int p1  = 0; p1 < list1.size(); p1++)
            {
                C = getInteraction(list0[p0],list1[p1], timeStamp);
                if (C != nullptr)
                {
                   is >> *C;
                }
            }
        }

    }
    else
    {
        std::vector<BaseParticle*> list0 = getDPMBase()->particleHandler.getObjectsById(id0);
        for (int p0 = 0; p0 < list0.size(); p0++)
        {
            C = getInteraction(list0[p0],getDPMBase()->wallHandler.getObjectById(id1), timeStamp);
            if (C != nullptr)
            {
                is >> *C;
            }
        }
    }
#else
    //this requires particleById/wallById to be set, which is done in the read() function
    if (idType == "particleIds")
        C = getInteraction(particleById[id0],particleById[id1],timeStamp);
    else
        C = getInteraction(particleById[id0],wallById[id1], timeStamp);
    is >> (*C);
#endif
    is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
}

void InteractionHandler::setWriteVTK(FileType fileType)
{
    writeVTK_ = fileType;
}

FileType InteractionHandler::getWriteVTK() const
{
    return writeVTK_;
}

double InteractionHandler::getLiquidBridgeVolume() const {
    double liquidVolume = 0;
    for (auto i : objects_) {
        auto j = dynamic_cast<LiquidMigrationWilletInteraction*>(i);
        if (j and !static_cast<BaseParticle*>(j->getP())->isMPIParticle()) liquidVolume += j->getLiquidBridgeVolume();
    }
    return getMPISum(liquidVolume);
};