BaseWall.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 "BaseWall.h"
#include "DPMBase.h"
#include "WallHandler.h"
#include "Species/ParticleSpecies.h"

BaseWall::BaseWall()
{
    handler_ = nullptr;
    logger(DEBUG, "BaseWall::BaseWall() finished");
}

BaseWall::BaseWall(const BaseWall& w)
: BaseInteractable(w)
{
    handler_ = w.handler_;
    logger(DEBUG, "BaseWall::BaseWall(const BaseWall &p) finished");
}

BaseWall::~BaseWall()
{
    logger(DEBUG, "BaseWall::~BaseWall() finished");
}

void BaseWall::clear()
{
    logger(DEBUG, "BaseWall::clear(), this function shouldn't be called");
}

void BaseWall::read(std::istream& is)
{
    BaseInteractable::read(is);
}

void BaseWall::write(std::ostream& os) const
{
    BaseInteractable::write(os);
}

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

WallHandler* BaseWall::getHandler() const
{
    return handler_;
}

void BaseWall::setIndSpecies(unsigned int indSpecies)
{
    if (handler_ != nullptr)
    {
        setSpecies(getHandler()->getDPMBase()->speciesHandler.getObject(getIndSpecies()));
    }
    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 BaseWall::setSpecies(const ParticleSpecies* species)
{
    BaseInteractable::setSpecies(species);

    //set pointer to the handler, which is needed to retrieve species information
    if (getHandler() == nullptr)
    {
        SpeciesHandler* sH = species->getHandler();
        DPMBase* dB = sH->getDPMBase();
        if (dB != nullptr)
            setHandler(&dB->wallHandler);
    }
}

// returns the point intersecting a wall (x-p).n=0 and a line x=p0+t(p1-p0)
bool BaseWall::getLinePlaneIntersect (Vec3D& intersect, const Vec3D& p0, const Vec3D& p1, const Vec3D& n, const Vec3D& p)
{
    // t = (p-p0).n / (p1-p0).n
    //first compute the denominator
    Mdouble denominator = Vec3D::dot(p1-p0,n);
    if (fabs(denominator)>=1e-10)
    {
        Mdouble t = Vec3D::dot(p-p0,n)/denominator;
        if (t<1+1e-12&&t>-1e-12)
            intersect = p0 + t * (p1-p0);
            return true;
    } 
    return false;
}

//checks if point is in wall (if close to the wall, the point is assumed out of the wall)
bool BaseWall::isInsideWallVTK (const Vec3D& point, const Vec3D& normal, const Vec3D& position) const {
    return Vec3D::dot(position-point,normal)<-1e-12;
}

void BaseWall::projectOntoWallVTK (Vec3D& point0, const Vec3D& point1, const Vec3D& normal, const Vec3D& position) const {
    Vec3D dPoint = point1-point0;
    point0 += Vec3D::dot(position-point0,normal)/Vec3D::dot(dPoint,normal)*dPoint;
}

void BaseWall::intersectVTK(std::vector<Vec3D> &points, const Vec3D normal, const Vec3D position) const
{
    // find first point in Wall
    std::vector<Vec3D>::iterator firstIn=points.begin(); 
    for (;firstIn!=points.end(); firstIn++) {
        //stop if points[first] is in domain
        if (isInsideWallVTK(*firstIn,normal,position)) {
            break;
        }
    }

    //if all points are out of the wall
    if (firstIn==points.end()) {
        logger(DEBUG, "BaseWall::intersectVTK: all points out of wall");
        return;
    }

    // find first point out of the wall after firstIn
    std::vector<Vec3D>::iterator firstOut=firstIn+1; 
    for (;firstOut!=points.end(); firstOut++) {
        if (!isInsideWallVTK(*firstOut,normal,position)) {
            break;
        }
    }

    //if all points are in the wall
    if (firstOut==points.end()&&firstIn==points.begin()) {
        logger(DEBUG, "BaseWall::intersectVTK: points completely in wall; removing points");
        points.clear();
        return;
    }
    
    //if the sequence starts with a point out of the wall
    //Several cases have to be distinguished, multiple points in the wall the last point in or out of the wall: ooiiioo, ooioo, ooiii, or ooi
    //In addition, we add the case iiioo and ioo
    if (firstIn!=points.begin() || !isInsideWallVTK(points.back(),normal,position)) {
        // remove unnessesary points in the wall 
        // ooiiioo -> ooiioo, ooiii -> ooii, iiioo -> iioo
        if (firstOut-firstIn>2) {
            logger(DEBUG, "BaseWall::intersectVTK: remove unnessesary points in the wall");
            //necessary reset of firstOut, as erase invalidates iterators after erase point
            points.erase(firstIn+1,firstOut-1); //note: erase does not delete the last point
            firstOut = firstIn+2;
        }
        
        // if there is only one point in the wall, make it two
        // ooioo -> ooiioo, ooi -> ooii
        if (firstOut==firstIn+1) {
            logger(DEBUG, "BaseWall::intersectVTK: there is only one point in the wall, make it two");
            //necessary reset of firstIn, firstOut, as insert invalidates iterators
            unsigned in = firstIn-points.begin();
            points.insert(firstIn+1,*firstIn);
            firstIn = points.begin()+in;//necessary, unless capacity is set right
            firstOut = firstIn+2;
        }

        // three cases remain: ooiioo, ooii, iioo

        //move both points onto the surface of the wall
        if (firstIn!=points.begin()) {
            logger(DEBUG, "BaseWall::intersectVTK: move first point onto the surface of the wall");
            projectOntoWallVTK(*firstIn, *(firstIn-1), normal, position);
        } else {
            logger(DEBUG, "BaseWall::intersectVTK: move first point (at the beginning of the list) onto the surface of the wall");
            projectOntoWallVTK(points.front(), points.back(), normal, position);
        }
        
        if (firstOut!=points.end()) {
            logger(DEBUG, "BaseWall::intersectVTK: move second point onto the surface of the wall");
            projectOntoWallVTK(*(firstOut-1), *firstOut, normal, position);
        } else {
            logger(DEBUG, "BaseWall::intersectVTK: move second point (at the end of the list) onto the surface of the wall");
            projectOntoWallVTK(points.back(), points.front(), normal, position);
        }
    //if sequence starts and ends with a point in the wall: iiiooiii
    } else { 
        logger(DEBUG, "BaseWall::intersectVTK: sequence starts and ends with a point in the wall");
        
        // find first point in wall after firstOut
        for (firstIn=firstOut+1;firstIn!=points.end(); firstIn++) {
            if (isInsideWallVTK(*firstIn,normal,position)) {
                break;
            }
        }

        // remove unnessesary points in the wall 
        // iiiooiii -> iooi
        points.erase(firstIn+1,points.end());
        points.erase(points.begin(),firstOut-1); //note: erase does not delete the last point //note iterators are invalid now
        
        //move both points onto the surface of the wall: iooi
        projectOntoWallVTK(points.front(), *(points.begin()+1), normal, position);
        projectOntoWallVTK(points.back(), *(points.end()-2), normal, position);
    }
}

std::vector<BaseInteraction*> BaseWall::getInteractionWith(BaseParticle* p, Mdouble timeStamp, InteractionHandler* interactionHandler)
{
    Mdouble distance;
    Vec3D normal;
    std::vector<BaseInteraction*> interactions;
    if (getDistanceAndNormal(*p,distance,normal))
    {
        BaseInteraction* c = interactionHandler->getInteraction(p, this, timeStamp);
        c->setNormal(-normal);
        c->setDistance(distance);
        c->setOverlap(p->getRadius() - distance);
        c->setContactPoint(p->getPosition()-(p->getRadius() - 0.5 * c->getOverlap()) * c->getNormal());
        interactions.push_back(c);
    }
    return interactions;
}

void BaseWall::writeVTK (VTKContainer &vtk) const
{
    logger(WARN,"Wall % (%) cannot has no vtk writer defined",getIndex(),getName());
}

void BaseWall::addToVTK (const std::vector<Vec3D>& points, VTKContainer& vtk)
{
    if (points.size()!=0) {
        //all all values in myPoints to points
        vtk.points.insert(vtk.points.end(), points.begin(), points.end());

        // create one cell object containing all indices of the added points (created a triangle strip connecting these points)
        std::vector<double> cell;
        cell.reserve(vtk.points.size()+1);
        cell.push_back(vtk.points.size()-1);
        for (unsigned i=vtk.points.size()-points.size(); i<vtk.points.size(); i++) {
            cell.push_back(i);
        }

        //add this triangle strip to the vtk file
        vtk.triangleStrips.push_back(cell);
    }
}