diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index 949463d66847775c975ed0425ec70f7b159ff493..d356d40abd4a2b125008f67fbd1b4e248ae97a4d 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -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"]
+ shapes = ["Sphere", "HalfSpace", "CylindricalBoundary", "Box"]
ps = data.ParticleStorage()
ch = data.ContactHistory()
diff --git a/src/core/math/Rot3.h b/src/core/math/Rot3.h
index 71f8d629c6be292f6efd0ee2cfa7f71a86027ab2..b0ab7ecdbfd2c0a53a660b0009c59d70a92441e3 100644
--- a/src/core/math/Rot3.h
+++ b/src/core/math/Rot3.h
@@ -57,6 +57,7 @@ public:
void rotate(const Vector3<Type>& rot);
void rotate(const Vector3<Type>& axis, const real_t& angle);
+ void rotate(const Rot3<Type>& rot);
private:
Quaternion<Type> quat_;
Matrix3<Type> mat_;
@@ -115,6 +116,14 @@ inline void Rot3<Type>::rotate(const Vector3<Type>& axis, const real_t& angle)
}
}
+template< typename Type > // floating point type
+inline void Rot3<Type>::rotate(const Rot3<Type>& rot)
+{
+ quat_ = rot.quat_ * quat_;
+ mat_ = quat_.toRotationMatrix();
+ WALBERLA_ASSERT_FLOAT_EQUAL( mat_.getDeterminant(), real_t(1), "Corrupted rotation matrix determinant" );
+}
+
template< typename Type > // floating point type
inline std::ostream& operator<<( std::ostream& os, const Rot3<Type>& r )
{
diff --git a/src/mesa_pd/collision_detection/AnalyticCollisionFunctions.h b/src/mesa_pd/collision_detection/AnalyticCollisionFunctions.h
index 62f4f30e5bf8dcbb0b65de3176b16013ef676466..bcf3baf63b2d70d05a9992929cd3f52b4cac67af 100644
--- a/src/mesa_pd/collision_detection/AnalyticCollisionFunctions.h
+++ b/src/mesa_pd/collision_detection/AnalyticCollisionFunctions.h
@@ -34,9 +34,9 @@ bool detectSphereSphereCollision( const Vec3& pos1,
const real_t& radius1,
const Vec3& pos2,
const real_t& radius2,
- Vec3& contactPoint,
- Vec3& contactNormal,
- real_t& penetrationDepth,
+ Vec3& contactPoint,
+ Vec3& contactNormal,
+ real_t& penetrationDepth,
const real_t& contactThreshold)
{
contactNormal = pos1 - pos2;
@@ -53,14 +53,128 @@ bool detectSphereSphereCollision( const Vec3& pos1,
return false;
}
+inline
+bool detectSphereBoxCollision( const Vec3& pos1,
+ const real_t& radius1,
+ const Vec3& pos2,
+ const Vec3& edgeLength2,
+ const Rot3& rot2,
+ Vec3& contactPoint,
+ Vec3& contactNormal,
+ real_t& penetrationDepth,
+ const real_t& contactThreshold)
+{
+ //WALBERLA_ABORT(s << b);
+ const Vec3 l( real_t(0.5)*edgeLength2 ); // Box half side lengths
+ const Mat3& R( rot2.getMatrix() ); // Rotation of the box
+
+ bool outside( false );
+
+ // Calculating the distance between the sphere and the box
+ const Vec3 d( pos1 - pos2 );
+
+ // Calculating the the sphere/box distance in the box's frame of reference
+ Vec3 p( d[0]*R[0] + d[1]*R[3] + d[2]*R[6],
+ d[0]*R[1] + d[1]*R[4] + d[2]*R[7],
+ d[0]*R[2] + d[1]*R[5] + d[2]*R[8] );
+
+ // Projection of the x-component
+ if ( p[0] < -l[0] ) { p[0] = -l[0]; outside = true; }
+ else if( p[0] > l[0] ) { p[0] = l[0]; outside = true; }
+
+ // Projection of the y-component
+ if ( p[1] < -l[1] ) { p[1] = -l[1]; outside = true; }
+ else if( p[1] > l[1] ) { p[1] = l[1]; outside = true; }
+
+ // Projection of the z-component
+ if ( p[2] < -l[2] ) { p[2] = -l[2]; outside = true; }
+ else if( p[2] > l[2] ) { p[2] = l[2]; outside = true; }
+
+ // Special treatment if the sphere's center of mass is inside the box
+ // In this case, a contact at the sphere's center of mass with a normal away from
+ // the closest face of the box is generated.
+ if( !outside )
+ {
+ real_t dist( std::fabs(p[0]) - l[0] );
+ size_t face( 0 );
+
+ for( size_t i=1; i<3; ++i ) {
+ const real_t tmp( std::fabs(p[i]) - l[i] );
+ if( dist < tmp ) {
+ dist = tmp;
+ face = i;
+ }
+ }
+
+ // Calculating the normal direction of the contact
+ Vec3 n;
+ n[face] = math::sign( p[face] );
+ n = R * n;
+
+ // Calculating the contact distance
+ dist -= radius1;
+
+ // Creating a single contact between the sphere and the box
+ contactPoint = pos1;
+ contactNormal = n;
+ penetrationDepth = dist;
+ return true;
+ }
+
+ const Vec3 q( R * p ); // Transformation from the projection to the global world frame
+ const Vec3 n( d - q ); // Normal direction of the contact (pointing from the box to the sphere)
+
+ const real_t dist( n.length() - radius1 ); // Distance between the sphere and the box
+
+ // Creating a single contact between the sphere and the box
+ if( dist < contactThreshold )
+ {
+ contactPoint = pos2+q;
+ contactNormal = n.getNormalizedOrZero();
+ if (contactNormal.sqrLength() < real_t(0.1))
+ {
+ contactNormal = Vec3(1,0,0);
+ }
+ penetrationDepth = dist;
+ return true;
+ }
+ return false;
+}
+
+inline
+bool detectSphereCylindricalBoundaryCollision( const Vec3& pos1,
+ const real_t& radius1,
+ const Vec3& pos2,
+ const real_t& radius2,
+ const Vec3& axis2,
+ Vec3& contactPoint,
+ Vec3& contactNormal,
+ real_t& penetrationDepth,
+ const real_t& contactThreshold)
+{
+ WALBERLA_ASSERT_FLOAT_EQUAL( axis2.sqrLength(), real_t(1));
+ WALBERLA_ASSERT_GREATER( radius2, radius1 );
+
+ const Vec3 dist = (pos1 - pos2) - ((pos1 - pos2) * axis2) * axis2;
+ const real_t effRadius = radius2 - radius1;
+ if( effRadius * effRadius - dist.sqrLength() < contactThreshold ) {
+ contactNormal = -dist.getNormalized();
+ penetrationDepth = effRadius - dist.length();
+ contactPoint = ( pos1 - ( radius1 + penetrationDepth ) * contactNormal );
+
+ return true;
+ }
+ return false;
+}
+
inline
bool detectSphereHalfSpaceCollision( const Vec3& pos1,
const real_t& radius1,
const Vec3& pos2,
const Vec3& normal2,
- Vec3& contactPoint,
- Vec3& contactNormal,
- real_t& penetrationDepth,
+ Vec3& contactPoint,
+ Vec3& contactNormal,
+ real_t& penetrationDepth,
const real_t& contactThreshold)
{
/**
@@ -86,7 +200,7 @@ bool detectSphereHalfSpaceCollision( const Vec3& pos1,
}
-} //namespace collision_detection
} //namespace analytic
+} //namespace collision_detection
} //namespace mesa_pd
} //namespace walberla
diff --git a/src/mesa_pd/collision_detection/AnalyticContactDetection.h b/src/mesa_pd/collision_detection/AnalyticContactDetection.h
index dd5656bf571fcc53d8824d1dca58fdeadddc4a00..cbc5b34d4168d4368b5fb232cb6b2dccef4251ba 100644
--- a/src/mesa_pd/collision_detection/AnalyticContactDetection.h
+++ b/src/mesa_pd/collision_detection/AnalyticContactDetection.h
@@ -23,6 +23,8 @@
#include <mesa_pd/collision_detection/AnalyticCollisionFunctions.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/HalfSpace.h>
#include <mesa_pd/data/shape/Sphere.h>
@@ -72,6 +74,34 @@ public:
const data::Sphere& geo2,
Accessor& ac);
+ template <typename Accessor>
+ bool operator()( const size_t idx1,
+ const size_t idx2,
+ const data::Sphere& s,
+ const data::Box& bx,
+ Accessor& ac );
+
+ template <typename Accessor>
+ bool operator()( const size_t idx1,
+ const size_t idx2,
+ const data::Box& bx,
+ const data::Sphere& s,
+ Accessor& ac);
+
+ template <typename Accessor>
+ bool operator()( const size_t idx1,
+ const size_t idx2,
+ const data::Sphere& s,
+ const data::CylindricalBoundary& cb,
+ Accessor& ac );
+
+ template <typename Accessor>
+ bool operator()( const size_t idx1,
+ const size_t idx2,
+ const data::CylindricalBoundary& cb,
+ const data::Sphere& s,
+ Accessor& ac);
+
template <typename Accessor>
bool operator()( const size_t idx1,
const size_t idx2,
@@ -131,6 +161,68 @@ inline bool AnalyticContactDetection::operator()( const size_t idx1,
getContactThreshold());
}
+template <typename Accessor>
+inline bool AnalyticContactDetection::operator()( const size_t idx1,
+ const size_t idx2,
+ const data::Sphere& s,
+ const data::Box& bx,
+ Accessor& ac )
+{
+ getIdx1() = idx1;
+ getIdx2() = idx2;
+ return analytic::detectSphereBoxCollision(ac.getPosition(getIdx1()),
+ s.getRadius(),
+ ac.getPosition(getIdx2()),
+ bx.getEdgeLength(),
+ ac.getRotation(getIdx2()),
+ getContactPoint(),
+ getContactNormal(),
+ getPenetrationDepth(),
+ getContactThreshold());
+}
+
+template <typename Accessor>
+inline bool AnalyticContactDetection::operator()( const size_t idx1,
+ const size_t idx2,
+ const data::Box& bx,
+ const data::Sphere& s,
+ Accessor& ac)
+{
+ return operator()(idx2, idx1, s, bx, ac);
+}
+
+template <typename Accessor>
+inline bool AnalyticContactDetection::operator()( const size_t idx1,
+ const size_t idx2,
+ const data::Sphere& s,
+ const data::CylindricalBoundary& cb,
+ Accessor& ac )
+{
+ WALBERLA_ASSERT_UNEQUAL(idx1, idx2, "colliding with itself!");
+
+ getIdx1() = idx1;
+ getIdx2() = idx2;
+ return analytic::detectSphereCylindricalBoundaryCollision(ac.getPosition(getIdx1()),
+ s.getRadius(),
+ ac.getPosition(getIdx2()),
+ cb.getRadius(),
+ cb.getAxis(),
+ getContactPoint(),
+ getContactNormal(),
+ getPenetrationDepth(),
+ getContactThreshold());
+}
+
+template <typename Accessor>
+inline bool AnalyticContactDetection::operator()( const size_t idx1,
+ const size_t idx2,
+ const data::CylindricalBoundary& cb,
+ const data::Sphere& s,
+ Accessor& ac)
+{
+ return operator()(idx2, idx1, s, cb, ac);
+}
+
template <typename Accessor>
inline bool AnalyticContactDetection::operator()( const size_t idx1,
const size_t idx2,
diff --git a/src/mesa_pd/collision_detection/BroadPhase.h b/src/mesa_pd/collision_detection/BroadPhase.h
index 07fe49801d983c29e35988e8437b389b6de8387a..4bad17b24d100f93c28e79b0934cb7e4db090ef2 100644
--- a/src/mesa_pd/collision_detection/BroadPhase.h
+++ b/src/mesa_pd/collision_detection/BroadPhase.h
@@ -20,6 +20,8 @@
#pragma once
+#include <mesa_pd/data/Flags.h>
+
namespace walberla {
namespace mesa_pd {
namespace collision_detection {
diff --git a/src/mesa_pd/data/ShapeStorage.h b/src/mesa_pd/data/ShapeStorage.h
index 867098b379a89cc24bb46f75c87deeee4b6c0d81..4ab850ff7ac86d53bb05d0ff21082c2f34465466 100644
--- a/src/mesa_pd/data/ShapeStorage.h
+++ b/src/mesa_pd/data/ShapeStorage.h
@@ -31,6 +31,8 @@
#include <mesa_pd/data/shape/BaseShape.h>
#include <mesa_pd/data/shape/Sphere.h>
#include <mesa_pd/data/shape/HalfSpace.h>
+#include <mesa_pd/data/shape/CylindricalBoundary.h>
+#include <mesa_pd/data/shape/Box.h>
#include <core/Abort.h>
#include <core/debug/Debug.h>
@@ -58,6 +60,8 @@ struct ShapeStorage
};
//Make sure that no two different shapes have the same unique identifier!
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!" );
template <typename ShapeT, typename... Args>
size_t ShapeStorage::create(Args&&... args)
@@ -75,6 +79,8 @@ ReturnType ShapeStorage::singleDispatch( ParticleStorage& ps, size_t idx, func&
{
case Sphere::SHAPE_TYPE : return f(ps, idx, *static_cast<Sphere*>(shapes[ps.getShapeID(idx)].get()));
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()));
default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idx)]->getShapeType() << ") could not be determined!");
}
}
@@ -100,6 +106,16 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
*static_cast<Sphere*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<HalfSpace*>(shapes[ps.getShapeID(idy)].get()));
+ case CylindricalBoundary::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<Sphere*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idy)].get()));
+ case Box::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<Sphere*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<Box*>(shapes[ps.getShapeID(idy)].get()));
default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
}
case HalfSpace::SHAPE_TYPE :
@@ -115,6 +131,66 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
*static_cast<HalfSpace*>(shapes[ps.getShapeID(idx)].get()),
idy,
*static_cast<HalfSpace*>(shapes[ps.getShapeID(idy)].get()));
+ case CylindricalBoundary::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<HalfSpace*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idy)].get()));
+ case Box::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<HalfSpace*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<Box*>(shapes[ps.getShapeID(idy)].get()));
+ default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
+ }
+ case CylindricalBoundary::SHAPE_TYPE :
+ switch (shapes[ps.getShapeID(idy)]->getShapeType())
+ {
+ case Sphere::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<Sphere*>(shapes[ps.getShapeID(idy)].get()));
+ case HalfSpace::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<HalfSpace*>(shapes[ps.getShapeID(idy)].get()));
+ case CylindricalBoundary::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idy)].get()));
+ case Box::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<Box*>(shapes[ps.getShapeID(idy)].get()));
+ default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
+ }
+ case Box::SHAPE_TYPE :
+ switch (shapes[ps.getShapeID(idy)]->getShapeType())
+ {
+ case Sphere::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<Box*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<Sphere*>(shapes[ps.getShapeID(idy)].get()));
+ case HalfSpace::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<Box*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<HalfSpace*>(shapes[ps.getShapeID(idy)].get()));
+ case CylindricalBoundary::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<Box*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idy)].get()));
+ case Box::SHAPE_TYPE : return f(ps,
+ idx,
+ *static_cast<Box*>(shapes[ps.getShapeID(idx)].get()),
+ idy,
+ *static_cast<Box*>(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!");
diff --git a/src/mesa_pd/data/shape/Box.h b/src/mesa_pd/data/shape/Box.h
new file mode 100644
index 0000000000000000000000000000000000000000..7a3de3a9a761fffe61cc479dd47fa24df7239125
--- /dev/null
+++ b/src/mesa_pd/data/shape/Box.h
@@ -0,0 +1,65 @@
+//======================================================================================================================
+//
+// 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 Box.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/shape/BaseShape.h>
+
+#include <core/math/Constants.h>
+
+namespace walberla {
+namespace mesa_pd {
+namespace data {
+
+class Box : public BaseShape
+{
+public:
+ explicit Box(const Vec3& edgeLength)
+ : BaseShape(Box::SHAPE_TYPE)
+ , edgeLength_(edgeLength)
+ {}
+
+ const Vec3& getEdgeLength() const { return edgeLength_; }
+
+ 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
+
+private:
+ Vec3 edgeLength_; ///< edge length of the box
+};
+
+inline
+void Box::updateMassAndInertia(const real_t density)
+{
+ const real_t m = getVolume() * density;
+ const Mat3 I = Mat3::makeDiagonalMatrix(
+ edgeLength_[1]*edgeLength_[1] + edgeLength_[2]*edgeLength_[2] ,
+ edgeLength_[0]*edgeLength_[0] + edgeLength_[2]*edgeLength_[2] ,
+ edgeLength_[0]*edgeLength_[0] + edgeLength_[1]*edgeLength_[1] ) * (m / static_cast<real_t>( 12 ));
+
+ getInvMass() = real_t(1.0) / m;
+ getInvInertiaBF() = I.getInverse();
+}
+
+} //namespace data
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/src/mesa_pd/data/shape/CylindricalBoundary.h b/src/mesa_pd/data/shape/CylindricalBoundary.h
new file mode 100644
index 0000000000000000000000000000000000000000..dc2f3e67cc816b4ac275e2df7befcf80ec8ac624
--- /dev/null
+++ b/src/mesa_pd/data/shape/CylindricalBoundary.h
@@ -0,0 +1,64 @@
+//======================================================================================================================
+//
+// 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 CylindricalBoundary.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/shape/BaseShape.h>
+
+#include <core/math/Constants.h>
+
+namespace walberla {
+namespace mesa_pd {
+namespace data {
+
+class CylindricalBoundary : public BaseShape
+{
+public:
+ explicit CylindricalBoundary(const real_t& radius, const Vec3& axis)
+ : BaseShape(CylindricalBoundary::SHAPE_TYPE)
+ , radius_(radius)
+ , axis_(axis)
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL(axis.sqrLength(), real_t(1.0), "Axis has to be normalized!");
+ }
+
+ const real_t& getRadius() const { return radius_; }
+ const Vec3& getAxis() const { return axis_; }
+
+ real_t getVolume() const override { return std::numeric_limits<real_t>::infinity(); };
+ void updateMassAndInertia(const real_t density) override;
+
+ static const int SHAPE_TYPE = 2; ///< Unique shape type identifier for cylindrical boundaries.\ingroup mesa_pd_shape
+
+private:
+ real_t radius_; ///< radius of the cylinder
+ Vec3 axis_; ///< axis of the cylinder
+};
+
+inline
+void CylindricalBoundary::updateMassAndInertia(const real_t /*density*/)
+{
+ getInvMass() = real_t(0.0);
+ getInvInertiaBF() = Mat3(real_t(0.0));
+}
+
+} //namespace data
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/src/mesa_pd/kernel/DoubleCast.h b/src/mesa_pd/kernel/DoubleCast.h
index 4d0b09ab6dfde3f571cc194fb23bdb95f4e600d3..5b4cee0e212f0cdb259be97fea78b6b1ffc2312d 100644
--- a/src/mesa_pd/kernel/DoubleCast.h
+++ b/src/mesa_pd/kernel/DoubleCast.h
@@ -31,6 +31,8 @@
#include <mesa_pd/data/shape/BaseShape.h>
#include <mesa_pd/data/shape/Sphere.h>
#include <mesa_pd/data/shape/HalfSpace.h>
+#include <mesa_pd/data/shape/CylindricalBoundary.h>
+#include <mesa_pd/data/shape/Box.h>
#include <core/Abort.h>
#include <core/debug/Debug.h>
@@ -76,6 +78,16 @@ auto DoubleCast::operator()( size_t idx, size_t idy, Accessor& ac, func& f, Args
*static_cast<Sphere*>(ac.getShape(idx)),
*static_cast<HalfSpace*>(ac.getShape(idy)),
std::forward<Args>(args)...);
+ case CylindricalBoundary::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<Sphere*>(ac.getShape(idx)),
+ *static_cast<CylindricalBoundary*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
+ case Box::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<Sphere*>(ac.getShape(idx)),
+ *static_cast<Box*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
default : WALBERLA_ABORT("Shape type (" << ac.getShape(idy)->getShapeType() << ") could not be determined!");
}
case HalfSpace::SHAPE_TYPE :
@@ -91,6 +103,66 @@ auto DoubleCast::operator()( size_t idx, size_t idy, Accessor& ac, func& f, Args
*static_cast<HalfSpace*>(ac.getShape(idx)),
*static_cast<HalfSpace*>(ac.getShape(idy)),
std::forward<Args>(args)...);
+ case CylindricalBoundary::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<HalfSpace*>(ac.getShape(idx)),
+ *static_cast<CylindricalBoundary*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
+ case Box::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<HalfSpace*>(ac.getShape(idx)),
+ *static_cast<Box*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
+ default : WALBERLA_ABORT("Shape type (" << ac.getShape(idy)->getShapeType() << ") could not be determined!");
+ }
+ case CylindricalBoundary::SHAPE_TYPE :
+ switch (ac.getShape(idy)->getShapeType())
+ {
+ case Sphere::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<CylindricalBoundary*>(ac.getShape(idx)),
+ *static_cast<Sphere*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
+ case HalfSpace::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<CylindricalBoundary*>(ac.getShape(idx)),
+ *static_cast<HalfSpace*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
+ case CylindricalBoundary::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<CylindricalBoundary*>(ac.getShape(idx)),
+ *static_cast<CylindricalBoundary*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
+ case Box::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<CylindricalBoundary*>(ac.getShape(idx)),
+ *static_cast<Box*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
+ default : WALBERLA_ABORT("Shape type (" << ac.getShape(idy)->getShapeType() << ") could not be determined!");
+ }
+ case Box::SHAPE_TYPE :
+ switch (ac.getShape(idy)->getShapeType())
+ {
+ case Sphere::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<Box*>(ac.getShape(idx)),
+ *static_cast<Sphere*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
+ case HalfSpace::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<Box*>(ac.getShape(idx)),
+ *static_cast<HalfSpace*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
+ case CylindricalBoundary::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<Box*>(ac.getShape(idx)),
+ *static_cast<CylindricalBoundary*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
+ case Box::SHAPE_TYPE : return f(idx,
+ idy,
+ *static_cast<Box*>(ac.getShape(idx)),
+ *static_cast<Box*>(ac.getShape(idy)),
+ std::forward<Args>(args)...);
default : WALBERLA_ABORT("Shape type (" << ac.getShape(idy)->getShapeType() << ") could not be determined!");
}
default : WALBERLA_ABORT("Shape type (" << ac.getShape(idx)->getShapeType() << ") could not be determined!");
diff --git a/src/mesa_pd/kernel/SingleCast.h b/src/mesa_pd/kernel/SingleCast.h
index 0514d5a6519fcdc26754d42c08fbf18d49f1940f..316fbb62a847fb4fce90b181c263d8e53ff90cf1 100644
--- a/src/mesa_pd/kernel/SingleCast.h
+++ b/src/mesa_pd/kernel/SingleCast.h
@@ -31,6 +31,8 @@
#include <mesa_pd/data/shape/BaseShape.h>
#include <mesa_pd/data/shape/Sphere.h>
#include <mesa_pd/data/shape/HalfSpace.h>
+#include <mesa_pd/data/shape/CylindricalBoundary.h>
+#include <mesa_pd/data/shape/Box.h>
#include <core/Abort.h>
#include <core/debug/Debug.h>
@@ -64,6 +66,8 @@ auto SingleCast::operator()( size_t idx, Accessor& ac, func& f, Args&&... args )
{
case Sphere::SHAPE_TYPE : return f(idx, *static_cast<Sphere*>(ac.getShape(idx)), std::forward<Args>(args)...);
case HalfSpace::SHAPE_TYPE : return f(idx, *static_cast<HalfSpace*>(ac.getShape(idx)), std::forward<Args>(args)...);
+ case CylindricalBoundary::SHAPE_TYPE : return f(idx, *static_cast<CylindricalBoundary*>(ac.getShape(idx)), std::forward<Args>(args)...);
+ case Box::SHAPE_TYPE : return f(idx, *static_cast<Box*>(ac.getShape(idx)), std::forward<Args>(args)...);
default : WALBERLA_ABORT("Shape type (" << ac.getShape(idx)->getShapeType() << ") could not be determined!");
}
}
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index 1938034f1336141423ab5b677e64d5aeec8621a2..9abfc2b20472f2decf915335dc799689b9d5ecfd 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -4,6 +4,9 @@
#
###################################################################################################
+waLBerla_compile_test( NAME MESA_PD_COLLISIONDETECTION_AnalyticCollisionFunctions FILES collision_detection/AnalyticCollisionFunctions.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_COLLISIONDETECTION_AnalyticCollisionFunctions )
+
waLBerla_compile_test( NAME MESA_PD_COLLISIONDETECTION_AnalyticContactDetection FILES collision_detection/AnalyticContactDetection.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_COLLISIONDETECTION_AnalyticContactDetection )
diff --git a/tests/mesa_pd/collision_detection/AnalyticCollisionFunctions.cpp b/tests/mesa_pd/collision_detection/AnalyticCollisionFunctions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1dc914aaea5572e16e06f75cbfeea0430fd02a51
--- /dev/null
+++ b/tests/mesa_pd/collision_detection/AnalyticCollisionFunctions.cpp
@@ -0,0 +1,186 @@
+//======================================================================================================================
+//
+// 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 AnalyticCollisionFunctions.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <mesa_pd/collision_detection/AnalyticCollisionFunctions.h>
+#include <mesa_pd/data/DataTypes.h>
+
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+
+#include <cmath>
+#include <iostream>
+
+namespace walberla {
+namespace mesa_pd {
+
+void checkSphereSphereCollision( )
+{
+ using namespace walberla::mesa_pd::collision_detection::analytic;
+
+ const real_t radius = real_t(0.5);
+ auto dir = Vec3(real_t(1), real_t(2), real_t(3)).getNormalized();
+ real_t shift = real_c(0.75);
+
+ Vec3 contactPoint;
+ Vec3 contactNormal;
+ real_t penetrationDepth;
+
+ //check two spheres in contact
+ WALBERLA_CHECK( detectSphereSphereCollision( Vec3(real_t(0), real_t(0), real_t(0)),
+ radius,
+ dir * shift,
+ radius,
+ contactPoint,
+ contactNormal,
+ penetrationDepth,
+ real_t(0) ) );
+ WALBERLA_CHECK_FLOAT_EQUAL( contactPoint, dir * shift * real_t(0.5));
+ WALBERLA_CHECK_FLOAT_EQUAL( contactNormal, -dir );
+ WALBERLA_CHECK_FLOAT_EQUAL( penetrationDepth, shift - real_t(1) );
+
+ //no collision
+ WALBERLA_CHECK( !detectSphereSphereCollision( Vec3(real_t(0), real_t(0), real_t(0)),
+ radius,
+ dir * real_t(1.1),
+ radius,
+ contactPoint,
+ contactNormal,
+ penetrationDepth,
+ real_t(0) ) );
+}
+
+void checkSphereCylindricalBoundaryCollision( )
+{
+ using namespace walberla::mesa_pd::collision_detection::analytic;
+
+ const real_t radius = real_t(0.5);
+ auto pos = Vec3(real_t(6), real_t(-2), real_t(2));
+ auto dir = Vec3(real_t(1), real_t(1), real_t(1)).getNormalized();
+
+ Vec3 contactPoint;
+ Vec3 contactNormal;
+ real_t penetrationDepth;
+
+ //check two spheres in contact
+ WALBERLA_CHECK( detectSphereCylindricalBoundaryCollision( pos,
+ radius,
+ Vec3(real_t(0), real_t(0), real_t(0)),
+ real_t(5),
+ dir,
+ contactPoint,
+ contactNormal,
+ penetrationDepth,
+ real_t(0) ) );
+ auto nearestPoint = dot(pos, dir) * dir;
+ WALBERLA_CHECK_FLOAT_EQUAL( contactPoint, (pos - nearestPoint).getNormalized() * real_t(5) + nearestPoint);
+ WALBERLA_CHECK_FLOAT_EQUAL( contactNormal, (nearestPoint - pos).getNormalized() );
+ WALBERLA_CHECK_FLOAT_EQUAL( penetrationDepth, -((pos - nearestPoint).length() + radius - real_t(5)) );
+
+ //no collision
+ WALBERLA_CHECK( !detectSphereCylindricalBoundaryCollision( Vec3(real_t(0), real_t(0), real_t(0)),
+ radius,
+ Vec3(real_t(0), real_t(0), real_t(0)),
+ real_t(5),
+ Vec3(real_t(1), real_t(1), real_t(1)).getNormalized(),
+ contactPoint,
+ contactNormal,
+ penetrationDepth,
+ real_t(0) ) );
+}
+
+void checkSphereBoxCollision( )
+{
+ using namespace walberla::mesa_pd::collision_detection::analytic;
+
+ Vec3 contactPoint;
+ Vec3 contactNormal;
+ real_t penetrationDepth;
+
+ //check two spheres in contact
+ WALBERLA_CHECK( detectSphereBoxCollision( Vec3(real_t(2), real_t(2), real_t(2)),
+ real_c(std::sqrt(3)) + real_t(0.1),
+ Vec3(real_t(0), real_t(0), real_t(0)),
+ Vec3(real_t(2), real_t(2), real_t(2)),
+ Rot3(),
+ contactPoint,
+ contactNormal,
+ penetrationDepth,
+ real_t(0) ) );
+ WALBERLA_CHECK_FLOAT_EQUAL( contactPoint, Vec3(real_t(1), real_t(1), real_t(1)));
+ WALBERLA_CHECK_FLOAT_EQUAL( contactNormal, Vec3(real_t(2), real_t(2), real_t(2)).getNormalized() );
+ WALBERLA_CHECK_FLOAT_EQUAL( penetrationDepth, -real_t(0.1) );
+}
+
+void checkSphereHalfSpaceCollision( )
+{
+ using namespace walberla::mesa_pd::collision_detection::analytic;
+
+ const real_t radius = real_t(1.0);
+ auto dir = Vec3(real_t(1), real_t(2), real_t(3)).getNormalized();
+ real_t shift = real_c(0.75);
+
+ Vec3 contactPoint;
+ Vec3 contactNormal;
+ real_t penetrationDepth;
+
+ //check sphere - half space contact
+ WALBERLA_CHECK( detectSphereHalfSpaceCollision( Vec3(real_t(0), real_t(0), real_t(0)),
+ radius,
+ dir * shift,
+ -dir,
+ contactPoint,
+ contactNormal,
+ penetrationDepth,
+ real_t(0) ) );
+ WALBERLA_CHECK_FLOAT_EQUAL( contactPoint, dir * shift * real_t(1.0));
+ WALBERLA_CHECK_FLOAT_EQUAL( contactNormal, -dir );
+ WALBERLA_CHECK_FLOAT_EQUAL( penetrationDepth, shift - real_t(1) );
+
+
+ auto pos = Vec3(shift, real_t(0), real_t(0));
+ WALBERLA_CHECK( detectSphereHalfSpaceCollision( Vec3(real_t(0), real_t(0), real_t(0)),
+ radius,
+ pos,
+ -Vec3(real_t(1), real_t(0), real_t(0)),
+ contactPoint,
+ contactNormal,
+ penetrationDepth,
+ real_t(0) ) );
+ WALBERLA_CHECK_FLOAT_EQUAL( contactPoint, pos);
+ WALBERLA_CHECK_FLOAT_EQUAL( contactNormal, -Vec3(real_t(1), real_t(0), real_t(0)) );
+ WALBERLA_CHECK_FLOAT_EQUAL( penetrationDepth, pos[0] - real_t(1) );
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
+int main( int argc, char ** argv )
+{
+ walberla::Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+ walberla::mpi::MPIManager::instance()->useWorldComm();
+
+ walberla::mesa_pd::checkSphereSphereCollision();
+ walberla::mesa_pd::checkSphereCylindricalBoundaryCollision();
+ walberla::mesa_pd::checkSphereBoxCollision();
+ walberla::mesa_pd::checkSphereHalfSpaceCollision();
+
+ return EXIT_SUCCESS;
+}
diff --git a/tests/mesa_pd/kernel/DoubleCast.cpp b/tests/mesa_pd/kernel/DoubleCast.cpp
index ec429e4edf5a0aebdd30b5a8af89943138457ba1..4e2c11cdc050541f58d2109820e47996b7594b3b 100644
--- a/tests/mesa_pd/kernel/DoubleCast.cpp
+++ b/tests/mesa_pd/kernel/DoubleCast.cpp
@@ -48,6 +48,11 @@ public:
class ShapeTester
{
public:
+ void operator()(const size_t /*idx0*/, const size_t /*idx1*/, data::BaseShape& /*hs*/, data::BaseShape& /*hs*/)
+ {
+ WALBERLA_ABORT("Not implemented!");
+ }
+
void operator()(const size_t idx0, const size_t idx1, data::HalfSpace& /*hs*/, data::HalfSpace& /*hs*/)
{
WALBERLA_CHECK_EQUAL(idx0, 0);
diff --git a/tests/mesa_pd/kernel/SingleCast.cpp b/tests/mesa_pd/kernel/SingleCast.cpp
index c8ea556da8ce4f5bef07a993fffe47e5494b632c..cf137df8c68077d14a6c18f5dd464823ef383223 100644
--- a/tests/mesa_pd/kernel/SingleCast.cpp
+++ b/tests/mesa_pd/kernel/SingleCast.cpp
@@ -48,6 +48,11 @@ public:
class ShapeTester
{
public:
+ void operator()(const size_t /*idx*/, data::BaseShape& /*hs*/)
+ {
+ WALBERLA_ABORT("Not implemented!");
+ }
+
void operator()(const size_t idx, data::HalfSpace& /*hs*/)
{
WALBERLA_CHECK_EQUAL(idx, 0);