MercuryBase.cc

Go to the documentation of this file.

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

#include <limits>
#include <string.h>

#include "MercuryBase.h"
#include "Math/Helpers.h"
#include "Particles/BaseParticle.h"
#include "Walls/BaseWall.h"

MercuryBase::MercuryBase()
{
    constructor();
    logger(DEBUG, "MercuryBase::MercuryBase() constructor finished");
}

MercuryBase::~MercuryBase()
{
    if (grid != nullptr)
    {
        delete grid;
        grid = nullptr;
    }
    logger(DEBUG, "MercuryBase::~MercuryBase() destructor finished.");
}

MercuryBase::MercuryBase(const MercuryBase& mercuryBase)
{
    grid = nullptr;
    gridNeedsUpdate_ = true;
    
    hGridMethod_ = mercuryBase.hGridMethod_;
    hGridDistribution_ = mercuryBase.hGridDistribution_;
    
    currentMaxRelativeDisplacement_ = mercuryBase.currentMaxRelativeDisplacement_;
    totalCurrentMaxRelativeDisplacement_ = mercuryBase.totalCurrentMaxRelativeDisplacement_;
    
    updateEachTimeStep_ = mercuryBase.updateEachTimeStep_;
    hGridMaxLevels_ = mercuryBase.hGridMaxLevels_;
    hGridCellOverSizeRatio_ = mercuryBase.hGridCellOverSizeRatio_;
    
    logger(DEBUG, "HGRID_base(HGrid_base& other) constructor finished.");
}

void MercuryBase::constructor()
{
    grid = nullptr;
    gridNeedsUpdate_ = true;
    hGridMaxLevels_ = 3;
    hGridCellOverSizeRatio_ = 1.0;
    updateEachTimeStep_ = true;
    hGridDistribution_ = EXPONENTIAL;
    hGridMethod_ = TOPDOWN;
}

void MercuryBase::hGridActionsBeforeTimeLoop()
{
}

void MercuryBase::read(std::istream& is)
{
    DPMBase::read(is);
    
    std::stringstream line(std::stringstream::in | std::stringstream::out);
    helpers::getLineFromStringStream(is, line);
    
    std::string dummy;
    
    line >> dummy >> dummy
            >> dummy >> hGridMaxLevels_
            >> dummy >> hGridCellOverSizeRatio_;
}

void MercuryBase::write(std::ostream& os, bool writeAllParticles) const
{
    DPMBase::write(os, writeAllParticles);
    if (grid != nullptr)
    {
        os << "NUM_BUCKETS " << grid->getNumberOfBuckets() << " "
                << "hGridMaxLevels " << getHGridMaxLevels() << " "
                << "cellOverSizeRatio " << getHGridCellOverSizeRatio() << std::endl;
    }
    else
    {
        os << "NUM_BUCKETS " << 0.0 << " "
                << "hGridMaxLevels " << getHGridMaxLevels() << " "
                << "cellOverSizeRatio " << getHGridCellOverSizeRatio() << std::endl;
    }
}

Mdouble MercuryBase::getHGridCurrentMaxRelativeDisplacement() const
{
    return currentMaxRelativeDisplacement_;
}

Mdouble MercuryBase::getHGridTotalCurrentMaxRelativeDisplacement() const
{
    return totalCurrentMaxRelativeDisplacement_;
}

void MercuryBase::setHGridUpdateEachTimeStep(bool updateEachTimeStep)
{
    updateEachTimeStep_ = updateEachTimeStep;
}

bool MercuryBase::getHGridUpdateEachTimeStep() const
{
    return updateEachTimeStep_;
}

void
MercuryBase::hGridRebuild()
{

    std::vector<Mdouble> cellSizes;

    const Mdouble minParticleInteractionRadius = getHGridTargetMinInteractionRadius();
    const Mdouble maxParticleInteractionRadius = getHGridTargetMaxInteractionRadius();
    if(minParticleInteractionRadius == 0.0 || minParticleInteractionRadius == maxParticleInteractionRadius)
    {
        //this case is executed if the particleHandler is empty (minParticleInteractionRadius == 0)
        //or if the particle distribution is monodispersed. 
        //nextafter(d,std::numeric_limits<Mdouble>::max()) chooses the smallest 
        // Mdouble that is bigger than d.
        const Mdouble maxCellSize = nextafter(2.0 * maxParticleInteractionRadius * getHGridCellOverSizeRatio(), std::numeric_limits<Mdouble>::max());
        cellSizes.push_back(maxCellSize);
        if(getHGridMaxLevels() != 1)
        {
            logger(DEBUG, "While rebuilding the hgrid: the number of levels was "
                    "set to one, as the particle distribution is monodispersed");
        }
    }
    else
    {
        switch(getHGridDistribution())
        {
        case LINEAR:
        {
            const Mdouble minCellSize = nextafter(2.0 * minParticleInteractionRadius * getHGridCellOverSizeRatio(), 0.0);
            const Mdouble maxCellSize = nextafter(2.0 * maxParticleInteractionRadius * getHGridCellOverSizeRatio(), std::numeric_limits<Mdouble>::max());
            //std::cout << "HGrid: using a linear cell size distribution from " << minCellSize << " to " << maxCellSize << " over " << getHGridMaxLevels() << " levels" << std::endl;
            for(unsigned int i = 0; i + 1 < getHGridMaxLevels(); i++)
            {
                cellSizes.push_back(minCellSize + (maxCellSize - minCellSize) 
                * (static_cast<Mdouble>(i + 1)) / getHGridMaxLevels());
            }
            //The last cell is added separately because in some cases accuracy was lost when calculating it.
            cellSizes.push_back(maxCellSize);
            break;
        }
        case EXPONENTIAL:
        {
            const Mdouble minCellSize = nextafter(2.0 * minParticleInteractionRadius * getHGridCellOverSizeRatio(), 0.0);
            const Mdouble maxCellSize = nextafter(2.0 * maxParticleInteractionRadius * getHGridCellOverSizeRatio(), std::numeric_limits<Mdouble>::max());
            //std::cout << "HGrid: using an exponential cell size distribution from " << minCellSize << " to " << maxCellSize << " over " << getHGridMaxLevels() << " levels" << std::endl;
            for(unsigned int i = 0; i + 1 < getHGridMaxLevels(); i++)
            {
                cellSizes.push_back(minCellSize 
                * std::pow(maxCellSize / minCellSize, static_cast<Mdouble>(i + 1) 
                / getHGridMaxLevels()));
            }
            //The last cell is added separately because in some cases accuracy was lost when calculating it.
            cellSizes.push_back(maxCellSize);
            break;
        }
        case USER:
        {
            for(unsigned int i = 0; i < getHGridMaxLevels(); i++)
            {
                cellSizes.push_back(userHGridCellSize(i));
            }
            break;
        }
        case OLDHGRID:
        {
            const Mdouble minCellSize = nextafter(2.0 * minParticleInteractionRadius * getHGridCellOverSizeRatio(), 0.0);

            //std::cout<<"HGrid: using the old HGrid cell size distribution starting from " <<minCellSize<<std::endl;
            for(unsigned int i = 0; i < getHGridMaxLevels(); i++)
            {
                cellSizes.push_back(minCellSize * std::pow(2, i));
            }
            break;
        }
        }
    }

    if(grid != nullptr)
    {
        delete grid;
    }

    grid = new HGrid(getHGridTargetNumberOfBuckets(), getHGridCellOverSizeRatio(), cellSizes);

    for(BaseParticle* const p: particleHandler)
    {
        hGridInsertParticle(p);
        //\todo{This is really ugly fix to force the particle to update}
        (p)->setHGridX(9999);
        hGridUpdateParticle(p);
    }
    gridNeedsUpdate_ = false;
}

void MercuryBase::hGridInsertParticle(BaseParticle *obj)
{
    if (grid != nullptr)
    {
        grid->insertParticleToHgrid(obj);
    }
}

#ifdef CONTACT_LIST_HGRID
void MercuryBase::broadPhase(BaseParticle *i UNUSED)
{
}
#else

void MercuryBase::broadPhase(BaseParticle *i)
{
    hGridFindOneSidedContacts(i);
}
#endif

void MercuryBase::hGridActionsBeforeTimeStep()
{
    static int stepsBeforeUpdate = 0;
    if (hGridNeedsRebuilding())
    {
        logger(VERBOSE, "HGrid needs rebuilding for void HGRID_base::hGridActionsBeforeTimeStep()");
        hGridRebuild();
        totalCurrentMaxRelativeDisplacement_ = 0;
        stepsBeforeUpdate = 0;
    }
    else
    {
#ifndef CONTACT_LIST_HGRID
        getHGrid()->clearBucketIsChecked();
#endif
        if (getHGridUpdateEachTimeStep() || getHGridTotalCurrentMaxRelativeDisplacement() >= getHGridCellOverSizeRatio() - 1)
        {
#ifndef CONTACT_LIST_HGRID
            getHGrid()->clearFirstBaseParticleInBucket();
#endif
            totalCurrentMaxRelativeDisplacement_ = 0;
            for (BaseParticle* const p : particleHandler)
            {
                hGridUpdateParticle(p);
            }
            stepsBeforeUpdate = 0;
        }
        else
        {
            stepsBeforeUpdate++;
        }
    }

}

void MercuryBase::hGridUpdateMove(BaseParticle* iP, Mdouble move)
{
    const Mdouble currentRelativeDisplacement = move / (getHGrid()->getCellSize(iP->getHGridLevel()));
    if (currentRelativeDisplacement > currentMaxRelativeDisplacement_)
    {
        currentMaxRelativeDisplacement_ = currentRelativeDisplacement;
    }
}

void MercuryBase::hGridActionsBeforeIntegration()
{
    currentMaxRelativeDisplacement_ = 0.0;
}

void MercuryBase::hGridActionsAfterIntegration()
{
    totalCurrentMaxRelativeDisplacement_ += 2.0 * currentMaxRelativeDisplacement_;
}

HGrid* MercuryBase::getHGrid()
{
    return grid;
}

const HGrid* MercuryBase::getHGrid() const
{
    return grid;
}

Mdouble MercuryBase::userHGridCellSize(unsigned int level)
{
    logger(WARN, "In Mdouble MercuryBase::userHGridCellSize(unsigned int level) with level= %", level);
    logger(WARN, "If you want to use user defined HGrid cell sizes, this function should be redefined");
    return 0.0;
}

bool MercuryBase::readNextArgument(int& i, int argc, char *argv[])
{
    if (!strcmp(argv[i], "-hGridMaxLevels"))
    {
        setHGridMaxLevels(static_cast<unsigned int>(atoi(argv[i + 1])));
    }
    else if (!strcmp(argv[i], "-cellOverSizeRatio"))
    {
        setHGridCellOverSizeRatio(atof(argv[i + 1]));
    }
    else
    {
        return DPMBase::readNextArgument(i, argc, argv); //if argv[i] is not found, check the commands in MD
    }
    return true; //returns true if argv[i] is found
}

HGridMethod MercuryBase::getHGridMethod() const
{
    return hGridMethod_;
}

void MercuryBase::setHGridMethod(HGridMethod hGridMethod)
{
    hGridMethod_ = hGridMethod;
}

HGridDistribution MercuryBase::getHGridDistribution() const
{
    return hGridDistribution_;
}

void MercuryBase::setHGridDistribution(HGridDistribution hGridDistribution)
{
    if (hGridDistribution_ != hGridDistribution)
    {
        gridNeedsUpdate_ = true;
        hGridDistribution_ = hGridDistribution;
    }
}

Mdouble MercuryBase::getHGridCellOverSizeRatio() const
{
    return hGridCellOverSizeRatio_;
}

void MercuryBase::setHGridCellOverSizeRatio(Mdouble hGridCellOverSizeRatio)
{
    //If the hGridCellOverSizeRatio changes significantly, assign the given parameter.
    if (mathsFunc::isEqual(hGridCellOverSizeRatio_, hGridCellOverSizeRatio, 1e-10))
    {
        gridNeedsUpdate_ = true;
        hGridCellOverSizeRatio_ = hGridCellOverSizeRatio;
    }
}

void MercuryBase::setHGridMaxLevels(unsigned int hGridMaxLevels)
{
    if (hGridMaxLevels_ != hGridMaxLevels)
    {
        gridNeedsUpdate_ = true;
        hGridMaxLevels_ = hGridMaxLevels;
    }
}

unsigned int MercuryBase::getHGridMaxLevels() const
{
    return hGridMaxLevels_;
}

bool MercuryBase::hGridNeedsRebuilding()
{
    if (grid == nullptr)
    {
        logger(VERBOSE, "HGrid needs updating, because there is no grid.");
        return true;
    }
    else if (gridNeedsUpdate_)
    {
        logger(VERBOSE, "HGrid needs updating, because some of its initialisation parameters have changed.");
        return true;
    }
    else if (grid->getNeedsRebuilding())
    {
        logger(VERBOSE, "HGrid needs updating, because said so by the grid itself");
        return true;
    }
    else if (getHGrid()->getNumberOfBuckets() > 10 * getHGridTargetNumberOfBuckets() || 10 * getHGrid()->getNumberOfBuckets() < getHGridTargetNumberOfBuckets())
    {
        logger(VERBOSE, "HGrid needs updating, because of number of buckets, current = %, target = %.", grid->getNumberOfBuckets(), particleHandler.getNumberOfObjects());
        return true;
    }
    else if (particleHandler.getLargestParticle() != nullptr && 2.0 * particleHandler.getLargestParticle()->getInteractionRadius() > getHGrid()->getCellSizes().back() * grid->getCellOverSizeRatio())
    {
        logger(VERBOSE, "HGrid needs updating, because of maximum cell size, current = %, required = %.", grid->getCellSizes().back() * hGridCellOverSizeRatio_, particleHandler.getLargestParticle()->getInteractionRadius());
        return true;
    }
    else
    {
        //std::cout<<"HGrid does not need updating, because of number of buckets, current="<<grid->NUM_BUCKETS<<" target="<<particleHandler.getNumberOfObjects()<<std::endl;
        //std::cout<<"HGrid does not need updating, because of maximum cell size, current="<<grid->cellSizes_.back()*grid->cellOverSizeRatio_<<" required="<<2.0*particleHandler.getLargestParticle()->getInteractionRadius()<<std::endl;
        return false;
    }
}

unsigned int MercuryBase::getHGridTargetNumberOfBuckets() const
{
    const unsigned int nParticles = particleHandler.getNumberOfObjects();
    if (nParticles > 10)
    {
        return nParticles;
    }
    else
    {
        return 10;
    }
}

Mdouble MercuryBase::getHGridTargetMinInteractionRadius() const
{
    if (particleHandler.getNumberOfObjects() == 0)
    {
        return 0.0;
    }
    else
    {
        return particleHandler.getSmallestParticle()->getInteractionRadius() * getHGridCellOverSizeRatio();
    }
}

Mdouble MercuryBase::getHGridTargetMaxInteractionRadius() const
{
    if (particleHandler.getNumberOfObjects() == 0)
    {
        return 0.0;
    }
    else
    {
        return particleHandler.getLargestParticle()->getInteractionRadius() * getHGridCellOverSizeRatio();
    }
}
bool MercuryBase::checkParticleForInteraction(const BaseParticle& p)
{
#ifdef MERCURY_USE_MPI
    int localInteraction = checkParticleForInteractionLocal(p);
    
    MPIContainer& communicator = MPIContainer::Instance();
    int numberOfProcessors = communicator.getNumberOfProcessors();
    
    //The root gathers all values and computes the global value
    int *interactionList = nullptr;
    if (communicator.getProcessorID() == 0)
    {
        interactionList = new int [numberOfProcessors];
    }

    //Gather all local values
    communicator.gather(localInteraction,interactionList);

    //Compute the global value
    int globalInteraction = 1;
    if (communicator.getProcessorID() == 0)
    {
        for (int i = 0; i<numberOfProcessors; i++)
        {
            if (interactionList[i] == 0)
            {
                globalInteraction = 0;
                break;
            }
        }
    }
    //The root now tells the other processors what the global value for the interaction is
    communicator.broadcast(globalInteraction);

    //Convert the result back to bool
    bool interaction = globalInteraction;
#else
    bool interaction = checkParticleForInteractionLocalPeriodic(p);
#endif   
    return interaction;
}

bool MercuryBase::checkParticleForInteractionLocal(const BaseParticle& p)
{
    Mdouble distance;
    Vec3D normal;

    //Check if it has no collision with walls
    for (BaseWall* w :wallHandler)
    {
        if (w->getDistanceAndNormal(p, distance, normal))
        {
            logger(VERBOSE, "Collision with wall %.", *w);
            return false;
        }
        else
        {
            logger(VERBOSE, "No collision with wall %.", *w);
        }
    }
    
    //Check if it has no collision with other particles
    if (hGridHasParticleContacts(&p))
    {
        logger(VERBOSE, "Collision with particle.");
        return false;
    }
    else
    {
        logger(VERBOSE, "No collision with particles.");
    }
    return true;
}


void MercuryBase::hGridInfo() const
{
    logger(INFO,"hGrid parameters:" );
    logger(INFO, "  method %", hGridMethod_);
    logger(INFO, "  distribution %", hGridDistribution_);
    logger(VERBOSE, "Mdouble currentMaxRelativeDisplacement_=%", currentMaxRelativeDisplacement_);
    logger(VERBOSE, "Mdouble totalCurrentMaxRelativeDisplacement_=%", totalCurrentMaxRelativeDisplacement_);
    logger(INFO,"  needsUpdate %", gridNeedsUpdate_);
    logger(INFO,"  updateEachTimeStep %", updateEachTimeStep_);
    logger(INFO, "  maxLevels %", hGridMaxLevels_);
    logger(INFO, "  cellOverSizeRatio %", hGridCellOverSizeRatio_);
    std::cout << "hGrid parameters:" << std::endl;
    std::cout << "  method " << hGridMethod_ << std::endl; 
    std::cout << "  distribution " << hGridDistribution_ << std::endl;
    if (grid != nullptr) 
    {
        grid->info();
    }
    else
    {
        logger(INFO,"  grid does not exist yet");
    }
}