BaseInteraction.cc

Go to the documentation of this file.

//Copyright (c) 2013-2020, 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 "DPMBase.h"

BaseInteraction::BaseInteraction(BaseInteractable* P, BaseInteractable* I, unsigned timeStamp)
        : BaseObject()
{
    P_ = P;
    I_ = I;
    normal_.setZero();
    contactPoint_.setNaN();
    overlap_ = 0;
    timeStamp_ = timeStamp;
    species_ = nullptr;
    force_.setZero();
    torque_.setZero();
    if (P->getSpecies()->getHandler()->getDPMBase()->getInteractionFile().getFileType() == FileType::ONE_FILE)
    {
        writeInteraction(P->getSpecies()->getHandler()->getDPMBase()->getInteractionFile().getFstream(), true);
    }
    lagrangeMultiplier_ = 0;

#ifdef DEBUG_CONSTRUCTOR
    std::cout<<"BaseInteraction::BaseInteraction() finished"<<std::endl;
#endif
}


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

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_;
    lagrangeMultiplier_ = p.lagrangeMultiplier_;
    contactPoint_ = p.contactPoint_;
    
//    InteractionHandler* handler_;
//    BaseInteractable* P_;
//    BaseInteractable* I_;
//    unsigned int identificationP_;
//    int identificationI_;
//    bool isWallInteraction_;
//    Vec3D contactPoint_;
//    Vec3D relativeVelocity_;
//    Mdouble normalRelativeVelocity_;
//    Mdouble absoluteNormalForce_;
//    Mdouble distance_;
//    Vec3D force_;
//    Vec3D torque_;
//    unsigned timeStamp_;
//    Vec3D normal_;
//    Mdouble overlap_;
//    const BaseSpecies* species_;
//    Mdouble lagrangeMultiplier_;
//    unsigned multiContactIdentifier_;

}

BaseInteraction::~BaseInteraction()
{
#if MERCURY_USE_MPI
    if (P_ == nullptr)
    {
        logger(DEBUG,"Destroying a fictitious interaction used in MPI transmissions");
    }
    else
    {
#endif
    logger.assert(P_ != nullptr, "Trying to destroy an interaction with P_ = nullptr");
    logger.assert(I_ != nullptr, "Trying to destroy an interaction with I_ = nullptr");
    File& interactionFile = getHandler()->getDPMBase()->getInteractionFile();
    if (interactionFile.getFileType() == FileType::ONE_FILE)
    {
        writeInteraction(interactionFile.getFstream(), false);
    }
    
    P_->removeInteraction(this);
    I_->removeInteraction(this);
#if MERCURY_USE_MPI
    }
#endif
}


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_;
    os << " contactPoint " << contactPoint_;
    os << " force " << force_;
    os << " torque " << torque_;
    //\todo add information that can recreate the contact information (necessary for CG)
    //  os <<" timeStamp "<<timeStamp_<< " contactPoint " << contactPoint_ << " overlap " << overlap_ << " force " << force_ << " torque " << torque_;
}

void BaseInteraction::read(std::istream& is)
{
    //the rest gets read by the interaction handler
    std::string dummy;
    helpers::readOptionalVariable(is,"contactPoint",contactPoint_);
    is >> dummy >> force_;
    is >> 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(unsigned 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
    handler_->addObject(C);
}

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

void BaseInteraction::integrate(Mdouble timeStep UNUSED)
{

}

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

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

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

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

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

Vec3D BaseInteraction::getIP() const
{
    return getNormal() * getDistance();
}

Vec3D BaseInteraction::getIC() const
{
    return getNormal() * getDistance() + getContactPoint() - getP()->getPosition();
}

Vec3D BaseInteraction::getCP() const
{
    return getP()->getPosition() - getContactPoint();
}

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

    // do not write fstat output if the force is an internal bond
    if (isBonded()) return;

    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);
    
    if (PParticle != nullptr && !PParticle->isFixed())
    {
        os << time << " " << P_->getIndex()
           << " " << static_cast<int>((IParticle == nullptr ? (-I_->getIndex() - 1) : I_->getIndex()))
           << " " << getContactPoint()
           << " " << getOverlap()
           << " " << tangentialOverlap
           << " " << scalarNormalForce
           << " " << scalarTangentialForce
           << " " << (IParticle == nullptr ? -normal_ : normal_)
           << " " << (IParticle == nullptr ? -tangential : tangential) << std::endl;
    }
    if (IParticle != nullptr && !IParticle->isFixed() && IParticle->getPeriodicFromParticle() == nullptr)
    {
        os << time << " " << I_->getIndex()
           << " " << P_->getIndex()
           << " " << getContactPoint()
           << " " << 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()
{
}

void BaseInteraction::writeInteraction(std::ostream& os, bool created) const
{
    if (created)
    {
        os << "started ";
    }
    else
    {
        os << "ended   ";
    }
    
    if (dynamic_cast<BaseParticle*>(I_) != nullptr)
    {
        os << " particleIds     " << P_->getId() << " " << I_->getId() << " timeStamp ";
    }
    else
    {
        os << " particleWallIds " << P_->getId() << " " << I_->getId() << " timeStamp ";
    }
    
    if (created)
    {
        os << timeStamp_;
    }
    else
    {
        os << P_->getSpecies()->getHandler()->getDPMBase()->getTime();
    }
    
    os << std::endl;
    
}


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
{
    Mdouble invEffectiveRadius = getP()->getCurvature(contactPoint_) + getI()->getCurvature(contactPoint_);
    logger.assert(invEffectiveRadius>0,
                  "getEffectiveRadius(): interaction % at % has infinite effective radius",getId(), getContactPoint());
    return 1.0 / invEffectiveRadius;
}

Mdouble BaseInteraction::getEffectiveMass() const
{
    Mdouble invEffectiveMass = getP()->getInvMass() + getI()->getInvMass();
    logger.assert(invEffectiveMass>0,
            "getEffectiveMass(): interaction % at % has infinite effective mass",getId(), getContactPoint());
    return 1.0 / invEffectiveMass;
}

void BaseInteraction::actionsAfterTimeStep()
{
}

void BaseInteraction::gatherContactStatistics()
{
    auto* IParticle = dynamic_cast<BaseParticle*>(I_);
    auto* PParticle = static_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->isFixed())
    {
        getHandler()->getDPMBase()->gatherContactStatistics(
                P_->getIndex(),
                static_cast<int>((IParticle == nullptr ? (-I_->getIndex() - 1) : I_->getIndex())),
                centre,
                getOverlap(),
                tangentialOverlap,
                scalarNormalForce,
                scalarTangentialForce,
                (IParticle == nullptr ? -normal_ : normal_),
                (IParticle == nullptr ? -tangential : tangential));
    }
    if (IParticle != nullptr && !IParticle->isFixed() && IParticle->getPeriodicFromParticle() == nullptr)
    {
        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>();
}


void BaseInteraction::createMPIType()
{
}

void BaseInteraction::getMPIInteraction(void* historyDataArray, unsigned int index) const
{
}

void BaseInteraction::setMPIInteraction(void* historyDataArray, unsigned int index, const bool resetPointers)
{
}

void* BaseInteraction::createMPIInteractionDataArray(unsigned int numberOfInteractions) const
{
    logger(ERROR, "BaseInteraction::createMPIInteractionDataArray should never be called");
    void* historyArray;
    return historyArray;
}

void BaseInteraction::deleteMPIInteractionDataArray(void* dataArray)
{
    logger(WARN, "Why on earth is this function called?");
}

void
BaseInteraction::getInteractionDetails(void* interactionDataArray, unsigned int index, unsigned int& identificationP,
                                       unsigned int& identificationI, bool& isWallInteraction, unsigned& timeStamp)
{
    logger(ERROR, "Something went wrong, this function should not be called");
}


void BaseInteraction::setIdentificationP(unsigned int identification)
{
    identificationP_ = identification;
}

void BaseInteraction::setIdentificationI(int identification)
{
    identificationI_ = identification;
}

unsigned int BaseInteraction::getIdentificationP()
{
    return identificationP_;
}

int BaseInteraction::getIdentificationI()
{
    return identificationI_;
}

void BaseInteraction::setWallInteraction(bool flag)
{
    isWallInteraction_ = flag;
}

bool BaseInteraction::isWallInteraction()
{
    return isWallInteraction_;
}

void BaseInteraction::setBasicMPIInteractionValues(int P, int I, unsigned timeStamp, Vec3D force, Vec3D torque,
                                                   bool isWallInteraction, const bool resetPointers)
{
    this->setIdentificationP(P);
    this->setIdentificationI(I);
    this->setTimeStamp(timeStamp);
    this->setForce(force);
    this->setTorque(torque);
    this->setWallInteraction(isWallInteraction);
    if (resetPointers)
    {
        this->I_ = nullptr;
        this->P_ = nullptr;
    }
}

//note centre is unused
void BaseInteraction::setFStatData(std::fstream& fstat, BaseParticle* P, BaseParticle* I)
{
    Mdouble overlap, tangentialOverlap, scalarNormalForce, scalarTangentialForce;
    Vec3D centre, normal, tangential;
    fstat >> centre >> overlap >> tangentialOverlap >> scalarNormalForce >> scalarTangentialForce >> normal
          >> tangential;
    const Vec3D force = scalarNormalForce * normal + scalarTangentialForce * tangential;
    setForce(force);
    setNormal(normal);
    setOverlap(overlap);
    const Mdouble radius = P->getRadius();
    const Vec3D branch = (radius - 0.5 * getOverlap()) * getNormal();
    setContactPoint(P->getPosition() - branch);
    setDistance(radius + I->getRadius() - getOverlap());
}

void BaseInteraction::setFStatData(std::fstream& fstat, BaseParticle* P, BaseWall* I)
{
    Mdouble overlap, tangentialOverlap, scalarNormalForce, scalarTangentialForce;
    Vec3D centre, normal, tangential;
    fstat >> centre >> overlap >> tangentialOverlap >> scalarNormalForce >> scalarTangentialForce >> normal
          >> tangential;
    const Vec3D force = scalarNormalForce * normal + scalarTangentialForce * tangential;
    //note walls are defined different than particles (with an extra minus)
    setForce(-force);
    setNormal(-normal);
    setOverlap(overlap);
    const Mdouble radius = P->getRadius();
    const Vec3D branch = (radius - 0.5 * getOverlap()) * getNormal();
    setContactPoint(P->getPosition() - branch);
    setDistance(radius - getOverlap());
}