Mercury3DClump.h

Go to the documentation of this file.

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#include "Mercury3D.h"
#include "Walls/InfiniteWall.h"
#include "Species/LinearViscoelasticFrictionSpecies.h"
#include "Particles/ClumpParticle.h"
#include "Boundaries/PeriodicBoundary.h"
#include <stdlib.h>
 
class Mercury3Dclump : public Mercury3D
{
public:
    explicit  Mercury3Dclump()
    {
 
    }
 
    void computeInternalForce(BaseParticle* P1, BaseParticle* P2) override
    {
        // Quit if:
        // 1) at least one particle is master
        // 2) pebbles of the same clump
        if (P1->isClump() || P2->isClump())
            return;
        if (P1->isPebble() && P2->isPebble() && P1->getClump() == P2->getClump())
            return;
 
        Mercury3D::computeInternalForce(P1, P2);
    }
 
    void computeForcesDueToWalls(BaseParticle* pI, BaseWall* w) override
    {
        // No interaction between walls and master particles
        if (pI->isClump())
            return;
 
        Mercury3D::computeForcesDueToWalls(pI, w);
    }
 
    void computeAllForces() override
    {
        Mercury3D::computeAllForces();
 
        #pragma omp for schedule(dynamic)
        for (BaseParticle* p : particleHandler)
        {
            if ((p->isPebble()) && (p->getPeriodicFromParticle() == nullptr))
            {
                p->getClump()->addForce(p->getForce());
                Vec3D lever = p->getPosition()- p->getClump()->getPosition();
 
                // Patch to fix lever - we do not create ghost masters, that's why this workaround is needed
                // Should affect all boundaries derived from PeriodicBoundary and do not affect other types
                for (auto b : boundaryHandler)
                {
                    auto pb = dynamic_cast<PeriodicBoundary*>(b);
                    if (pb){
                        Vec3D shift = pb->getShift();
                        if (lever.getLength() > (lever - shift).getLength()) lever -= shift;
                        if (lever.getLength() > (lever + shift).getLength()) lever += shift;
                    }
                }
                // End patch
 
                Vec3D torque = p->getTorque();
                torque += Vec3D::cross(lever, p->getForce());
                p->getClump()->addTorque(torque);
            }
        }
    }
 
    bool checkClumpForInteraction(BaseParticle& particle)
    {
        ClumpParticle* mp = dynamic_cast<ClumpParticle*>(&particle);
        if (mp->isClump()){
            for (int i = 0; i< mp->NPebble(); i++){
                for (BaseParticle *q: particleHandler) {
                    if (!q->isClump()) { // check every slave vs every non-master intersection
                        if (Vec3D::getDistanceSquared(q->getPosition(), mp->getPebblePositions()[i])
                            < q->getRadius() + mp->getPebbleRadii()[i]) {
                            return false;
                        }
                    }
                }
            }
        }
        return true;
    }
 
    bool checkClumpForInteractionPeriodic(BaseParticle& particle)
    {
        // Note that this implementation only check for interaction with particles
        bool NP = true;
        // Periodic case
        for (BaseBoundary* b : boundaryHandler)
        {
            PeriodicBoundary* pb = dynamic_cast<PeriodicBoundary*>(b);
            if (pb) NP = false;
            if (pb && particleHandler.getNumberOfObjects() > 0 )
            {
                ClumpParticle* mp = dynamic_cast<ClumpParticle*>(&particle);
                if (mp->isClump()){
                    for (int i = 0; i< mp->NPebble(); i++){
                        for (BaseParticle *q: particleHandler) {
                            if (!q->isClump()) { // check every slave vs every non-master intersection
                                if (Vec3D::getDistanceSquared(q->getPosition(), mp->getPebblePositions()[i])
                                    < (q->getMaxInteractionRadius() + mp->getPebbleRadii()[i])
                                    * (q->getMaxInteractionRadius() + mp->getPebbleRadii()[i]) ) {
                                    return false;
                                }
                                BaseParticle *q1 = q->copy(); // check every slave vs non-master ghost intersection
                                pb->shiftPosition(q1);
                                Mdouble dist2 = Vec3D::getDistanceSquared(q1->getPosition(), mp->getPebblePositions()[i]);
                                Mdouble dist02 = (q1->getMaxInteractionRadius() + mp->getPebbleRadii()[i])
                                        * (q1->getMaxInteractionRadius() + mp->getPebbleRadii()[i]);
                                delete q1;
                                if (dist2 < dist02) {
                                    return false;
                                }
                            }
                        }
                    }
                }
            }
        }
        // Non-periodic case
        if ((NP)&&(particleHandler.getNumberOfObjects() > 0 ))
        {
            ClumpParticle* mp = dynamic_cast<ClumpParticle*>(&particle);
            if (mp->isClump()){
                for (int i = 0; i< mp->NPebble(); i++){
                    for (BaseParticle *q: particleHandler) {
                        if (!q->isClump()) { // check every slave vs every non-master intersection
                            if (Vec3D::getDistanceSquared(q->getPosition(), mp->getPebblePositions()[i])
                                < (q->getRadius() + mp->getPebbleRadii()[i]) * (q->getRadius() + mp->getPebbleRadii()[i])) {
                                return false;
                            }
                        }
                    }
                }
            }
        }
        return true;
    }
};