BaseInteraction.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 "BaseInteraction.h"
#include "InteractionHandler.h"
#include "Particles/BaseParticle.h"
#include "Species/BaseSpecies.h"
#include "DPMBase.h"
#include<iomanip>
#include<fstream>

BaseInteraction::BaseInteraction(BaseInteractable* P, BaseInteractable* I, Mdouble timeStamp)
    : BaseObject()
{
    P_ = P;
    P->addInteraction(this);
    I_ = I;
    I->addInteraction(this);
    normal_.setZero();
    overlap_ = 0;
    timeStamp_ = timeStamp;
    species_ = 0;
    force_.setZero();
    torque_.setZero();
#ifdef DEBUG_CONSTRUCTOR
    std::cout<<"BaseInteraction::BaseInteraction() finished"<<std::endl;
#endif
}

BaseInteraction::BaseInteraction(const BaseInteraction &p)
    : BaseObject(p)
{
    P_ = p.P_;
    I_ = p.I_;
    normal_ = p.normal_;
    overlap_ = p.overlap_;
    force_ = p.force_;
    torque_ = p.torque_;
    species_ = p.species_;
    timeStamp_ = p.timeStamp_;
}

BaseInteraction::BaseInteraction()
        : BaseObject()
{
    P_ = nullptr;
    I_ = nullptr;
    normal_.setZero();
    overlap_ = 0;
    timeStamp_ = 0;
    species_ = 0;
    force_.setZero();
    torque_.setZero();
}

BaseInteraction::~BaseInteraction()
{
    P_->removeInteraction(this);
    I_->removeInteraction(this);
}

void BaseInteraction::write(std::ostream& os) const
{
    os << getName();
    if (dynamic_cast<BaseParticle*>(I_) != nullptr)
        os << " particleIds " << P_->getId() << " " << I_->getId();
    else
        os << " particleWallIds " << P_->getId() << " " << I_->getId();
    os <<" timeStamp "<<timeStamp_<< " force " << force_ << " torque " << torque_;
}

void BaseInteraction::read(std::istream& is)
{
    //the rest gets read by the interaction handler
    std::string dummy;
    is >> dummy >> force_ >> dummy >> torque_;
}

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

void BaseInteraction::setNormal(Vec3D normal)
{
    normal_ = normal;
}

void BaseInteraction::setDistance(Mdouble distance)
{
    distance_ = distance;
}

void BaseInteraction::setOverlap(Mdouble overlap)
{
    overlap_ = overlap;
}

/*
 * \details set the contact point between the two interactable objects involved
 * \param[in] contactPoint  Vec3D vector which will become the new contact point.
 */
void BaseInteraction::setContactPoint(Vec3D contactPoint)
{
    contactPoint_ = contactPoint;
}

void BaseInteraction::setTimeStamp(Mdouble timeStamp)
{
    timeStamp_ = timeStamp;
}

void BaseInteraction::setHandler(InteractionHandler* handler)
{
    handler_ = handler;
}

InteractionHandler* BaseInteraction::getHandler() const
{
    return handler_;
}

void BaseInteraction::removeFromHandler()
{
    handler_->removeObject(getIndex());
}

void BaseInteraction::copySwitchPointer(const BaseInteractable* original, BaseInteractable* ghost) const
{
    //Copy the interaction of ghost
    BaseInteraction* C = copy();
    //Finds out which of P and I is that the particle from whom the ghost is begin created.
    //The object being interacted with is set to P
    if (C->getP() == original)
    {
        //Reverse some force history
        C->reverseHistory();
        //Set the P to the original particle 
        C->P_ =C->getI();
    }
    //The new ghost particle is set to I in the interaction. 
    C->I_ = ghost;

    //Add the the interaction to both original and the ghost
    C->getP()->addInteraction(C);
    C->getI()->addInteraction(C);
    handler_->addObject(C);
}

const Vec3D& BaseInteraction::getForce() const
{
    return force_;
}

const Vec3D& BaseInteraction::getTorque() const
{
    return torque_;
}

const Vec3D& BaseInteraction::getNormal() const
{
    return normal_;
}

const Vec3D& BaseInteraction::getContactPoint() const
{
    return contactPoint_;
}

Mdouble BaseInteraction::getOverlap() const
{
    return overlap_;
}

Mdouble BaseInteraction::getContactRadius() const
{
    return sqrt(2.0*getEffectiveRadius()*getOverlap());
}

BaseInteractable* BaseInteraction::getP()
{
    return P_;
}

BaseInteractable* BaseInteraction::getI()
{
    return I_;
}

const BaseInteractable* BaseInteraction::getP() const
{
    return P_;
}

const BaseInteractable* BaseInteraction::getI() const
{
    return I_;
}

Mdouble BaseInteraction::getTimeStamp() const
{
    return timeStamp_;
}

void BaseInteraction::integrate(Mdouble timeStep UNUSED)
{

}

void BaseInteraction::setSpecies(BaseSpecies *species)
{
    species_ = species;
}

void BaseInteraction::setP(BaseInteractable* P)
{
    P_->removeInteraction(this);
    P_=P;
    P_->addInteraction(this);
}

void BaseInteraction::setI(BaseInteractable* I)
{
    I_->removeInteraction(this);
    I_=I;
    I_->addInteraction(this);
}

void BaseInteraction::writeToFStat(std::ostream& os) const
{
    BaseParticle* IParticle = dynamic_cast<BaseParticle*>(I_);
    BaseParticle* PParticle = dynamic_cast<BaseParticle*>(P_);

    Vec3D tangentialForce = getTangentialForce();
    Mdouble tangentialOverlap = getTangentialOverlap();

    Mdouble scalarNormalForce = Vec3D::dot(force_, getNormal());
    Mdouble scalarTangentialForce = tangentialForce.getLength();
    Vec3D tangential;
    if (scalarTangentialForce!=0.0)
        tangential = tangentialForce/scalarTangentialForce;
    else
        tangential = Vec3D(0.0,0.0,0.0);

    Vec3D centre;
    if (IParticle!=0)
        centre = getP()->getPosition() - normal_ * (PParticle->getRadius() + IParticle->getRadius() - overlap_)/2.0;
        //centre = 0.5 * (getP()->getPosition() + getI()->getPosition());
    else
        centre = getP()->getPosition() - normal_ * (PParticle->getRadius() - overlap_);

    if (PParticle!=0 && !PParticle->isFixed())
    {
        os << timeStamp_
            << " " << P_->getIndex()
            << " " << static_cast<int>((IParticle==0?(-I_->getIndex()-1):I_->getIndex()))
            << " " << centre
            << " " << getOverlap()
            << " " << tangentialOverlap
            << " " << scalarNormalForce
            << " " << scalarTangentialForce
            << " " << (IParticle==0?-normal_:normal_)
            << " " << (IParticle==0?-tangential:tangential) << std::endl;
    }
    if (IParticle!=0 && !IParticle->isFixed() && IParticle->getPeriodicFromParticle()==0)
    {
        os << timeStamp_
            << " " << I_->getIndex()
            << " " << P_->getIndex()
            << " " << centre
            << " " << getOverlap()
            << " " << tangentialOverlap
            << " " << scalarNormalForce
            << " " << scalarTangentialForce
            << " " << -normal_
            << " " << -tangential << std::endl;
    }

}

/*
 * \details Returns the distance between the two interactable objects involved
 *          in the interaction.
 * \return  Mdouble which is the distance between the two interacting objects.
 */
Mdouble BaseInteraction::getDistance() const
{
    return distance_;
}

Mdouble BaseInteraction::getTangentialOverlap() const
{
    return 0;
}

const Vec3D BaseInteraction::getTangentialForce() const
{
    return Vec3D(0.0,0.0,0.0);
}

const Vec3D& BaseInteraction::getRelativeVelocity() const
{
    return relativeVelocity_;
}

Mdouble BaseInteraction::getNormalRelativeVelocity() const
{
    return normalRelativeVelocity_;
}

Mdouble BaseInteraction::getAbsoluteNormalForce() const
{
    return absoluteNormalForce_;
}

void BaseInteraction::addForce(Vec3D force)
{
    force_+=force;
}

void BaseInteraction::addTorque(Vec3D torque)
{
    torque_+=torque;
}

void BaseInteraction::setForce(Vec3D force)
{
    force_=force;
}

void BaseInteraction::setTorque(Vec3D torque)
{
    torque_=torque;
}

void BaseInteraction::setRelativeVelocity(Vec3D relativeVelocity)
{
    relativeVelocity_=relativeVelocity;
}

void BaseInteraction::setNormalRelativeVelocity(Mdouble normalRelativeVelocity)
{
    normalRelativeVelocity_=normalRelativeVelocity;
}

void BaseInteraction::setAbsoluteNormalForce(Mdouble absoluteNormalForce)
{
    absoluteNormalForce_ = absoluteNormalForce;
}

const BaseSpecies* BaseInteraction::getBaseSpecies() const
{
    return species_;
}

void BaseInteraction::computeForce()
{}

Mdouble BaseInteraction::getElasticEnergy() const
{
    return 0.0;
}

void BaseInteraction::reverseHistory()
{
}

unsigned int BaseInteraction::getMultiContactIdentifier() const
{
    return multiContactIdentifier_;
}

void BaseInteraction::setMultiContactIdentifier(unsigned int multiContactIdentifier)
{
    multiContactIdentifier_ = multiContactIdentifier;
}


void BaseInteraction::rotateHistory(Matrix3D& rotationMatrix)
{
    contactPoint_=rotationMatrix*contactPoint_;
    relativeVelocity_=rotationMatrix*relativeVelocity_;
    force_=rotationMatrix*force_;
    torque_=rotationMatrix*torque_;
    normal_=rotationMatrix*normal_;
}

Mdouble BaseInteraction::getEffectiveRadius() const
{
    const BaseParticle* PParticle = dynamic_cast<const BaseParticle*>(getP());
    const BaseParticle* IParticle = dynamic_cast<const BaseParticle*>(getI());
    if (PParticle==nullptr)
    {
        std::cerr << "BaseInteraction::getEffectiveCorrectedRadius(): first interactable P is not a particle" << std::endl;
        exit(-1);
    }
    //Compute the reduced diameter
    if (IParticle==nullptr) //if particle-wall
    {
        return PParticle->getRadius();
    }
    else
    {
        Mdouble radiusP = PParticle->getRadius();
        Mdouble radiusI = IParticle->getRadius();
        return radiusP * radiusI / (radiusP + radiusI);
    }
}

Mdouble BaseInteraction::getEffectiveMass() const
{
    const BaseParticle* PParticle = dynamic_cast<const BaseParticle*>(getP());
    const BaseParticle* IParticle = dynamic_cast<const BaseParticle*>(getI());
    if (PParticle==nullptr)
    {
        std::cerr << "BaseInteraction::getEffectiveCorrectedRadius(): first interactable P is not a particle" << std::endl;
        exit(-1);
    }
    //Compute the reduced diameter
    if (IParticle==nullptr) //if particle-wall
    {
        return PParticle->getMass();
    }
    else
    {
        Mdouble massP = PParticle->getMass();
        Mdouble massI = IParticle->getMass();
        return massP * massI / (massP + massI);
    }
}

Mdouble BaseInteraction::getEffectiveCorrectedRadius()
{
    BaseParticle* PParticle = dynamic_cast<BaseParticle*>(getP());
    BaseParticle* IParticle = dynamic_cast<BaseParticle*>(getI());
    if (PParticle==nullptr)
    {
        std::cerr << "BaseInteraction::getEffectiveCorrectedRadius(): first interactable P is not a particle" << std::endl;
        exit(-1);
    }
    //Compute the reduced diameter
    if (IParticle==nullptr) //if particle-wall
    {
        return PParticle->getRadius() - 0.5*getOverlap();
    }
    else
    {
        Mdouble radiusP = PParticle->getRadius() - 0.5*getOverlap();
        Mdouble radiusI = IParticle->getRadius() - 0.5*getOverlap();
        return radiusP * radiusI / (radiusP + radiusI);
    }
}

void BaseInteraction::actionsAfterTimeStep()
{
}

void BaseInteraction::gatherContactStatistics()
{
    BaseParticle* IParticle = dynamic_cast<BaseParticle*>(I_);
    BaseParticle* PParticle = dynamic_cast<BaseParticle*>(P_);

    Vec3D tangentialForce = getTangentialForce();
    Mdouble tangentialOverlap = getTangentialOverlap();

    Mdouble scalarNormalForce = Vec3D::dot(force_, getNormal());
    Mdouble scalarTangentialForce = tangentialForce.getLength();
    Vec3D tangential;
    if (scalarTangentialForce!=0.0)
        tangential = tangentialForce/scalarTangentialForce;
    else
        tangential = Vec3D(0.0,0.0,0.0);

    Vec3D centre;
    if (IParticle!=nullptr)
        centre = getP()->getPosition() - normal_ * (PParticle->getRadius() + IParticle->getRadius() - overlap_)/2.0;
    else
        centre = getP()->getPosition() - normal_ * (PParticle->getRadius() - overlap_);

    if (PParticle!=nullptr && !PParticle->isFixed())
    {
        getHandler()->getDPMBase()->gatherContactStatistics(
            P_->getIndex(),
            static_cast<int>((IParticle==0?(-I_->getIndex()-1):I_->getIndex())),
            centre,
            getOverlap(),
            tangentialOverlap,
            scalarNormalForce,
            scalarTangentialForce,
            (IParticle==0?-normal_:normal_),
            (IParticle==0?-tangential:tangential));
    }
    if (IParticle!= nullptr && !IParticle->isFixed() && IParticle->getPeriodicFromParticle()==0)
    {
        getHandler()->getDPMBase()->gatherContactStatistics(
            I_->getIndex(),
            static_cast<int>(P_->getIndex()),
            centre,
            getOverlap(),
            tangentialOverlap,
            scalarNormalForce,
            scalarTangentialForce,
            -normal_,
            -tangential);

    }
}

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

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

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

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