InteractionHandler.cc

Go to the documentation of this file.

//Copyright (c) 2013-2014, The MercuryDPM Developers Team. All rights reserved.
//For the list of developers, see <http://www.MercuryDPM.org/Team>.
//
//Redistribution and use in source and binary forms, with or without
//modification, are permitted provided that the following conditions are met:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name MercuryDPM nor the
//    names of its contributors may be used to endorse or promote products
//    derived from this software without specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
//ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
//WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
//DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
//LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
//ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
//(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
//SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include <iostream>
#include <iomanip>
#include <string>
#include <stdlib.h>
#include "Particles/BaseParticle.h"
#include "InteractionHandler.h"
#include "SpeciesHandler.h"
#include "Species/BaseSpecies.h"
#include "Interactions/BaseInteraction.h"
#include "DPMBase.h"

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)
{
    //Puts the particle in the Particle list
    BaseHandler<BaseInteraction>::addObject(I);
    //set the particleHandler pointer
    I->setHandler(this);
}

BaseInteraction* InteractionHandler::getExistingInteraction(BaseInteractable* P, BaseInteractable* 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 (BaseInteraction* i :  P->getInteractions())
    {
        if (i->getI() == I)
        {
            return i;
        }
    }
    return nullptr;
}

BaseInteraction* InteractionHandler::getInteraction(BaseInteractable* P, BaseInteractable* I, Mdouble timeStamp)
{
    BaseSpecies* species = getDPMBase()->speciesHandler.getMixedObject(P->getIndSpecies(), I->getIndSpecies());
    
    //std::cout << "Trying to reconnect to BaseInteraction between P=" << P->getId() << " and " << I->getId() << std::endl;
    BaseInteraction* C = getExistingInteraction(P, I);
    if (C == nullptr)
    {
        C = species->getNewInteraction(P, I, timeStamp);
        addObject(C);
        //std::cout << "Creating new interaction with index=" << getLastObject()->getIndex() << " id=" << getLastObject()->getId() << std::endl;
    }
    
    //set timeStamp
    C->setTimeStamp(timeStamp);
    
    //set species of collision
    C->setSpecies(species);
    
    return C;
}

BaseInteraction*
InteractionHandler::getInteraction(BaseInteractable* P, BaseInteractable* I, Mdouble 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(Mdouble lastTimeStep)
{
    //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() <= lastTimeStep)
        {
            getObject(id)->actionsOnErase();
            removeObject(id);
            --id;
        }
    }
}

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

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

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

void InteractionHandler::write(std::ostream& os) const
{
    os << "Interactions " << getNumberOfObjects() << std::endl;
    for (BaseInteraction* i : objects_)
        os << (*i) << std::endl;
}

void InteractionHandler::readObject(std::istream& is)
{
    std::string type, dummy, idType;
    unsigned int id0, id1;
    Mdouble timeStamp;
    
    std::stringstream line(std::stringstream::in | std::stringstream::out);
    helpers::getLineFromStringStream(is, line);
    
    line >> type >> idType >> id0 >> id1 >> dummy >> timeStamp;
    logger(VERBOSE, "InteractionHandler::readObject(is): reading type % % %", type, id0, id1);
    BaseInteraction* C;
    if (idType.compare("particleIds") == 0) {
//        C =  getDPMBase()->particleHandler.getObjectById(id1)->getInteractionWith(getDPMBase()->particleHandler.getObjectById(id0), timeStamp, this)[0];
        C = getInteraction(getDPMBase()->particleHandler.getObjectById(id0),
                           getDPMBase()->particleHandler.getObjectById(id1), timeStamp);
    }
    else {
//        C =  getDPMBase()->wallHandler.getObjectById(id1)->getInteractionWith(getDPMBase()->particleHandler.getObjectById(id0), timeStamp, this)[0];
        C = getInteraction(getDPMBase()->particleHandler.getObjectById(id0),
                           getDPMBase()->wallHandler.getObjectById(id1), timeStamp);
    }
    line >> (*C);
}

void InteractionHandler::writeVTK() const
{
    if (getWriteVTK()==FileType::NO_FILE)
        return;

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

    std::fstream file;
    file.open(fileName.c_str(), std::ios::out);
    if (file.fail())
    {
        logger(WARN, "Error in writeToFile: file % could not be opened", fileName);
    }

    file << "<?xml version=\"1.0\"?>\n\n";
    file << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
    file << "<UnstructuredGrid>\n";
    file << "<Piece NumberOfPoints=\"" << this->getNumberOfObjects() << "\" NumberOfCells=\"" << 0 << "\">\n";
    file << "<Points>\n";
    file << "  <DataArray type=\"Float32\" Name=\"Position\" NumberOfComponents=\"3\" format=\"ascii\">\n";
    for (const auto &p: *this)
    {
        file << '\t' << p->getContactPoint() << '\n';
    }
    file << "  </DataArray>\n";
    file << "</Points>\n";
    file << "<PointData  Vectors=\"vector\">\n";
    file << "  <DataArray type=\"Float32\" Name=\"Normal\" NumberOfComponents=\"3\" format=\"ascii\">\n";
    // Add velocity
    for (const auto &p: *this)
    {
        file << '\t' << p->getNormal() << '\n';
    }
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"Float32\" Name=\"Overlap\" format=\"ascii\">\n";
    // Add radius
    for (const auto &p: *this)
    {
        file << '\t' << p->getOverlap() << '\n';
    }
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"Float32\" Name=\"ContactRadius\" format=\"ascii\">\n";
    // Add radius
    for (const auto &p: *this)
    {
        file << '\t' << p->getContactRadius() << '\n';
    }
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"Float32\" Name=\"Force\" format=\"ascii\">\n";
    // Add species type
    for (const auto &p: *this)
    {
        file << '\t' << p->getForce() << '\n';
    }
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"Float32\" Name=\"TangentialOverlap\" format=\"ascii\">\n";
    // Add species type
    for (const auto &p: *this)
    {
        file << '\t' << p->getTangentialOverlap() << '\n';
    }
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"Float32\" Name=\"Torque\" format=\"ascii\">\n";
    // Add species type
    for (const auto &p: *this)
    {
        file << '\t' << p->getTorque() << '\n';
    }
    file << "  </DataArray>\n";
    //check if this type of Interaction has extra fields
    if (getLastObject()!= nullptr)
    {
        for (unsigned i = 0; i<getLastObject()->getNumberOfFieldsVTK(); i++)
        {
            file << "  <DataArray type=\"" << getLastObject()->getTypeVTK(i) << "\" Name=\"" << getLastObject()->getNameVTK(i) << "\" format=\"ascii\">\n";
            // Add species type
            for (const auto &p: *this)
            {
                for (auto f : p->getFieldVTK(i))
                    file << '\t' << f << '\n';
            }
            file << "  </DataArray>\n";
        }
    }
    file << "</PointData>\n";
    file << "<Cells>\n";
    file << "  <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n";
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n";
    file << "  </DataArray>\n";
    file << "  <DataArray type=\"UInt8\"  Name=\"types\" format=\"ascii\">\n";
    file << "  </DataArray>\n";
    file << "</Cells>\n";
    file << "</Piece>\n";
    file << "</UnstructuredGrid>\n";
    file << "</VTKFile>\n";
    file.close();
}

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

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