AxisymmetricIntersectionOfWalls.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 "AxisymmetricIntersectionOfWalls.h"
#include "Particles/BaseParticle.h"
#include "WallHandler.h"
#include "DPMBase.h"

AxisymmetricIntersectionOfWalls::AxisymmetricIntersectionOfWalls()
{
    logger(DEBUG, "AxisymmetricIntersectionOfWalls() finished");
}

AxisymmetricIntersectionOfWalls::AxisymmetricIntersectionOfWalls(const AxisymmetricIntersectionOfWalls &other)
        : IntersectionOfWalls(other)
{
    logger(DEBUG, "AxisymmetricIntersectionOfWalls(const AxisymmetricIntersectionOfWalls &p) finished");
}

AxisymmetricIntersectionOfWalls::AxisymmetricIntersectionOfWalls(Vec3D position, Vec3D orientation, std::vector<AxisymmetricIntersectionOfWalls::normalAndPosition> walls, const ParticleSpecies* species)
    : IntersectionOfWalls(walls, species)
{
    setPosition(position);
    setOrientation(orientation);
}

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

AxisymmetricIntersectionOfWalls& AxisymmetricIntersectionOfWalls::operator=(const AxisymmetricIntersectionOfWalls& other)
{
    if (this == &other)
    {
        return *this;
    }
    else
    {
        return *(other.copy());
    }
}

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

bool AxisymmetricIntersectionOfWalls::getDistanceAndNormal(const BaseParticle& p, Mdouble& distance, Vec3D& normalReturn) const
{
    //transform to axisymmetric coordinates
    //move the coordinate system to the axis origin, so pOrigin=(xhat,yhat,zhat)
    Vec3D pOrigin = p.getPosition() -getPosition(); 
    Mdouble normal = Vec3D::dot(pOrigin, getOrientation());
    //tangential is the projection into the (xhat,yhat) plane
    Vec3D tangentialUnitVector = pOrigin - normal * getOrientation();
    Mdouble tangential = tangentialUnitVector.getLength();
    if (tangential!=0.0)
        tangentialUnitVector /= tangential;
    else //in this case the tangential vector is irrelevant
        logger(WARN, "Warning: Particle % is exactly on the symmetry axis of wall %", p.getIndex(), getIndex());
    Vec3D transformedPosition = Vec3D(tangential, 0.0, normal); //now P=(r,phi,zhat) is cylindrical
    Vec3D transformedNormal;
    //determine wall distance, normal and contact in axissymmetric coordinates
    //and transform from axisymmetric coordinates
    if (!IntersectionOfWalls::getDistanceAndNormal(transformedPosition, p.getWallInteractionRadius(), distance, transformedNormal))
    {
        //if not in contact
        return false;
    }
    else
    {
        //if in contact
        normalReturn = transformedNormal.Z * getOrientation() + transformedNormal.X * tangentialUnitVector;
        return true;
    }
}

void AxisymmetricIntersectionOfWalls::read(std::istream& is)
{
    IntersectionOfWalls::read(is);
}

void AxisymmetricIntersectionOfWalls::write(std::ostream& os) const
{
    IntersectionOfWalls::write(os);
}

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

void AxisymmetricIntersectionOfWalls::convertLimits (Vec3D& min, Vec3D& max) const {
        //define the box in the rz plane 
        double x[2] = {min.X, max.X};
        double y[2] = {min.Y, max.Y};
        double z[2] = {min.Z, max.Z};
        Vec3D o = getOrientation();
        min.X = 0;
        Mdouble maxXSquared = 0; //to be set
        min.Y = 0;
        max.Y = 0.001;
        min.Z = Vec3D::dot(max,o); //to be set
        max.Z = min.Z; //to be set
        for (unsigned i=0; i<2; i++)
        for (unsigned j=0; j<2; j++)
        for (unsigned k=0; k<2; k++) {
            Vec3D p = Vec3D(x[i],y[j],z[k])-getPosition();
            double Z = Vec3D::dot(p,o);
            double XSquared = Vec3D::getLengthSquared(p-Z*o);
            if (XSquared>maxXSquared) maxXSquared=XSquared;
            if (Z<min.Z) min.Z=Z;
            if (Z>max.Z) max.Z=Z;
        }
        max.X = sqrt(maxXSquared);
}

void AxisymmetricIntersectionOfWalls::writeVTK (VTKContainer& vtk) const
{
    for (auto wall=wallObjects_.begin(); wall!=wallObjects_.end(); wall++)
    {
        //plot each of the intersecting walls
        std::vector<Vec3D> myPoints;
        
        //first create a slice of non-rotated wall in the xz plane, 0<y<1            
        Vec3D min = getHandler()->getDPMBase()->getMin();
        Vec3D max = getHandler()->getDPMBase()->getMax();
        convertLimits(min, max);

        //create the basic slice for the first wall using the InfiniteWall routine
        wall->createVTK (myPoints,min,max);

        //create intersections with the other walls, similar to the IntersectionOfWalls routine
        for (auto other=wallObjects_.begin(); other!=wallObjects_.end(); other++)
        {
            if (other!=wall)
            {
                intersectVTK(myPoints, -other->getNormal(), other->getPosition());
            }
        }

        //only keep the y=0 values
        std::vector<Vec3D> rzVec;
        for (auto& p: myPoints)
        {   
            if (p.Y==0)
            {
                rzVec.push_back(p);
            }
        }
        if (rzVec.size()==0)
            return;
        
        //create 60 points on the unit circle
        unsigned nr = 60;
        struct XY {
            double X;
            double Y;
        };
        std::vector<XY> xyVec;
        for (unsigned ir=0; ir<nr; ir++) {
            Mdouble angle = 2.0 * constants::pi * ir/nr;
            xyVec.push_back({cos(angle),sin(angle)});
        }
        
        //now create rings of points on the axisym. shape
        unsigned nPoints = vtk.points.size();
        Vec3D p;
        for (auto rz : rzVec)
        {
            for (auto xy : xyVec)
            {
                p = getPosition();
                p.X += rz.X * xy.X;
                p.Y += rz.X * xy.Y;
                p.Z += rz.Z;
                vtk.points.push_back(p);
            }
        }

        //finally create the connectivity matri to plot shell-like triangle strips.
        unsigned nz = rzVec.size();
        unsigned nCells = vtk.triangleStrips.size();
        vtk.triangleStrips.reserve(nCells+(nz-1));
        for (unsigned iz=0; iz<nz-1; iz++) {
            std::vector<double> cell;
            cell.reserve(2*nr+2);
            for (unsigned ir=0; ir<nr; ir++) {
                cell.push_back(nPoints+ir+iz*nr);
                cell.push_back(nPoints+ir+(iz+1)*nr);
            }
            cell.push_back(nPoints+iz*nr);
            cell.push_back(nPoints+(iz+1)*nr);
            vtk.triangleStrips.push_back(cell);
        }
    }

}