BaseParticle.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 <Species/LinearViscoelasticSpecies.h>
#include "InteractionHandler.h"
#include "Particles/BaseParticle.h"
#include "Interactions/BaseInteraction.h"
#include "Species/ParticleSpecies.h"
#include "ParticleHandler.h"
#include "DPMBase.h"

BaseParticle::BaseParticle()
{
    handler_ = nullptr;
    displacement_.setZero();
    radius_ = 1.0;
    mass_ = invMass_ = 1.0;
    inertia_ = invInertia_ = 1.0;
    HGridNextObject_ = nullptr;

    periodicFromParticle_ = nullptr;
#ifdef CONTACT_LIST_HGRID
    firstPossibleContact = nullptr;
#endif
    HGridNextObject_ = nullptr;
    HGridPrevObject_ = nullptr;
    HGridLevel_ = 99999;
    HGridX_ = 99999;
    HGridY_ = 99999;
    HGridZ_ = 99999;

    logger(DEBUG, "BaseParticle::BaseParticle() finished");
}

BaseParticle::BaseParticle(const BaseParticle &p)
: BaseInteractable(p)
{
    handler_ = nullptr;
    displacement_ = p.displacement_;
    radius_ = p.radius_;
    mass_ = p.getMass();
    invMass_ = p.getInvMass();
    inertia_ = p.getInertia();
    invInertia_ = p.getInvInertia();

    HGridNextObject_ = nullptr;
    HGridPrevObject_ = nullptr;
    HGridX_ = 99999;
    HGridY_ = 99999;
    HGridZ_ = 99999;
    HGridLevel_ = p.HGridLevel_;

    periodicFromParticle_ = p.periodicFromParticle_;
#ifdef CONTACT_LIST_HGRID
    firstPossibleContact = nullptr;
#endif
    logger(DEBUG, "BaseParticle::BaseParticle(BaseParticle &p) finished");
}

BaseParticle::~BaseParticle()
{
    if (getHandler() != nullptr)
    {
        getHandler()->checkExtremaOnDelete(this);
    }
    logger(DEBUG, "BaseParticle::~BaseParticle() of particle % finished.", getId());
}

BaseParticle* BaseParticle::copy() const
{
    return new BaseParticle(*this);
}

Mdouble BaseParticle::getVolume() const
{
    if (handler_ == nullptr)
    {
        logger(ERROR, "[BaseParticle::getVolume] no particle handler specified");
        return 0;
    }
    switch (getParticleDimensions())
    {
    case 3:
        return (4.0 / 3.0 * constants::pi * radius_ * radius_ * radius_);
    case 2:
        return (constants::pi * radius_ * radius_);
    case 1:
        return (2.0 * radius_);
    default:
        logger(ERROR, "[BaseParticle::getVolume] dimension of the particle is not set");
        return 0;
    }
}

void BaseParticle::fixParticle()
{
    mass_ = 1e20;
    invMass_ = 0.0;
    inertia_ = 1e20;
    invInertia_ = 0.0;
    setVelocity(Vec3D(0.0, 0.0, 0.0));
    setAngularVelocity(Vec3D(0.0, 0.0, 0.0));
}

bool BaseParticle::isFixed() const
{
    return (invMass_ == 0.0);
}

void BaseParticle::unfix()
{
    invMass_ = 1.0;
    getSpecies()->computeMass(this);
}

void BaseParticle::write(std::ostream& os) const
{
    BaseInteractable::write(os);
    os << " radius " << radius_
        << " invMass " << invMass_
        << " invInertia " << invInertia_;
}

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

void BaseParticle::read(std::istream& is)
{
    BaseInteractable::read(is);
    std::string dummy;
    is >> dummy >> radius_ >> dummy >> invMass_ >> dummy >> invInertia_;
    if (invMass_ != 0.0)
        mass_ = 1.0 / invMass_;
    else
        mass_ = 1e20;
    if (invInertia_ != 0.0)
        inertia_ = 1.0 / invInertia_;
    else
        inertia_ = 1e20;
}

void BaseParticle::oldRead(std::istream& is)
{

    unsigned int indSpecies;
    Vec3D orientation;
    Vec3D position;
    is >> invMass_ >> invInertia_ >> indSpecies;
    setPosition(position);
    setOrientation(orientation);
    setIndSpecies(indSpecies);
    if (invMass_ != 0.0)
        mass_ = 1.0 / invMass_;
    else
        mass_ = 1e20;
    if (invInertia_ != 0.0)
        inertia_ = 1.0 / invInertia_;
    else
        inertia_ = 1e20;
}

void BaseParticle::printHGrid(std::ostream& os) const
{
    os << "Particle( HGRID_Level:" << HGridLevel_
        << ", HGRID_x:" << HGridX_
        << ", HGRID_y:" << HGridY_
        << ", HGRID_z:" << HGridZ_
        << ")";
}

unsigned int BaseParticle::getHGridLevel() const
{
    return HGridLevel_;
}

BaseParticle* BaseParticle::getHGridNextObject() const
{
    return HGridNextObject_;
}

BaseParticle* BaseParticle::getHGridPrevObject() const
{
    return HGridPrevObject_;
}

#ifdef CONTACT_LIST_HGRID

PossibleContact* BaseParticle::getFirstPossibleContact() const
{
    return firstPossibleContact;
}
#endif

int BaseParticle::getHGridX() const
{
    return HGridX_;
}

int BaseParticle::getHGridY() const
{
    return HGridY_;
}

int BaseParticle::getHGridZ() const
{
    return HGridZ_;
}

Mdouble BaseParticle::getInertia() const
{
    return inertia_;
}

Mdouble BaseParticle::getInvInertia() const
{
    return invInertia_;
}

Mdouble BaseParticle::getInvMass() const
{
    return invMass_;
}

Mdouble BaseParticle::getKineticEnergy() const
{
    return 0.5 * mass_ * getVelocity().getLengthSquared();
}

Mdouble BaseParticle::getMass() const
{
    return mass_;
}

BaseParticle* BaseParticle::getPeriodicFromParticle() const
{
    return periodicFromParticle_;
}

Mdouble BaseParticle::getRadius() const
{
    return radius_;
}

Mdouble BaseParticle::getInteractionRadius() const
{
    return radius_ + getSpecies()->getInteractionDistance() * 0.5;
}

Mdouble BaseParticle::getWallInteractionRadius() const
{
    return radius_ + getSpecies()->getInteractionDistance();
}

const Vec3D& BaseParticle::getDisplacement() const
{
    return displacement_;
}

const Vec3D& BaseParticle::getPreviousPosition() const
{
    return previousPosition_;
}

const Vec3D BaseParticle::getDisplacement2(Mdouble xmin, Mdouble xmax, Mdouble ymin, Mdouble ymax, Mdouble zmin, Mdouble zmax, Mdouble t) const
{
    Vec3D disp = getPosition() - getPreviousPosition();
    if (xmax > xmin && fabs(disp.X) > .5 * (xmax - xmin))
    {
        if (disp.X > 0)
            disp.X -= xmax - xmin;
        else
            disp.X += xmax - xmin;
    }
    if (ymax > ymin && fabs(disp.Y) > .5 * (ymax - ymin))
    {
        if (disp.Y > 0)
            disp.Y -= ymax - ymin;
        else
            disp.Y += ymax - ymin;
    }
    if (zmax > zmin && fabs(disp.Z) > .5 * (zmax - zmin))
    {
        if (disp.Z > 0)
            disp.Z -= zmax - zmin;
        else
            disp.Z += zmax - zmin;
    }
    disp /= t;
    return disp;
}

void BaseParticle::setInertia(const Mdouble newInertia)
{
    if (newInertia >= 0)
    {
        inertia_ = newInertia;
        invInertia_ = 1.0 / newInertia;
    }
    else
    {
        logger(ERROR, "Error in set_inertia (%)", newInertia);
    }
}

void BaseParticle::setInfiniteInertia()
{
    inertia_ = 1e20;
    invInertia_ = 0;
} //> i.e. no rotations

void BaseParticle::setPeriodicFromParticle(BaseParticle* p)
{
    periodicFromParticle_ = p;
}

void BaseParticle::setHGridX(const int x)
{
    HGridX_ = x;
}

void BaseParticle::setHGridY(const int y)
{
    HGridY_ = y;
}

void BaseParticle::setHGridZ(const int z)
{
    HGridZ_ = z;
}

void BaseParticle::setHGridLevel(const unsigned int level)
{
    HGridLevel_ = level;
}

void BaseParticle::setHGridNextObject(BaseParticle* p)
{
    HGridNextObject_ = p;
}

void BaseParticle::setHGridPrevObject(BaseParticle* p)
{
    HGridPrevObject_ = p;
}

#ifdef CONTACT_LIST_HGRID

void BaseParticle::setFirstPossibleContact(PossibleContact* PC)
{
    firstPossibleContact = PC;
}
#endif

void BaseParticle::setRadius(const Mdouble radius)
{
    radius_ = radius;
    if (getHandler())
    {
        getSpecies()->computeMass(this);
        getHandler()->checkExtrema(this);
    }
}

void BaseParticle::setMass(const Mdouble mass)
{
    logger(WARN, "WARNING: Do not use particle->setMass, instead use "
           "particleSpecies->computeMass, since this function can cause "
           "inconsistencies between the mass, density and radius of this particle!");
    if (mass >= 0.0)
    {
        if (invMass_ != 0.0) //InvMass=0 is a flag for a fixed particle
        {
            mass_ = mass;
            invMass_ = 1.0 / mass;
        }
    }
    else
    {
        logger(ERROR, "Error in BaseParticle::setMass, the given mass to be set must be positive.");
    }
}

void BaseParticle::setMassForP3Statistics(const Mdouble mass)
{
    if(mass >= 0.0)
    {
        if(invMass_ != 0.0) //InvMass=0 is a flag for a fixed particle
        {
            mass_ = mass;
            invMass_ = 1.0 / mass;
        }
    }
    else
    {
        logger(ERROR, "Error in BaseParticle::setMass, the given mass to be set must be positive.");
    }
} 

void BaseParticle::setDisplacement(const Vec3D& disp)
{
    displacement_ = disp;
}

void BaseParticle::setPreviousPosition(const Vec3D& pos)
{
    previousPosition_ = pos;
}

void BaseParticle::movePrevious(const Vec3D& posMove)
{
    previousPosition_ += posMove;
}

void BaseParticle::accelerate(const Vec3D& vel)
{
    addVelocity(vel);
}

void BaseParticle::angularAccelerate(const Vec3D& angVel)
{
    addAngularVelocity(angVel);
}

void BaseParticle::addDisplacement(const Vec3D& addDisp)
{
    displacement_ += addDisp;
}

void BaseParticle::setHandler(ParticleHandler* handler)
{
    handler_ = handler;
    setSpecies(getHandler()->getDPMBase()->speciesHandler.getObject(getIndSpecies()));
}

ParticleHandler* BaseParticle::getHandler() const
{
    return handler_;
}

std::vector<BaseInteraction*> BaseParticle::getInteractionWith(BaseParticle *P, Mdouble timeStamp, InteractionHandler* interactionHandler)
{
    //get the normal (from P away from the contact)
    Vec3D branchVector = P->getPosition() - getPosition();
    //Get the square of the distance between particle i and particle j
    Mdouble distanceSquared = Vec3D::getLengthSquared(branchVector);
    Mdouble sumOfInteractionRadii = P->getInteractionRadius() + getInteractionRadius();
    std::vector<BaseInteraction*> interactions;
    if (distanceSquared < (sumOfInteractionRadii * sumOfInteractionRadii))
    {
        BaseInteraction* C = interactionHandler->getInteraction(P, this, timeStamp);
        Mdouble distance = std::sqrt(distanceSquared);
        C->setNormal(branchVector / distance);
        C->setOverlap(P->getRadius() + getRadius() - distance);
        C->setDistance(distance);
        C->setContactPoint(P->getPosition() - (P->getRadius() - 0.5 * C->getOverlap()) * C->getNormal());
        //Mdouble ratio=P->getRadius()/(getRadius()+P->getRadius());
        //C->setContactPoint(P->getPosition() - (P->getRadius() - ratio * C->getOverlap()) * C->getNormal());
        interactions.push_back(C);
    }
    return interactions;
}




void BaseParticle::integrateBeforeForceComputation(double time, double timeStep)
{
    if (getInvMass() == 0.0)
    {
        BaseInteractable::integrateBeforeForceComputation(time, timeStep);
    }
    else
    {
        accelerate(getForce() * getInvMass() * 0.5 * timeStep);
        setDisplacement(getVelocity() * timeStep);
        move(getDisplacement());
        getHandler()->getDPMBase()->hGridUpdateMove(this, getDisplacement().getLength());
        if (getHandler()->getDPMBase()->getRotation())
        {
            angularAccelerate(getTorque() * getInvInertia() * 0.5 * timeStep);
            rotate(getAngularVelocity() * timeStep);
        }
    }
}

void BaseParticle::integrateAfterForceComputation(double time, double timeStep)
{
    if (getInvMass() == 0.0)
    {
        BaseInteractable::integrateAfterForceComputation(time, timeStep);
    }
    else
    {
        accelerate(getForce() * getInvMass() * 0.5 * timeStep);
        if (getHandler()->getDPMBase()->getRotation())
        {
            angularAccelerate(getTorque() * getInvInertia() * 0.5 * timeStep);
        }
    }
}

unsigned int BaseParticle::getParticleDimensions() const
{
    return getHandler()->getDPMBase()->getParticleDimensions();
}

void BaseParticle::setIndSpecies(unsigned int indSpecies)
{
    if (handler_ != nullptr)
    {
        setSpecies(handler_->getDPMBase()->speciesHandler.getObject(indSpecies));
    }
    else
    {
        BaseInteractable::setIndSpecies(indSpecies);
        logger(ERROR, "setIndSpecies called on a particle with no particle handler.\n"
               "Therefore I can't request the given species from the species handler.\n"
               " PartID = %", getId());
    }
}

void BaseParticle::setSpecies(const ParticleSpecies* species)
{
    BaseInteractable::setSpecies(species);
    //set pointer to the ParticleHandler handler_, which is needed to retrieve 
    //species information
    //\todo maybe these if statements should throw warnings
    if (handler_ == nullptr)
    {
        SpeciesHandler* sH = species->getHandler();
        DPMBase* dB = sH->getDPMBase();
        if (dB != nullptr)
        {
            setHandler(&dB->particleHandler);
        }
    }
}

unsigned BaseParticle::getNumberOfFieldsVTK() const {
    return 0;
}

std::string BaseParticle::getTypeVTK(unsigned i) const {
    return "";
}

std::string BaseParticle::getNameVTK(unsigned i) const {
    return "";
}

std::vector<Mdouble> BaseParticle::getFieldVTK(unsigned i) const {
    return std::vector<Mdouble>();
}