Commit c654b324 authored by Sebastian Eibl's avatar Sebastian Eibl

added ellipsoid, gjk, epa

parent c47b6cb9
......@@ -33,13 +33,8 @@ public:
, ss_(ss)
{}
const walberla::real_t& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
walberla::real_t& getInvMassRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
void setInvMass(const size_t p_idx, const walberla::real_t& v) { ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass() = v;}
const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
auto& getInvInertiaBFRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
void setInvInertiaBF(const size_t p_idx, const Mat3& v) { ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF() = v;}
data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)].get();}
private:
......
......@@ -27,7 +27,7 @@ if __name__ == '__main__':
os.makedirs(args.path + "/src/mesa_pd/mpi/notifications", exist_ok = True)
os.makedirs(args.path + "/src/mesa_pd/vtk", exist_ok = True)
shapes = ["Sphere", "HalfSpace", "CylindricalBoundary", "Box"]
shapes = ["Sphere", "HalfSpace", "CylindricalBoundary", "Box", "Ellipsoid"]
ps = data.ParticleStorage()
ch = data.ContactHistory()
......
......@@ -351,7 +351,7 @@ void check_float_equal( const T & lhs, const U & rhs,
<< "File: " << filename << ":" << line << '\n'
<< "Expression: " << lhsExpression << " == " << rhsExpression << '\n'
//<< "ULP: " << distance << '\n'
<< "Values: " << std::setw(length) << std::setfill(' ') << lhsExpression << " = ";
<< "Values: " << std::setw(length) << std::setfill(' ') << lhsExpression << " = ";
printValue( ss, lhs ) << '\n';
......
......@@ -158,7 +158,7 @@ private:
* \brief Returns a random floating point number of type REAL in the range [min,max) (max excluded!)
*/
//**********************************************************************************************************************
template< typename REAL >
template< typename REAL = real_t >
REAL realRandom( const REAL min = REAL(0), const REAL max = REAL(1), std::mt19937 & generator = internal::getGenerator() )
{
static_assert( std::numeric_limits<REAL>::is_specialized && !std::numeric_limits<REAL>::is_integer, "Floating point type required/expected!" );
......
......@@ -84,6 +84,8 @@ inline real_t round( real_t a );
//
// The sign function only works for signed built-in data types. The attempt to use unsigned data
// types or user-defined class types will result in a compile time error.
//
// http://stackoverflow.com/questions/1903954/is-there-a-standard-sign-function-signum-sgn-in-c-c
*/
template< typename T >
inline const T sign( T a )
......
......@@ -46,6 +46,8 @@ namespace collision_detection {
class AnalyticContactDetection
{
public:
virtual ~AnalyticContactDetection() {}
size_t& getIdx1() {return idx1_;}
size_t& getIdx2() {return idx2_;}
Vec3& getContactPoint() {return contactPoint_;}
......@@ -116,7 +118,35 @@ public:
const data::Sphere& s,
Accessor& ac);
private:
template <typename Accessor>
bool operator()( const size_t idx1,
const size_t idx2,
const data::HalfSpace& p0,
const data::HalfSpace& p1,
Accessor& ac);
template <typename Accessor>
bool operator()( const size_t idx1,
const size_t idx2,
const data::CylindricalBoundary& p0,
const data::HalfSpace& p1,
Accessor& ac);
template <typename Accessor>
bool operator()( const size_t idx1,
const size_t idx2,
const data::HalfSpace& p0,
const data::CylindricalBoundary& p1,
Accessor& ac);
template <typename Accessor>
bool operator()( const size_t idx1,
const size_t idx2,
const data::CylindricalBoundary& p0,
const data::CylindricalBoundary& p1,
Accessor& ac);
protected:
size_t idx1_;
size_t idx2_;
Vec3 contactPoint_;
......@@ -254,6 +284,46 @@ inline bool AnalyticContactDetection::operator()( const size_t idx1,
return operator()(idx2, idx1, s, p, ac);
}
template <typename Accessor>
inline bool AnalyticContactDetection::operator()( const size_t /*idx1*/,
const size_t /*idx2*/,
const data::HalfSpace& /*geo1*/,
const data::HalfSpace& /*geo2*/,
Accessor& /*ac*/)
{
WALBERLA_ABORT("Collision between two half spaces is not defined!")
}
template <typename Accessor>
inline bool AnalyticContactDetection::operator()( const size_t /*idx1*/,
const size_t /*idx2*/,
const data::HalfSpace& /*geo1*/,
const data::CylindricalBoundary& /*geo2*/,
Accessor& /*ac*/)
{
WALBERLA_ABORT("Collision between half spaces and cylindrical boundary is not defined!")
}
template <typename Accessor>
inline bool AnalyticContactDetection::operator()( const size_t /*idx1*/,
const size_t /*idx2*/,
const data::CylindricalBoundary& /*geo1*/,
const data::HalfSpace& /*geo2*/,
Accessor& /*ac*/)
{
WALBERLA_ABORT("Collision between half spaces and cylindrical boundary is not defined!")
}
template <typename Accessor>
inline bool AnalyticContactDetection::operator()( const size_t /*idx1*/,
const size_t /*idx2*/,
const data::CylindricalBoundary& /*geo1*/,
const data::CylindricalBoundary& /*geo2*/,
Accessor& /*ac*/)
{
WALBERLA_ABORT("Collision between two cylindrical boundaries is not defined!")
}
inline
std::ostream& operator<<( std::ostream& os, const AnalyticContactDetection& ac )
{
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla is free software: you can
// redistribute it and/or modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation, either version 3 of
// the License, or (at your option) any later version.
//
// waLBerla is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License along
// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file GJK.h
//! \author Tobias Scharpff
//! \author Tobias Leemann
//
//======================================================================================================================
#pragma once
#include <mesa_pd/collision_detection/Support.h>
#include <mesa_pd/data/DataTypes.h>
#include <vector>
#include <core/Abort.h>
#include <core/math/Limits.h>
#include <core/math/Vector3.h>
namespace walberla {
namespace mesa_pd {
namespace collision_detection {
//*************************************************************************************************
/*!\brief Impelementation of the Gilbert-Johnson-Keerthi Algorithm.
*/
class GJK
{
public:
//**Constructor*********************************************************************************
/*! \name Constructor */
//@{
explicit GJK();
//@}
//**********************************************************************************************
//**Query functions*****************************************************************************
/*! \name Query functions */
//@{
real_t doGJK( const Support &geom1, const Support &geom2, Vec3& normal, Vec3& contactPoint );
bool doGJKmargin( const Support &geom1, const Support &geom2, const real_t margin = real_t(0));
//@}
//**********************************************************************************************
//**Get functions*******************************************************************************
/*! \name Get functions */
//@{
inline const std::vector<Vec3>& getSimplex() const { return simplex_; }
inline size_t getSimplexSize() const { return numPoints_; }
inline const std::vector<Vec3>& getSupportA() const { return supportA_; }
inline const std::vector<Vec3>& getSupportB() const { return supportB_; }
//@}
//**********************************************************************************************
private:
//**Utility functions***************************************************************************
/*! \name Utility functions */
//@{
bool simplex2(Vec3& d);
bool simplex3(Vec3& d);
bool simplex4(Vec3& d);
/// Checks if two vectors roughly point in the same direction
inline bool sameDirection ( const Vec3& vec1, const Vec3& vec2 ) const { return vec1 * vec2 > real_t(0.0); }
/// Checks if the length of a vector is zero or as close to zero that it can not be distinguished form zero
inline bool zeroLengthVector( const Vec3& vec ) const { return vec.sqrLength() < math::Limits<real_t>::fpuAccuracy(); }
real_t calcDistance ( Vec3& normal, Vec3& contactPoint );
inline const Vec3 putSupport(const Support &geom1,
const Support &geom2,
const Vec3& dir,
const real_t margin,
std::vector<Vec3> &simplex,
std::vector<Vec3> &supportA,
std::vector<Vec3> &supportB,
size_t index);
//@}
//**********************************************************************************************
//**Member variables****************************************************************************
/*! \name Member variables */
//@{
std::vector<Vec3> simplex_; ///< Container to hold the simplex.
std::vector<Vec3> supportA_; ///< Container to hold the support points generated in triangle mesh mA
std::vector<Vec3> supportB_; ///< Container to hold the support points generated in triangle mesh mB
unsigned char numPoints_; ///< Current number of points in the simplex.
Vec3 d_; ///< The next search direction.
//@}
//**********************************************************************************************
};
//*************************************************************************************************
} // namespace collision_detection
} // namespace mesa_pd
} // namespace walberla
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla is free software: you can
// redistribute it and/or modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation, either version 3 of
// the License, or (at your option) any later version.
//
// waLBerla is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License along
// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file GeneralContactDetection.h
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//
//======================================================================================================================
#pragma once
#include <mesa_pd/collision_detection/AnalyticContactDetection.h>
#include <mesa_pd/collision_detection/EPA.h>
#include <mesa_pd/collision_detection/GJK.h>
#include <mesa_pd/collision_detection/Support.h>
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/shape/BaseShape.h>
#include <mesa_pd/data/shape/Box.h>
#include <mesa_pd/data/shape/CylindricalBoundary.h>
#include <mesa_pd/data/shape/Ellipsoid.h>
#include <mesa_pd/data/shape/HalfSpace.h>
#include <mesa_pd/data/shape/Sphere.h>
#include <core/logging/Logging.h>
namespace walberla {
namespace mesa_pd {
namespace collision_detection {
/**
* Collision detection functor which uses recommended collision functions.
*
* Calculates and stores contact information between two particles.
* If a collision was successfully detected by the operator()
* contactPoint, contactNormal and penetrationDepth contain the
* contact information. If no collision was detected the values of
* these variables is undefined!
*/
class GeneralContactDetection : public AnalyticContactDetection
{
public:
using AnalyticContactDetection::operator();
template <typename GEO1_T, typename GEO2_T, typename Accessor>
bool operator()( const size_t idx1,
const size_t idx2,
const GEO1_T& geo1,
const GEO2_T& geo2,
Accessor& ac);
template <typename GEO2_T, typename Accessor>
bool operator()(const size_t idx1,
const size_t idx2,
const data::HalfSpace& geo1,
const GEO2_T& geo2,
Accessor& ac);
template <typename GEO1_T, typename Accessor>
bool operator()(const size_t idx1,
const size_t idx2,
const GEO1_T& geo1,
const data::HalfSpace& geo2,
Accessor& ac);
template <typename GEO2_T, typename Accessor>
bool operator()(const size_t idx1,
const size_t idx2,
const data::CylindricalBoundary& geo1,
const GEO2_T& geo2,
Accessor& ac);
template <typename GEO1_T, typename Accessor>
bool operator()(const size_t idx1,
const size_t idx2,
const GEO1_T& geo1,
const data::CylindricalBoundary& geo2,
Accessor& ac);
private:
bool collideGJKEPA(Support& geom0, Support& geom1);
};
template <typename GEO1_T, typename GEO2_T, typename Accessor>
bool GeneralContactDetection::operator()( const size_t idx1,
const size_t idx2,
const GEO1_T& geo1,
const GEO2_T& geo2,
Accessor& ac)
{
WALBERLA_ASSERT_UNEQUAL(idx1, idx2, "colliding with itself!");
//ensure collision order idx2 has to be larger than idx1
if (ac.getUid(idx2) < ac.getUid(idx1))
return operator()(idx2, idx1, geo2, geo1, ac);
getIdx1() = idx1;
getIdx2() = idx2;
Support e1(ac.getPosition(getIdx1()), ac.getRotation(getIdx1()), geo1);
Support e2(ac.getPosition(getIdx2()), ac.getRotation(getIdx2()), geo2);
return collideGJKEPA(e1, e2);
}
template <typename GEO2_T, typename Accessor>
bool GeneralContactDetection::operator()(const size_t idx1,
const size_t idx2,
const data::HalfSpace& geo1,
const GEO2_T& geo2,
Accessor& ac)
{
getIdx1() = idx1;
getIdx2() = idx2;
Support sup(ac.getPosition(idx2), ac.getRotation(idx2), geo2);
Vec3 support_dir = -geo1.getNormal();
// We now have a direction facing to the "wall".
// Compute support point of body b in this direction. This will be the furthest point overlapping.
Vec3 contactp = sup.support(support_dir);
real_t pdepth = contactp * geo1.getNormal() - ac.getPosition(idx1) * geo1.getNormal();
if(pdepth < contactThreshold_)
{ //We have a collision
contactNormal_ = support_dir;
penetrationDepth_ = pdepth;
contactPoint_ = contactp + real_t(0.5) * penetrationDepth_ * contactNormal_;
return true;
} else
{ //No collision
return false;
}
}
template <typename GEO1_T, typename Accessor>
bool GeneralContactDetection::operator()(const size_t idx1,
const size_t idx2,
const GEO1_T& geo1,
const data::HalfSpace& geo2,
Accessor& ac)
{
return operator()(idx2, idx1, geo2, geo1, ac);
}
template <typename GEO2_T, typename Accessor>
bool GeneralContactDetection::operator()(const size_t idx1,
const size_t idx2,
const data::CylindricalBoundary& geo1,
const GEO2_T& geo2,
Accessor& ac)
{
getIdx1() = idx1;
getIdx2() = idx2;
Support sup(ac.getPosition(idx2), ac.getRotation(idx2), geo2);
WALBERLA_CHECK_FLOAT_EQUAL(geo1.getAxis().sqrLength(), real_t(1));
auto d = ac.getPosition(idx2) - (dot((ac.getPosition(idx2) - ac.getPosition(idx1)), geo1.getAxis()) * geo1.getAxis() + ac.getPosition(idx1));
Vec3 farestP = sup.support(d);
auto d2 = farestP - (dot((farestP - ac.getPosition(idx1)), geo1.getAxis()) * geo1.getAxis() + ac.getPosition(idx1));
real_t dist = d2.sqrLength();
if(dist > geo1.getRadius() * geo1.getRadius())
{ //We have a collision
penetrationDepth_ = geo1.getRadius() - std::sqrt(dist);
normalize(d2);
contactNormal_ = d2;
contactPoint_ = farestP + d2 * penetrationDepth_ * real_t(0.5);
return true;
} else
{ //No collision
return false;
}
}
template <typename GEO1_T, typename Accessor>
bool GeneralContactDetection::operator()(const size_t idx1,
const size_t idx2,
const GEO1_T& geo1,
const data::CylindricalBoundary& geo2,
Accessor& ac)
{
return operator()(idx2, idx1, geo2, geo1, ac);
}
bool GeneralContactDetection::collideGJKEPA(Support& geom0, Support& geom1)
{
real_t margin = real_t(1e-4);
GJK gjk;
if(gjk.doGJKmargin(geom0, geom1, margin))
{
EPA epa;
epa.useSphereOptimization(false);
if (epa.doEPAmargin(geom0, geom1, gjk, contactNormal_, contactPoint_, penetrationDepth_, margin))
{
return true;
} else
{
return false;
}
} else
{
return false;
}
}
} //namespace collision_detection
} //namespace mesa_pd
} //namespace walberla
//======================================================================================================================
//
// This file is part of waLBerla. waLBerla is free software: you can
// redistribute it and/or modify it under the terms of the GNU General Public
// License as published by the Free Software Foundation, either version 3 of
// the License, or (at your option) any later version.
//
// waLBerla is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License along
// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file Support.h
//! \author Tobias Scharpff
//! \author Tobias Leemann
//
//======================================================================================================================
#pragma once
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/shape/BaseShape.h>
namespace walberla {
namespace mesa_pd {
namespace collision_detection {
class Support
{
public:
/// \attention Lifetime of shape has to surpass lifetime of Support!!!
/// No copy is stored.
Support(const Vec3& pos, const Rot3& rot, const data::BaseShape& shape)
: pos_(pos)
, rot_(rot)
, shape_(shape)
{}
/**
* \brief Estimates the point which is farthest in direction \a d.
*
* \param d The normalized search direction in world-frame coordinates.
* \return The support point in world-frame coordinates in direction a\ d.
*/
Vec3 support( Vec3 d ) const;
public:
Vec3 pos_;
Rot3 rot_;
const data::BaseShape& shape_;
};
inline
Vec3 Support::support( Vec3 d ) const
{
auto len = d.sqrLength();
if (math::equal(len, real_t(0)))
return Vec3(0,0,0);
d = rot_.getMatrix().getTranspose() * d.getNormalized(); //vectorFromWFtoBF
d = shape_.support(d);
return pos_ + rot_.getMatrix() * d; //vectorFromBFtoWF
}
} // namespace collision_detection
} // namespace mesa_pd
} // namespace walberla
......@@ -528,6 +528,8 @@ void swap(ParticleStorage::Particle lhs, ParticleStorage::Particle rhs)
std::swap(lhs.getNewContactHistoryRef(), rhs.getNewContactHistoryRef());
std::swap(lhs.getTemperatureRef(), rhs.getTemperatureRef());
std::swap(lhs.getHeatFluxRef(), rhs.getHeatFluxRef());
std::swap(lhs.getDvRef(), rhs.getDvRef());
std::swap(lhs.getDwRef(), rhs.getDwRef());
}
inline
......
......@@ -33,6 +33,7 @@
#include <mesa_pd/data/shape/HalfSpace.h>
#include <mesa_pd/data/shape/CylindricalBoundary.h>
#include <mesa_pd/data/shape/Box.h>
#include <mesa_pd/data/shape/Ellipsoid.h>
#include <core/Abort.h>
#include <core/debug/Debug.h>
......@@ -62,6 +63,7 @@ struct ShapeStorage
static_assert( Sphere::SHAPE_TYPE != HalfSpace::SHAPE_TYPE, "Shape types have to be different!" );
static_assert( Sphere::SHAPE_TYPE != CylindricalBoundary::SHAPE_TYPE, "Shape types have to be different!" );
static_assert( Sphere::SHAPE_TYPE != Box::SHAPE_TYPE, "Shape types have to be different!" );
static_assert( Sphere::SHAPE_TYPE != Ellipsoid::SHAPE_TYPE, "Shape types have to be different!" );
template <typename ShapeT, typename... Args>
size_t ShapeStorage::create(Args&&... args)
......@@ -81,6 +83,7 @@ ReturnType ShapeStorage::singleDispatch( ParticleStorage& ps, size_t idx, func&
case HalfSpace::SHAPE_TYPE : return f(ps, idx, *static_cast<HalfSpace*>(shapes[ps.getShapeID(idx)].get()));
case CylindricalBoundary::SHAPE_TYPE : return f(ps, idx, *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idx)].get()));
case Box::SHAPE_TYPE : return f(ps, idx, *static_cast<Box*>(shapes[ps.getShapeID(idx)].get()));
case Ellipsoid::SHAPE_TYPE : return f(ps, idx, *static_cast<Ellipsoid*>(shapes[ps.getShapeID(idx)].get()));
default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idx)]->getShapeType() << ") could not be determined!");
}
}
......@@ -116,6 +119,11 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
*static_cast<Sphere*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Box*>(shapes[ps.getShapeID(idy)].get()));
case Ellipsoid::SHAPE_TYPE : return f(ps,
idx,
*static_cast<Sphere*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
}
case HalfSpace::SHAPE_TYPE :
......@@ -141,6 +149,11 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
*static_cast<HalfSpace*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Box*>(shapes[ps.getShapeID(idy)].get()));
case Ellipsoid::SHAPE_TYPE : return f(ps,
idx,
*static_cast<HalfSpace*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
}
case CylindricalBoundary::SHAPE_TYPE :
......@@ -166,6 +179,11 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
*static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Box*>(shapes[ps.getShapeID(idy)].get()));
case Ellipsoid::SHAPE_TYPE : return f(ps,
idx,
*static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
}
case Box::SHAPE_TYPE :
......@@ -191,6 +209,41 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
*static_cast<Box*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Box*>(shapes[ps.getShapeID(idy)].get()));
case Ellipsoid::SHAPE_TYPE : return f(ps,
idx,
*static_cast<Box*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
}
case Ellipsoid::SHAPE_TYPE :
switch (shapes[ps.getShapeID(idy)]->getShapeType())
{
case Sphere::SHAPE_TYPE : return f(ps,
idx,
*static_cast<Ellipsoid*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Sphere*>(shapes[ps.getShapeID(idy)].get()));
case HalfSpace::SHAPE_TYPE : return f(ps,
idx,
*static_cast<Ellipsoid*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<HalfSpace*>(shapes[ps.getShapeID(idy)].get()));
case CylindricalBoundary::SHAPE_TYPE : return f(ps,
idx,
*static_cast<Ellipsoid*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idy)].get()));
case Box::SHAPE_TYPE : return f(ps,
idx,
*static_cast<Ellipsoid*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Box*>(shapes[ps.getShapeID(idy)].get()));
case Ellipsoid::SHAPE_TYPE : return f(ps,
idx,
*static_cast<Ellipsoid*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
}
default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idx)]->getShapeType() << ") could not be determined!");
......
......@@ -51,6 +51,8 @@ public:
const int& getShapeType() const {return shapeType_;}
virtual Vec3 support( const Vec3& /*d*/ ) const {WALBERLA_ABORT("Not implemented!");}
static const int INVALID_SHAPE = -1; ///< Unique *invalid* shape type identifier.\ingroup mesa_pd_shape
protected:
real_t mass_ = real_t(0); ///< mass
......
......@@ -41,7 +41,9 @@ public:
real_t getVolume() const override { return edgeLength_[0] * edgeLength_[1] * edgeLength_[2]; };
void updateMassAndInertia(const real_t density) override;
static const int SHAPE_TYPE = 3; ///< Unique shape type identifier for boxes.\ingroup mesa_pd_shape
Vec3 support( const Vec3& bfD ) const override;
constexpr static int SHAPE_TYPE = 3; ///< Unique shape type identifier for boxes.\ingroup mesa_pd_shape
private:
Vec3 edgeLength_; ///< edge length of the box
......@@ -56,8 +58,19 @@ void Box::updateMassAndInertia(const real_t density)
edgeLength_[0]*edgeLength_[0] + edgeLength_[2]*edgeLength_[2] ,
edgeLength_[0]*edgeLength_[0] + edgeLength_[1]*edgeLength_[1] ) * (m / static_cast<real_t>( 12 ));