Mercury3D.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 "Mercury3D.h"
#include "Particles/BaseParticle.h"

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

Mercury3D::Mercury3D(const Mercury3D& other)
        : DPMBase(other), MercuryBase(other)
{
    logger(DEBUG, "Mercury3D::Mercury3D(Mercury3D& other) copy constructor finished.");
}

Mercury3D::Mercury3D(const DPMBase& other)
        : DPMBase(other), MercuryBase()
{
    constructor();
    logger(DEBUG, "Mercury3D::Mercury3D(DPMBase& other) constructor finished");
}

void Mercury3D::constructor()
{
    setParticleDimensions(3);
    setSystemDimensions(3);
}

void Mercury3D::hGridFindContactsWithinTargetCell(int x, int y, int z, unsigned int l)
{
    HGrid* const hgrid = getHGrid();
    const unsigned int bucket = hgrid->computeHashBucketIndex(x, y, z, l);
    
    //Check if this function is already applied to this bucket
    if (hgrid->getBucketIsChecked(bucket))
    {
        return;
    }
    
    BaseParticle* p1 = hgrid->getFirstBaseParticleInBucket(bucket);
    while (p1 != nullptr)
    {
        BaseParticle* p2 = p1->getHGridNextObject();
        while (p2 != nullptr)
        {
            //Check if the BaseParticle* p1 and BaseParticle* p2 are really in the same cell (i.e. no hashing error has occurred)
            if ((p1->getHGridX() == p2->getHGridX()) && (p1->getHGridY() == p2->getHGridY()) && (p1->getHGridZ() == p2->getHGridZ()) && (p1->getHGridLevel() == p2->getHGridLevel()))
            {
                computeInternalForces(p1, p2);
            }
            p2 = p2->getHGridNextObject();
        }
        p1 = p1->getHGridNextObject();
    }
    hgrid->setBucketIsChecked(bucket);
}

void Mercury3D::hGridFindContactsWithTargetCell(int x, int y, int z, unsigned int l, BaseParticle *obj)
{
    
    //Check if the object is not in the same cell as being checked, CheckCell_current should handle these cases.
    if ((obj->getHGridX() == x) && (obj->getHGridY() == y) && (obj->getHGridZ() == z) && (obj->getHGridLevel() == l))
    {
        return;
    }

    HGrid* const hgrid = getHGrid();

    // Calculate the bucket
    const unsigned int bucket = hgrid->computeHashBucketIndex(x, y, z, l);

    // Loop through all objects in the bucket to find nearby objects
    BaseParticle* p = hgrid->getFirstBaseParticleInBucket(bucket);
    while (p != nullptr)
    {
        //This is the most-expensive function in most codes (the two checks in this function slows down granular jet by 15%). It is neccesary, for example: If two cells hash to the same bucket and a particle in one of these cells check for collisions with the other cell. Then due to the hashing collision it also gets all particles in it's own cell and thus generating false collisions.
        //Check if the BaseParticle *p really is in the target cell (i.e. no hashing error has occurred)
        if ((p->getHGridX() == x) && (p->getHGridY() == y) && (p->getHGridZ() == z) && (p->getHGridLevel() == l))
        {
            computeInternalForces(obj, p);
        }
        p = p->getHGridNextObject();
    }
}

void Mercury3D::hGridFindOneSidedContacts(BaseParticle* obj)
{
    HGrid* const hgrid = getHGrid();
    const unsigned int startLevel = obj->getHGridLevel();

    switch (getHGridMethod())
    {
        case BOTTOMUP:
            {
            int occupiedLevelsMask = hgrid->getOccupiedLevelsMask() >> obj->getHGridLevel();
            for (unsigned int level = startLevel; level < hgrid->getNumberOfLevels(); occupiedLevelsMask >>= 1, level++)
            {
                // If no objects in rest of grid, stop now
                if (occupiedLevelsMask == 0)
                {
                    break;
                }
                
                // If no objects at this level, go on to the next level
                if ((occupiedLevelsMask & 1) == 0)
                {
                    continue;
                }
                
                if (level == startLevel)
                {
                    const int x = obj->getHGridX();
                    const int y = obj->getHGridY();
                    const int z = obj->getHGridZ();
                    
                    hGridFindContactsWithinTargetCell(x, y, z, level);
                    hGridFindContactsWithTargetCell(x + 1, y - 1, z, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y, z, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y + 1, z, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y - 1, z + 1, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y, z + 1, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y + 1, z + 1, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y - 1, z - 1, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y, z - 1, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y + 1, z - 1, level, obj);
                    hGridFindContactsWithTargetCell(x, y + 1, z, level, obj);
                    hGridFindContactsWithTargetCell(x, y, z - 1, level, obj);
                    hGridFindContactsWithTargetCell(x, y + 1, z - 1, level, obj);
                    hGridFindContactsWithTargetCell(x, y + 1, z + 1, level, obj);
                }
                else
                {
                    const Mdouble inv_size = hgrid->getInvCellSize(level);
                    const int xs = static_cast<int>(std::floor((obj->getPosition().X - obj->getInteractionRadius()) * inv_size - 0.5));
                    const int xe = static_cast<int>(std::floor((obj->getPosition().X + obj->getInteractionRadius()) * inv_size + 0.5));
                    const int ys = static_cast<int>(std::floor((obj->getPosition().Y - obj->getInteractionRadius()) * inv_size - 0.5));
                    const int ye = static_cast<int>(std::floor((obj->getPosition().Y + obj->getInteractionRadius()) * inv_size + 0.5));
                    const int zs = static_cast<int>(std::floor((obj->getPosition().Z - obj->getInteractionRadius()) * inv_size - 0.5));
                    const int ze = static_cast<int>(std::floor((obj->getPosition().Z + obj->getInteractionRadius()) * inv_size + 0.5));
                    for (int x = xs; x <= xe; ++x)
                    {
                        for (int y = ys; y <= ye; ++y)
                        {
                            for (int z = zs; z <= ze; ++z)
                            {
                                hGridFindContactsWithTargetCell(x, y, z, level, obj);
                            }
                        }
                    }
                }
            }
            break;
        }
        case TOPDOWN:
            {
            int occupiedLevelsMask = hgrid->getOccupiedLevelsMask();
            for (unsigned int level = 0; level <= startLevel; occupiedLevelsMask >>= 1, level++)
            {
                // If no objects in rest of grid, stop now
                if (occupiedLevelsMask == 0)
                {
                    break;
                }
                
                // If no objects at this level, go on to the next level
                if ((occupiedLevelsMask & 1) == 0)
                {
                    continue;
                }
                
                if (level == startLevel)
                {
                    const int x = obj->getHGridX();
                    const int y = obj->getHGridY();
                    const int z = obj->getHGridZ();
                    
                    hGridFindContactsWithinTargetCell(x, y, z, level);
                    hGridFindContactsWithTargetCell(x + 1, y - 1, z, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y, z, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y + 1, z, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y - 1, z + 1, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y, z + 1, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y + 1, z + 1, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y - 1, z - 1, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y, z - 1, level, obj);
                    hGridFindContactsWithTargetCell(x + 1, y + 1, z - 1, level, obj);
                    hGridFindContactsWithTargetCell(x, y + 1, z, level, obj);
                    hGridFindContactsWithTargetCell(x, y, z - 1, level, obj);
                    hGridFindContactsWithTargetCell(x, y + 1, z - 1, level, obj);
                    hGridFindContactsWithTargetCell(x, y + 1, z + 1, level, obj);
                }
                else
                {
                    const Mdouble inv_size = getHGrid()->getInvCellSize(level);
                    const int xs = static_cast<int>(std::floor((obj->getPosition().X - obj->getInteractionRadius()) * inv_size - 0.5));
                    const int xe = static_cast<int>(std::floor((obj->getPosition().X + obj->getInteractionRadius()) * inv_size + 0.5));
                    const int ys = static_cast<int>(std::floor((obj->getPosition().Y - obj->getInteractionRadius()) * inv_size - 0.5));
                    const int ye = static_cast<int>(std::floor((obj->getPosition().Y + obj->getInteractionRadius()) * inv_size + 0.5));
                    const int zs = static_cast<int>(std::floor((obj->getPosition().Z - obj->getInteractionRadius()) * inv_size - 0.5));
                    const int ze = static_cast<int>(std::floor((obj->getPosition().Z + obj->getInteractionRadius()) * inv_size + 0.5));
                    for (int x = xs; x <= xe; ++x)
                    {
                        for (int y = ys; y <= ye; ++y)
                        {
                            for (int z = zs; z <= ze; ++z)
                            {
                                hGridFindContactsWithTargetCell(x, y, z, level, obj);
                            }
                        }
                    }
                }
            }
            break;
        }
    }
}

void Mercury3D::hGridUpdateParticle(BaseParticle* obj)
{
    HGrid* const hGrid=getHGrid();
    if (hGrid)
    {
        const unsigned int l = obj->getHGridLevel();
        const Mdouble inv_size = hGrid->getInvCellSize(l);

        int x = static_cast<int>(std::floor(obj->getPosition().X * inv_size));
        int y = static_cast<int>(std::floor(obj->getPosition().Y * inv_size));
        int z = static_cast<int>(std::floor(obj->getPosition().Z * inv_size));

#ifdef CONTACT_LIST_HGRID
        if(obj->getHGridX() != x || obj->getHGridY() != y || obj->getHGridZ() != z)
        {
            int bucket = hGrid->computeHashBucketIndex(x, y, z, l);

            //First the object has to be removed
            hGridRemoveParticle(obj);

            //Also remove all contact associated with it
            getPossibleContactList().remove_ParticlePosibleContacts(obj);

            //And now reinserted
            obj->setHGridNextObject(hGrid->getFirstBaseParticleInBucket(bucket));
            obj->setHGridPrevObject(nullptr);
            if(hGrid->getFirstBaseParticleInBucket(bucket))
            {
                hGrid->getFirstBaseParticleInBucket(bucket)->setHGridPrevObject(obj);
            }
            hGrid->setFirstBaseParticleInBucket(bucket,obj);

            obj->setHGridX(x);
            obj->setHGridY(y);
            obj->setHGridZ(z);
            InsertObjAgainstGrid(obj);
        }
#else
        const unsigned int bucket = hGrid->computeHashBucketIndex(x, y, z, l);

        obj->setHGridNextObject(hGrid->getFirstBaseParticleInBucket(bucket));
        obj->setHGridPrevObject(nullptr);
        if (hGrid->getFirstBaseParticleInBucket(bucket))
        {
            hGrid->getFirstBaseParticleInBucket(bucket)->setHGridPrevObject(obj);
        }

        hGrid->setFirstBaseParticleInBucket(bucket, obj);

        obj->setHGridX(x);
        obj->setHGridY(y);
        obj->setHGridZ(z);
#endif
    }
}

void Mercury3D::hGridRemoveParticle(BaseParticle* obj)
{
    HGrid* const hGrid=getHGrid();
    if (hGrid)
    {
        const unsigned int bucket = hGrid->computeHashBucketIndex(obj->getHGridX(), obj->getHGridY(), obj->getHGridZ(), obj->getHGridLevel());
        if (obj->getHGridPrevObject())
        {
            obj->getHGridPrevObject()->setHGridNextObject(obj->getHGridNextObject());
        }
        else
        {
            if (hGrid->getFirstBaseParticleInBucket(bucket) == obj)
            {
                hGrid->setFirstBaseParticleInBucket(bucket, obj->getHGridNextObject());
            }
        }
        
        if (obj->getHGridNextObject())
        {
            obj->getHGridNextObject()->setHGridPrevObject(obj->getHGridPrevObject());
        }
    }
}

bool Mercury3D::hGridHasContactsInTargetCell(int x, int y, int z, unsigned int l, const BaseParticle* obj) const
{
    // Loop through all objects in the bucket to find nearby objects
    const unsigned int bucket = getHGrid()->computeHashBucketIndex(x, y, z, l);
    
    const BaseParticle* p = getHGrid()->getFirstBaseParticleInBucket(bucket);
    while (p != nullptr)
    {
        if ((p->getHGridX() == x) && (p->getHGridY() == y) && (p->getHGridZ() == z) && (p->getHGridLevel() == l))
        {
            if (areInContact(obj, p))
            {
                return true;
            }
        }
        p = p->getHGridNextObject();
    }
    return false;
}

bool Mercury3D::hGridHasParticleContacts(const BaseParticle* obj)
{
    if (getHGrid() == nullptr || getHGrid()->getNeedsRebuilding())
    {
        logger(INFO, "HGrid needs rebuilding for \"bool Mercury3D::hGridHasParticleContacts(BaseParticle *obj)\"");
        hGridRebuild();
    }

    int occupiedLevelsMask = getHGrid()->getOccupiedLevelsMask();
    
    for (unsigned int level = 0; level < getHGrid()->getNumberOfLevels(); occupiedLevelsMask >>= 1, level++)
    {
        // If no objects in rest of grid, stop now
        if (occupiedLevelsMask == 0)
        {
            logger(VERBOSE, "Level % and higher levels are empty", level);
            break;
        }
        
        // If no objects at this level, go on to the next level
        if ((occupiedLevelsMask & 1) == 0)
        {
            logger(VERBOSE, "Level % is empty", level);
            continue;
        }

        const Mdouble inv_size = getHGrid()->getInvCellSize(level);
        const int xs = static_cast<int>(std::floor((obj->getPosition().X - obj->getInteractionRadius()) * inv_size - 0.5));
        const int xe = static_cast<int>(std::floor((obj->getPosition().X + obj->getInteractionRadius()) * inv_size + 0.5));
        const int ys = static_cast<int>(std::floor((obj->getPosition().Y - obj->getInteractionRadius()) * inv_size - 0.5));
        const int ye = static_cast<int>(std::floor((obj->getPosition().Y + obj->getInteractionRadius()) * inv_size + 0.5));
        const int zs = static_cast<int>(std::floor((obj->getPosition().Z - obj->getInteractionRadius()) * inv_size - 0.5));
        const int ze = static_cast<int>(std::floor((obj->getPosition().Z + obj->getInteractionRadius()) * inv_size + 0.5));
        
        logger(VERBOSE, "Level = % grid cells [%,%] x [%,%] x [%,%]", level, xs, xe, ys, ye, zs, ze);
        for (int x = xs; x <= xe; ++x)
        {
            for (int y = ys; y <= ye; ++y)
            {
                for (int z = zs; z <= ze; ++z)
                {
                    if (hGridHasContactsInTargetCell(x, y, z, level, obj))
                    {
                        return true;
                    }
                }
            }
        }
    } //end for level
    
    return false;
}

#ifdef CONTACT_LIST_HGRID

void Mercury3D::InsertCell(int x, int y, int z, unsigned int l, BaseParticle *obj)
{   
    // Loop through all objects in the bucket to find nearby objects
    unsigned int bucket = getHGrid()->computeHashBucketIndex(x,y,z,l);
    BaseParticle *p = getHGrid()->getFirstBaseParticleInBucket(bucket);

    while (p!=nullptr)
    {   
        if ((p->getHGridX() == x) && (p->getHGridY() == y) && (p->getHGridZ() == z) && (p->getHGridLevel() == l) && (obj!=p))
        {   
            getPossibleContactList().add_PossibleContact(obj,p);
        }
        p = p->getHGridNextObject();
    }
}

void Mercury3D::InsertObjAgainstGrid(BaseParticle *obj)
{   
    Mdouble inv_size;
    int occupiedLevelsMask_ = getHGrid()->getOccupiedLevelsMask();

    inv_size=getHGrid()->getInvCellSize(obj->getHGridLevel());

    double ownXMin = (obj->getHGridX() - 0.5) * getHGrid()->getCellSize(obj->getHGridLevel());
    double ownXMax = (obj->getHGridX() + 1.5) * getHGrid()->getCellSize(obj->getHGridLevel());
    double ownYMin = (obj->getHGridY() - 0.5) * getHGrid()->getCellSize(obj->getHGridLevel());
    double ownYMax = (obj->getHGridY() + 1.5) * getHGrid()->getCellSize(obj->getHGridLevel());
    double ownZMin = (obj->getHGridZ() - 0.5) * getHGrid()->getCellSize(obj->getHGridLevel());
    double ownZMax = (obj->getHGridZ() + 1.5) * getHGrid()->getCellSize(obj->getHGridLevel());

    for (int level = 0; level < getHGrid()->getNumberOfLevels(); occupiedLevelsMask_ >>= 1, level++)
    {   
        // If no objects in rest of grid, stop now
        if (occupiedLevelsMask_ == 0)
        {   
            break;
        }

        // If no objects at this level, go on to the next level
        if ((occupiedLevelsMask_ & 1) == 0)
        {   
            continue;
        }

        // Treat level as a third dimension coordinate
        inv_size = getHGrid()->getInvCellSize(level);

        int xs, xe, ys, ye, zs, ze;
        xs=static_cast<int>(std::floor(ownXMin * inv_size - 0.5));
        xe=static_cast<int>(std::floor(ownXMax * inv_size + 0.5));
        ys=static_cast<int>(std::floor(ownYMin * inv_size - 0.5));
        ye=static_cast<int>(std::floor(ownYMax * inv_size + 0.5));
        zs=static_cast<int>(std::floor(ownZMin * inv_size - 0.5));
        ze=static_cast<int>(std::floor(ownZMax * inv_size + 0.5));

        for(int x=xs; x<=xe; ++x)
        {   
            for(int y=ys; y<=ye; ++y)
            {   
                for(int z=zs; z<=ze; ++z)
                {   
                    InsertCell(x, y, z, level, obj);
                }
            }
        }
    } //end for level
}
#endif