diff --git a/python/mesa_pd/templates/common/ParticleFunctions.templ.h b/python/mesa_pd/templates/common/ParticleFunctions.templ.h
index 807a1291949b644f6f17b637d008c19688cd9467..599c8d58b45939d9ced44fee0326531f55191800 100644
--- a/python/mesa_pd/templates/common/ParticleFunctions.templ.h
+++ b/python/mesa_pd/templates/common/ParticleFunctions.templ.h
@@ -27,10 +27,6 @@
#pragma once
#include <mesa_pd/data/DataTypes.h>
-#include <mesa_pd/domain/IDomain.h>
-
-#include <core/logging/Logging.h>
-#include <core/mpi/MPIManager.h>
namespace walberla {
namespace mesa_pd {
@@ -44,6 +40,33 @@ inline Vec3 getVelocityAtWFPoint(const size_t p_idx, Accessor& ac, const Vec3& w
return ac.getLinearVelocity(p_idx) + cross(ac.getAngularVelocity(p_idx), ( wf_pt - ac.getPosition(p_idx) ));
}
+/**
+ * Transformations between world frame (WF) and body frame (BF) coordinates
+ */
+template <typename Accessor>
+inline Vec3 transformPositionFromWFtoBF(const size_t p_idx, Accessor& ac, const Vec3& positionWF)
+{
+ return ac.getRotation(p_idx).getMatrix().getTranspose() * ( positionWF - ac.getPosition(p_idx) );
+}
+
+template <typename Accessor>
+inline Vec3 transformVectorFromWFtoBF(const size_t p_idx, Accessor& ac, const Vec3& vectorWF)
+{
+ return ac.getRotation(p_idx).getMatrix().getTranspose() * vectorWF;
+}
+
+template <typename Accessor>
+inline Vec3 transformPositionFromBFtoWF(const size_t p_idx, Accessor& ac, const Vec3& positionBF)
+{
+ return ac.getPosition(p_idx) + ac.getRotation(p_idx).getMatrix() * positionBF;
+}
+
+template <typename Accessor>
+inline Vec3 transformVectorFromBFtoWF(const size_t p_idx, Accessor& ac, const Vec3& vectorBF)
+{
+ return ac.getRotation(p_idx).getMatrix() * vectorBF;
+}
+
/**
* Force is applied at the center of mass.
*/
diff --git a/src/mesa_pd/common/Contains.h b/src/mesa_pd/common/Contains.h
index 2b13cd351a59c5d4b58605b7e69380dfec657832..6f3bd1f1a4fa8b42da9fd69f7a1822ed634b7c63 100644
--- a/src/mesa_pd/common/Contains.h
+++ b/src/mesa_pd/common/Contains.h
@@ -23,6 +23,12 @@
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/Flags.h>
#include <mesa_pd/data/IAccessor.h>
+
+#include <mesa_pd/common/ParticleFunctions.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>
@@ -30,22 +36,44 @@ namespace walberla {
namespace mesa_pd {
/*
- * contains functionality
+ * "contains point" functionality
+ * can either be formulated in world frame coordinates (then the rotation of the geometry is not taken into account)
+ * or in body frame coordinates (BF) which requires the point to be first transformed
*/
-bool isPointInsideSphere(const Vector3<real_t>& point,
- const Vector3<real_t>& spherePosition, const real_t sphereRadius )
+bool isPointInsideSphere(const Vec3& point,
+ const Vec3& spherePosition, const real_t sphereRadius )
{
return !((point - spherePosition).sqrLength() > sphereRadius * sphereRadius);
}
-bool isPointInsideHalfSpace(const Vector3<real_t>& point,
- const Vector3<real_t>& halfSpacePosition, const Vector3<real_t>& halfSpaceNormal )
+bool isPointInsideHalfSpace(const Vec3& point,
+ const Vec3& halfSpacePosition, const Vec3& halfSpaceNormal )
{
return !((point - halfSpacePosition) * halfSpaceNormal > real_t(0));
}
-//TODO add ellipsoids
+bool isPointInsideBoxBF(const Vec3& pointBF,
+ const Vec3& edgeLengths )
+{
+ return std::fabs(pointBF[0]) <= real_t(0.5)*edgeLengths[0] &&
+ std::fabs(pointBF[1]) <= real_t(0.5)*edgeLengths[1] &&
+ std::fabs(pointBF[2]) <= real_t(0.5)*edgeLengths[2];
+}
+
+bool isPointInsideEllipsoidBF(const Vec3& pointBF,
+ const Vec3& semiAxes )
+{
+ return ( (pointBF[0] * pointBF[0])/(semiAxes[0] * semiAxes[0]) + (pointBF[1] * pointBF[1])/(semiAxes[1] * semiAxes[1])
+ + (pointBF[2] * pointBF[2])/(semiAxes[2] * semiAxes[2]) <= 1_r );
+}
+
+bool isPointInsideCylindricalBoundary(const Vec3& point,
+ const Vec3& cylindricalBoundaryPosition, const real_t radius, const Vec3& axis )
+{
+ const Vec3 distanceFromCylinderCenterLine = (point - cylindricalBoundaryPosition) - ((point - cylindricalBoundaryPosition) * axis) * axis;
+ return distanceFromCylinderCenterLine.sqrLength() >= radius*radius;
+}
struct ContainsPointFunctor
{
@@ -71,6 +99,30 @@ struct ContainsPointFunctor
return isPointInsideHalfSpace(point, ac.getPosition(particleIdx), halfSpace.getNormal() );
}
+ template<typename ParticleAccessor_T>
+ bool operator()(const size_t particleIdx, const mesa_pd::data::Box& box, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
+ {
+ static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
+
+ return isPointInsideBoxBF(mesa_pd::transformPositionFromWFtoBF(particleIdx, ac, point), box.getEdgeLength());
+ }
+
+ template<typename ParticleAccessor_T>
+ bool operator()(const size_t particleIdx, const mesa_pd::data::Ellipsoid& ellipsoid, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
+ {
+ static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
+
+ return isPointInsideEllipsoidBF(mesa_pd::transformPositionFromWFtoBF(particleIdx, ac, point), ellipsoid.getSemiAxes());
+ }
+
+ template<typename ParticleAccessor_T>
+ bool operator()(const size_t particleIdx, const mesa_pd::data::CylindricalBoundary& cylindricalBoundary, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
+ {
+ static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
+
+ return isPointInsideCylindricalBoundary(point, ac.getPosition(particleIdx), cylindricalBoundary.getRadius(), cylindricalBoundary.getAxis());
+ }
+
};
diff --git a/src/mesa_pd/common/ParticleFunctions.h b/src/mesa_pd/common/ParticleFunctions.h
index bb1c89880afbce2932545b36709808b01db619b8..52b043d1289ae7fdd3778512f5b93f4e697a827c 100644
--- a/src/mesa_pd/common/ParticleFunctions.h
+++ b/src/mesa_pd/common/ParticleFunctions.h
@@ -27,10 +27,6 @@
#pragma once
#include <mesa_pd/data/DataTypes.h>
-#include <mesa_pd/domain/IDomain.h>
-
-#include <core/logging/Logging.h>
-#include <core/mpi/MPIManager.h>
namespace walberla {
namespace mesa_pd {
@@ -44,6 +40,33 @@ inline Vec3 getVelocityAtWFPoint(const size_t p_idx, Accessor& ac, const Vec3& w
return ac.getLinearVelocity(p_idx) + cross(ac.getAngularVelocity(p_idx), ( wf_pt - ac.getPosition(p_idx) ));
}
+/**
+ * Transformations between world frame (WF) and body frame (BF) coordinates
+ */
+template <typename Accessor>
+inline Vec3 transformPositionFromWFtoBF(const size_t p_idx, Accessor& ac, const Vec3& positionWF)
+{
+ return ac.getRotation(p_idx).getMatrix().getTranspose() * ( positionWF - ac.getPosition(p_idx) );
+}
+
+template <typename Accessor>
+inline Vec3 transformVectorFromWFtoBF(const size_t p_idx, Accessor& ac, const Vec3& vectorWF)
+{
+ return ac.getRotation(p_idx).getMatrix().getTranspose() * vectorWF;
+}
+
+template <typename Accessor>
+inline Vec3 transformPositionFromBFtoWF(const size_t p_idx, Accessor& ac, const Vec3& positionBF)
+{
+ return ac.getPosition(p_idx) + ac.getRotation(p_idx).getMatrix() * positionBF;
+}
+
+template <typename Accessor>
+inline Vec3 transformVectorFromBFtoWF(const size_t p_idx, Accessor& ac, const Vec3& vectorBF)
+{
+ return ac.getRotation(p_idx).getMatrix() * vectorBF;
+}
+
/**
* Force is applied at the center of mass.
*/
diff --git a/src/mesa_pd/common/RayParticleIntersection.h b/src/mesa_pd/common/RayParticleIntersection.h
index d1c9ce912dd325b931ecfa5e52021e0fe4ccbdc1..b62306aa5dd6014d301ddbd6943a8a0f45b700b9 100644
--- a/src/mesa_pd/common/RayParticleIntersection.h
+++ b/src/mesa_pd/common/RayParticleIntersection.h
@@ -24,6 +24,7 @@
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/Flags.h>
#include <mesa_pd/data/IAccessor.h>
+#include <mesa_pd/data/shape/Ellipsoid.h>
#include <mesa_pd/data/shape/HalfSpace.h>
#include <mesa_pd/data/shape/Sphere.h>
#include <mesa_pd/kernel/SingleCast.h>
@@ -33,10 +34,12 @@ namespace mesa_pd {
/*
* ray - particle intersection ratio functionality
+ * can either be formulated in world frame coordinates (then the rotation of the geometry is not taken into account)
+ * or in body frame coordinates (BF) which requires the point to be first transformed
*/
-real_t raySphereIntersectionRatio( const Vector3<real_t> & rayOrigin, const Vector3<real_t> & rayDirection,
- const Vector3<real_t> & spherePosition, const real_t sphereRadius )
+real_t raySphereIntersectionRatio( const Vec3& rayOrigin, const Vec3& rayDirection,
+ const Vec3& spherePosition, const real_t sphereRadius )
{
WALBERLA_ASSERT( !isPointInsideSphere(rayOrigin, spherePosition, sphereRadius ), "rayOrigin: " << rayOrigin );
WALBERLA_ASSERT( isPointInsideSphere(rayOrigin + rayDirection, spherePosition, sphereRadius ), "rayOrigin + rayDirection: " << rayOrigin + rayDirection );
@@ -53,14 +56,14 @@ real_t raySphereIntersectionRatio( const Vector3<real_t> & rayOrigin, const Vect
real_t delta = ( lsqr > rsqr ) ? s - q : s + q;
delta /= dirLength;
- WALBERLA_ASSERT_GREATER_EQUAL( delta, real_t( 0 ) );
- WALBERLA_ASSERT_LESS_EQUAL( delta, real_t( 1 ) );
+ WALBERLA_ASSERT_GREATER_EQUAL( delta, 0_r );
+ WALBERLA_ASSERT_LESS_EQUAL( delta, 1_r );
return delta;
}
-real_t rayHalfSpaceIntersectionRatio( const Vector3<real_t> & rayOrigin, const Vector3<real_t> & rayDirection,
- const Vector3<real_t> & halfSpacePosition, const Vector3<real_t> & halfSpaceNormal)
+real_t rayHalfSpaceIntersectionRatio( const Vec3& rayOrigin, const Vec3& rayDirection,
+ const Vec3& halfSpacePosition, const Vec3& halfSpaceNormal)
{
WALBERLA_ASSERT( !isPointInsideHalfSpace( rayOrigin, halfSpacePosition, halfSpaceNormal ), "rayOrigin: " << rayOrigin );
WALBERLA_ASSERT( isPointInsideHalfSpace( rayOrigin + rayDirection, halfSpacePosition, halfSpaceNormal ), "rayOrigin + rayDirection: " << rayOrigin + rayDirection );
@@ -69,21 +72,49 @@ real_t rayHalfSpaceIntersectionRatio( const Vector3<real_t> & rayOrigin, const V
auto diff = halfSpacePosition - rayOrigin;
- WALBERLA_ASSERT_FLOAT_UNEQUAL(denom, real_t(0));
+ WALBERLA_ASSERT_FLOAT_UNEQUAL(denom, 0_r);
real_t delta = diff * halfSpaceNormal / denom;
- WALBERLA_ASSERT_GREATER_EQUAL( delta, real_t( 0 ) );
- WALBERLA_ASSERT_LESS_EQUAL( delta, real_t( 1 ) );
+ WALBERLA_ASSERT_GREATER_EQUAL( delta, 0_r );
+ WALBERLA_ASSERT_LESS_EQUAL( delta, 1_r );
return delta;
}
-//TODO add ellipsoids from pe_coupling/geometry/PeIntersectionRatio.cpp
+real_t rayEllipsoidIntersectionRatioBF( const Vec3& rayOriginBF, const Vec3& rayDirectionBF,
+ const Vec3& ellipsoidSemiAxes)
+{
+ WALBERLA_ASSERT( !isPointInsideEllipsoidBF( rayOriginBF, ellipsoidSemiAxes ), "rayOriginBF: " << rayOriginBF );
+ WALBERLA_ASSERT( isPointInsideEllipsoidBF( rayOriginBF + rayDirectionBF, ellipsoidSemiAxes ), "rayOriginBF + rayDirectionBF: " << rayOriginBF + rayDirectionBF );
+
+ Matrix3<real_t> M = Matrix3<real_t>::makeDiagonalMatrix(1_r/ellipsoidSemiAxes[0], 1_r/ellipsoidSemiAxes[1], 1_r/ellipsoidSemiAxes[2]);
+
+ Vec3 P_M = M*rayOriginBF;
+ Vec3 d_M = M*rayDirectionBF;
+
+ const real_t a = d_M*d_M;
+ const real_t b = 2_r*P_M*d_M;
+ const real_t c = P_M*P_M - 1_r;
+
+ const real_t discriminant = b*b - 4_r*a*c;
+
+ WALBERLA_ASSERT_GREATER_EQUAL(discriminant, 0_r, "No intersection possible!");
+ WALBERLA_ASSERT_FLOAT_UNEQUAL(a, 0_r);
+
+ const real_t root = std::sqrt(discriminant);
+ real_t delta = (-b - root) / (2_r * a);
+
+ WALBERLA_ASSERT_GREATER_EQUAL( delta, 0_r );
+ WALBERLA_ASSERT_LESS_EQUAL( delta - 1_r, math::Limits<real_t>::accuracy(), delta );
+
+ return std::min(std::max(delta, real_c(0)), real_c(1));
+
+}
template <typename ParticleAccessor_T>
real_t intersectionRatioBisection( const size_t particleIdx, const ParticleAccessor_T& ac,
- const Vector3<real_t>& rayOrigin, const Vector3<real_t>& rayDirection,
+ const Vec3& rayOrigin, const Vec3& rayDirection,
real_t epsilon)
{
mesa_pd::kernel::SingleCast singleCast;
@@ -100,8 +131,8 @@ real_t intersectionRatioBisection( const size_t particleIdx, const ParticleAcces
while( qDelta * qDelta * sqDirectionLength >= sqEpsilon )
{
- qDelta *= real_t( 0.5 );
- Vector3<real_t> p = rayOrigin + q * rayDirection;
+ qDelta *= 0.5_r;
+ Vec3 p = rayOrigin + q * rayDirection;
if( singleCast(particleIdx, ac, containsPointFctr, ac, p) )
{
q -= qDelta;
@@ -112,8 +143,8 @@ real_t intersectionRatioBisection( const size_t particleIdx, const ParticleAcces
}
}
- WALBERLA_ASSERT_GREATER_EQUAL( q, real_t( 0 ) );
- WALBERLA_ASSERT_LESS_EQUAL( q, real_t( 1 ) );
+ WALBERLA_ASSERT_GREATER_EQUAL( q, 0_r );
+ WALBERLA_ASSERT_LESS_EQUAL( q, 1_r );
return q;
}
@@ -121,7 +152,7 @@ real_t intersectionRatioBisection( const size_t particleIdx, const ParticleAcces
struct RayParticleIntersectionRatioFunctor
{
template<typename ParticleAccessor_T, typename Shape_T>
- real_t operator()(const size_t particleIdx, const Shape_T& /*shape*/, const ParticleAccessor_T& ac, const Vector3<real_t>& rayOrigin, const Vector3<real_t>& rayDirection, real_t epsilon )
+ real_t operator()(const size_t particleIdx, const Shape_T& /*shape*/, const ParticleAccessor_T& ac, const Vec3& rayOrigin, const Vec3& rayDirection, real_t epsilon )
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
@@ -129,7 +160,7 @@ struct RayParticleIntersectionRatioFunctor
}
template<typename ParticleAccessor_T>
- real_t operator()(const size_t particleIdx, const mesa_pd::data::Sphere& sphere, const ParticleAccessor_T& ac, const Vector3<real_t>& rayOrigin, const Vector3<real_t>& rayDirection, real_t /*epsilon*/ )
+ real_t operator()(const size_t particleIdx, const mesa_pd::data::Sphere& sphere, const ParticleAccessor_T& ac, const Vec3& rayOrigin, const Vec3& rayDirection, real_t /*epsilon*/ )
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
@@ -137,13 +168,21 @@ struct RayParticleIntersectionRatioFunctor
}
template<typename ParticleAccessor_T>
- real_t operator()(const size_t particleIdx, const mesa_pd::data::HalfSpace& halfSpace, const ParticleAccessor_T& ac, const Vector3<real_t>& rayOrigin, const Vector3<real_t>& rayDirection, real_t /*epsilon*/ )
+ real_t operator()(const size_t particleIdx, const mesa_pd::data::HalfSpace& halfSpace, const ParticleAccessor_T& ac, const Vec3& rayOrigin, const Vec3& rayDirection, real_t /*epsilon*/ )
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
return rayHalfSpaceIntersectionRatio(rayOrigin, rayDirection, ac.getPosition(particleIdx), halfSpace.getNormal() );
}
+ template<typename ParticleAccessor_T>
+ real_t operator()(const size_t particleIdx, const mesa_pd::data::Ellipsoid& ellipsoid, const ParticleAccessor_T& ac, const Vec3& rayOrigin, const Vec3& rayDirection, real_t /*epsilon*/ )
+ {
+ static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
+
+ return rayEllipsoidIntersectionRatioBF(mesa_pd::transformPositionFromWFtoBF(particleIdx, ac, rayOrigin), mesa_pd::transformVectorFromWFtoBF(particleIdx, ac, rayDirection), ellipsoid.getSemiAxes() );
+ }
+
};
} //namespace mesa_pd
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index e55d03d2543789625502526388ebb716159e0e67..d2bf597b215594b57fa9f568cc9306135172c448 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -34,6 +34,9 @@ waLBerla_execute_test( NAME MESA_PD_COLLISIONDETECTION_SphereSupport )
waLBerla_compile_test( NAME MESA_PD_COMMON_IntersectionRatio FILES common/IntersectionRatio.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_COMMON_IntersectionRatio )
+waLBerla_compile_test( NAME MESA_PD_COMMON_ContainsPoint FILES common/ContainsPoint.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_COMMON_ContainsPoint )
+
waLBerla_compile_test( NAME MESA_PD_ContactDetection FILES ContactDetection.cpp DEPENDS blockforest core pe)
waLBerla_execute_test( NAME MESA_PD_ContactDetection PROCESSES 8 )
diff --git a/tests/mesa_pd/common/ContainsPoint.cpp b/tests/mesa_pd/common/ContainsPoint.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..be5b07996aca428761943d5c30ea0ad3adad6a5e
--- /dev/null
+++ b/tests/mesa_pd/common/ContainsPoint.cpp
@@ -0,0 +1,230 @@
+//======================================================================================================================
+//
+// 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
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/math/all.h"
+
+#include "mesa_pd/common/Contains.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/data/DataTypes.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 "mesa_pd/kernel/SingleCast.h"
+
+namespace contains_point_test
+{
+using namespace walberla;
+using mesa_pd::Vec3;
+
+
+/*!\brief Tests the contains pint functionality implemented in mesa_pd/common/Contains.h
+ *
+ * Currently the following shapes are tested:
+ * - sphere
+ * - halfspace
+ * - box (default and rotated)
+ * - ellipsoid (default and rotated)
+ * - cylindrical boundary
+ */
+
+//////////
+// MAIN //
+//////////
+int main( int argc, char **argv )
+{
+ debug::enterTestMode();
+
+ mpi::Environment env( argc, argv );
+
+ auto ps = std::make_shared<mesa_pd::data::ParticleStorage>(1);
+ auto shapeStorage = std::make_shared<mesa_pd::data::ShapeStorage>();
+ using ParticleAccessor = mesa_pd::data::ParticleAccessorWithShape;
+ ParticleAccessor accessor(ps, shapeStorage);
+
+
+ mesa_pd::kernel::SingleCast singleCast;
+ mesa_pd::ContainsPointFunctor containsPointFctr;
+
+ ////////////
+ // SPHERE //
+ ////////////
+ {
+ real_t sphereRadius = 1_r;
+ auto sphereShape = shapeStorage->create<mesa_pd::data::Sphere>( sphereRadius );
+
+ Vec3 position(1_r, 0_r, 0_r);
+
+ mesa_pd::data::Particle&& p = *ps->create();
+ p.setPosition(position);
+ p.setShapeID(sphereShape);
+ auto idx = p.getIdx();
+
+ std::vector<Vec3> testPositions {position, position + Vec3(0,0,1.01_r*sphereRadius)};
+ std::vector<bool> shouldBeContained {true, false};
+
+ for(size_t i = 0; i < testPositions.size(); ++i)
+ {
+ bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
+ WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Sphere check, case " << i << ": Wrong containment info for position " << testPositions[i] );
+ }
+ }
+
+ ///////////////
+ // HALFSPACE //
+ ///////////////
+ {
+ Vec3 position(1_r, 0_r, 0_r);
+ Vec3 normal(0_r, 1_r, 1_r);
+
+ auto planeShape = shapeStorage->create<mesa_pd::data::HalfSpace>( normal.getNormalized() );
+
+ mesa_pd::data::Particle&& p = *ps->create(true);
+ p.setPosition(position);
+ p.setShapeID(planeShape);
+ auto idx = p.getIdx();
+
+ std::vector<Vec3> testPositions {position - normal, position + normal};
+ std::vector<bool> shouldBeContained {true, false};
+
+ for(size_t i = 0; i < testPositions.size(); ++i)
+ {
+ bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
+ WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Halfspace check, case " << i << ": Wrong containment info for position " << testPositions[i] );
+ }
+ }
+
+ /////////
+ // BOX //
+ /////////
+ {
+ Vec3 position(1_r, 0_r, 0_r);
+ Vec3 edgeLength{1_r};
+
+ auto boxShape = shapeStorage->create<mesa_pd::data::Box>( edgeLength );
+
+ mesa_pd::data::Particle&& p = *ps->create();
+ p.setPosition(position);
+ p.setShapeID(boxShape);
+ auto idx = p.getIdx();
+
+ std::vector<Vec3> testPositions {position, position + 0.51_r*Vec3(0_r,0_r,edgeLength[2])};
+ std::vector<bool> shouldBeContained {true, false};
+
+ for(size_t i = 0; i < testPositions.size(); ++i)
+ {
+ bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
+ WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Box check, case " << i << ": Wrong containment info for position " << testPositions[i] );
+ }
+
+ auto rotation = p.getRotation();
+ rotation.rotate( Vec3(1_r,0_r,0_r), math::pi / 4_r ); // rotate by 45° around x axis
+ p.setRotation(rotation);
+
+ shouldBeContained[1] = true;
+
+ for(size_t i = 0; i < testPositions.size(); ++i)
+ {
+ bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
+ WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Box rotation check, case " << i << ": Wrong containment info for position " << testPositions[i] );
+ }
+ }
+
+ ///////////////
+ // ELLIPSOID //
+ ///////////////
+ {
+ Vec3 position(1_r, 0_r, 0_r);
+ Vec3 semiAxes{2_r,1_r, 1_r};
+
+ auto ellipsoidShape = shapeStorage->create<mesa_pd::data::Ellipsoid>( semiAxes );
+
+ mesa_pd::data::Particle&& p = *ps->create();
+ p.setPosition(position);
+ p.setShapeID(ellipsoidShape);
+ auto idx = p.getIdx();
+
+ std::vector<Vec3> testPositions {position,
+ position + 0.9_r*Vec3(semiAxes[0],0_r,0_r),
+ position + 1.1_r*Vec3(0_r,semiAxes[1], 0_r),
+ position + 1.1_r*Vec3(0_r,0_r,semiAxes[2])};
+ std::vector<bool> shouldBeContained {true, true, false, false};
+
+ for(size_t i = 0; i < testPositions.size(); ++i)
+ {
+ bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
+ WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Ellipsoid check, case " << i << ": Wrong containment info for position " << testPositions[i] );
+ }
+
+ auto rotation = p.getRotation();
+ rotation.rotate( Vec3(0_r,1_r,0_r), math::pi / 2_r ); // rotate by 90° around y axis
+ p.setRotation(rotation);
+
+ shouldBeContained[1] = false;
+ shouldBeContained[3] = true;
+
+ for(size_t i = 0; i < testPositions.size(); ++i)
+ {
+ bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
+ WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Ellipsoid rotation check, case " << i << ": Wrong containment info for position " << testPositions[i] );
+ }
+ }
+
+ //////////////////////////
+ // CYLINDRICAL BOUNDARY //
+ //////////////////////////
+ {
+ Vec3 position(1_r, 0_r, 0_r);
+ Vec3 axis(0_r, 1_r, 1_r);
+ real_t radius = 1_r;
+
+ auto cylindricalBoundaryShape = shapeStorage->create<mesa_pd::data::CylindricalBoundary>( radius, axis.getNormalized() );
+
+ mesa_pd::data::Particle&& p = *ps->create(true);
+ p.setPosition(position);
+ p.setShapeID(cylindricalBoundaryShape);
+ auto idx = p.getIdx();
+
+ Vec3 orthogonalAxis(0_r, -1_r, 1_r);
+ std::vector<Vec3> testPositions {position, position + axis, position + 1.1_r*radius*orthogonalAxis.getNormalized()};
+ std::vector<bool> shouldBeContained {false, false, true};
+
+ for(size_t i = 0; i < testPositions.size(); ++i)
+ {
+ bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
+ WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Cylindrical boundary check, case " << i << ": Wrong containment info for position " << testPositions[i] );
+ }
+ }
+
+ return 0;
+
+}
+
+} //namespace contains_point_test
+
+int main( int argc, char **argv ){
+ contains_point_test::main(argc, argv);
+}
diff --git a/tests/mesa_pd/common/IntersectionRatio.cpp b/tests/mesa_pd/common/IntersectionRatio.cpp
index 71db12ae58ce1e05ea246c0cc9a3e670d8752f99..edb34250e266717a1ea60fa1c65cd9a43d185cb2 100644
--- a/tests/mesa_pd/common/IntersectionRatio.cpp
+++ b/tests/mesa_pd/common/IntersectionRatio.cpp
@@ -24,36 +24,26 @@
#include "core/math/all.h"
#include "mesa_pd/common/RayParticleIntersection.h"
-#include "mesa_pd/data/ParticleAccessor.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
#include "mesa_pd/data/ParticleStorage.h"
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/data/DataTypes.h"
+#include <mesa_pd/data/shape/Ellipsoid.h>
+#include <mesa_pd/data/shape/HalfSpace.h>
+#include <mesa_pd/data/shape/Sphere.h>
#include "mesa_pd/kernel/SingleCast.h"
namespace intersection_ratio_test
{
using namespace walberla;
+using mesa_pd::Vec3;
-
-class ParticleAccessorWithShape : public mesa_pd::data::ParticleAccessor
-{
-public:
- ParticleAccessorWithShape(std::shared_ptr<mesa_pd::data::ParticleStorage>& ps, std::shared_ptr<mesa_pd::data::ShapeStorage>& ss)
- : ParticleAccessor(ps)
- , ss_(ss)
- {}
-
- mesa_pd::data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
- std::shared_ptr<mesa_pd::data::ShapeStorage> ss_;
-};
-
-
-/*!\brief Tests the ray-particle intersection ratio functionality of the RPDUtility.h in the lbm_rpd_coupling module
+/*!\brief Tests the ray-particle intersection ratio functionality implemented in mesa_pd/common/RayParticleIntersection.h
*
* Currently the following shapes are tested:
* - sphere
- * - halfspace ( default and rotated )
+ * - halfspace
+ * - ellipsoid (default and rotated)
*
* Additionally, the default variant with the bisection line search is tested with the help of a sphere.
*/
@@ -69,10 +59,10 @@ int main( int argc, char **argv )
auto ps = std::make_shared<mesa_pd::data::ParticleStorage>(1);
auto shapeStorage = std::make_shared<mesa_pd::data::ShapeStorage>();
- using ParticleAccessor = ParticleAccessorWithShape;
+ using ParticleAccessor = mesa_pd::data::ParticleAccessorWithShape;
ParticleAccessor accessor(ps, shapeStorage);
- const real_t epsilon( real_t(1e-4) );
+ const real_t epsilon( 1e-4_r );
mesa_pd::kernel::SingleCast singleCast;
mesa_pd::RayParticleIntersectionRatioFunctor intersectionRatioFctr;
@@ -81,33 +71,33 @@ int main( int argc, char **argv )
// SPHERE //
////////////
{
- real_t sphereRadius = real_t(1);
+ real_t sphereRadius = 1_r;
auto sphereShape = shapeStorage->create<mesa_pd::data::Sphere>( sphereRadius );
- Vector3<real_t> position(real_t(1), real_t(0), real_t(0));
+ Vec3 position(1_r, 0_r, 0_r);
mesa_pd::data::Particle&& p = *ps->create();
p.setPosition(position);
p.setShapeID(sphereShape);
auto idx = p.getIdx();
- Vector3<real_t> pos1(real_t(-0.5), real_t(0), real_t(0));
- Vector3<real_t> dir1(real_t(1), real_t(0), real_t(0));
+ Vec3 pos1(-0.5_r, 0_r, 0_r);
+ Vec3 dir1(1_r, 0_r, 0_r);
real_t delta1 = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos1, dir1, epsilon );
- WALBERLA_CHECK_FLOAT_EQUAL(delta1, real_t(0.5), "Intersection ratio 1 with sphere wrong!");
+ WALBERLA_CHECK_FLOAT_EQUAL(delta1, 0.5_r, "Intersection ratio 1 with sphere wrong!");
- Vector3<real_t> pos2(real_t(1), real_t(1), real_t(1));
- Vector3<real_t> dir2(real_t(0), -real_t(1), -real_t(1));
+ Vec3 pos2(1_r, 1_r, 1_r);
+ Vec3 dir2(0_r, -1_r, -1_r);
real_t delta2 = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos2, dir2, epsilon );
- WALBERLA_CHECK_FLOAT_EQUAL(delta2, (std::sqrt(2) - real_t(1)) / std::sqrt(2), "Intersection ratio 2 with sphere wrong!");
+ WALBERLA_CHECK_FLOAT_EQUAL(delta2, (std::sqrt(2_r) - 1_r) / std::sqrt(2_r), "Intersection ratio 2 with sphere wrong!");
}
///////////////
// HALFSPACE //
///////////////
{
- Vector3<real_t> position(real_t(1), real_t(0), real_t(0));
- Vector3<real_t> normal(real_t(0), real_t(1), real_t(1));
+ Vec3 position(1_r, 0_r, 0_r);
+ Vec3 normal(0_r, 1_r, 1_r);
auto planeShape = shapeStorage->create<mesa_pd::data::HalfSpace>( normal.getNormalized() );
@@ -116,25 +106,58 @@ int main( int argc, char **argv )
p.setShapeID(planeShape);
auto idx = p.getIdx();
- Vector3<real_t> pos1(real_t(1), real_t(0.5), real_t(0.5));
- Vector3<real_t> dir1(real_t(0), -real_t(1), -real_t(1));
+ Vec3 pos1(1_r, 0.5_r, 0.5_r);
+ Vec3 dir1(0_r, -1_r, -1_r);
real_t delta1 = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos1, dir1, epsilon );
- WALBERLA_CHECK_FLOAT_EQUAL(delta1, real_t(0.5), "Intersection ratio 1 with half space wrong!");
+ WALBERLA_CHECK_FLOAT_EQUAL(delta1, 0.5_r, "Intersection ratio 1 with half space wrong!");
- Vector3<real_t> dir2(real_t(0), real_t(0), -real_t(2));
+ Vec3 dir2(0_r, 0_r, -2_r);
real_t delta2 = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos1, dir2, epsilon );
- WALBERLA_CHECK_FLOAT_EQUAL(delta2, real_t(0.5), "Intersection ratio 2 with half space wrong!");
+ WALBERLA_CHECK_FLOAT_EQUAL(delta2, 0.5_r, "Intersection ratio 2 with half space wrong!");
- Vector3<real_t> dir3(real_t(0), -real_t(3), real_t(0));
+ Vec3 dir3(0_r, -3_r, 0_r);
real_t delta3 = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos1, dir3, epsilon );
- WALBERLA_CHECK_FLOAT_EQUAL(delta3, real_t(1)/real_t(3), "Intersection ratio 3 with half space wrong!");
+ WALBERLA_CHECK_FLOAT_EQUAL(delta3, 1_r/3_r, "Intersection ratio 3 with half space wrong!");
}
- /////////////////////////
- // HALFSPACE (rotated) //
- /////////////////////////
+ ///////////////
+ // ELLIPSOID //
+ ///////////////
+ {
+ Vec3 semiAxes{1_r, 1_r, 2_r};
+ auto ellipsoidShape = shapeStorage->create<mesa_pd::data::Ellipsoid>( semiAxes );
+
+ Vec3 position(1_r, 0_r, 0_r);
+
+ mesa_pd::data::Particle&& p = *ps->create();
+ p.setPosition(position);
+ p.setShapeID(ellipsoidShape);
+ auto idx = p.getIdx();
+
+ Vec3 pos1(-0.5_r, 0_r, 0_r);
+ Vec3 dir1(1_r, 0_r, 0_r);
+ real_t delta1 = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos1, dir1, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta1, 0.5_r, "Intersection ratio 1 with ellipsoid wrong!");
+
+ Vec3 pos2(1_r, 1.5_r, 0_r);
+ Vec3 dir2(0_r, -1_r, 0_r);
+ real_t delta2 = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos2, dir2, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta2, 0.5_r, "Intersection ratio 2 with ellipsoid wrong!");
+
+ auto rotation = p.getRotation();
+ rotation.rotate( Vec3(1_r,0_r,0_r), math::pi / 2_r ); // rotate by 90° around x axis
+ p.setRotation(rotation);
+
+ real_t delta1Rot = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos1, dir1, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta1Rot, 0.5_r, "Intersection ratio 1 with rotated ellipsoid wrong!");
+
+ Vec3 pos2Rot(1_r, 0_r, -1.5_r);
+ Vec3 dir2Rot(0_r, 0_r, 1_r);
+ real_t delta2Rot = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos2Rot, dir2Rot, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta2Rot, 0.5_r, "Intersection ratio 2 with rotated ellipsoid wrong!");
+
+ }
- // removed because rotation is not supported by half space, see half space docu
///////////////////////////
// Bisection Line Search //
@@ -142,25 +165,25 @@ int main( int argc, char **argv )
{
// note: here tested with a sphere and therefore called explicitly because otherwise the sphere specialization would be selected
- real_t sphereRadius = real_t(1);
+ real_t sphereRadius = 1_r;
auto sphereShape = shapeStorage->create<mesa_pd::data::Sphere>( sphereRadius );
- Vector3<real_t> position(real_t(1), real_t(0), real_t(0));
+ Vec3 position(1_r, 0_r, 0_r);
mesa_pd::data::Particle&& p = *ps->create();
p.setPosition(position);
p.setShapeID(sphereShape);
auto idx = p.getIdx();
- Vector3<real_t> pos1(real_t(-0.5), real_t(0), real_t(0));
- Vector3<real_t> dir1(real_t(1), real_t(0), real_t(0));
+ Vec3 pos1(-0.5_r, 0_r, 0_r);
+ Vec3 dir1(1_r, 0_r, 0_r);
real_t delta1 = mesa_pd::intersectionRatioBisection(idx, accessor, pos1, dir1, epsilon);
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(delta1, real_t(0.5), epsilon, "Intersection ratio 1 with bisection line search for sphere wrong!");
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(delta1, 0.5_r, epsilon, "Intersection ratio 1 with bisection line search for sphere wrong!");
- Vector3<real_t> pos2(real_t(1), real_t(1), real_t(1));
- Vector3<real_t> dir2(real_t(0), -real_t(1), -real_t(1));
+ Vec3 pos2(1_r, 1_r, 1_r);
+ Vec3 dir2(0_r, -1_r, -1_r);
real_t delta2 = mesa_pd::intersectionRatioBisection(idx, accessor, pos2, dir2, epsilon);
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(delta2, (std::sqrt(2) - real_t(1)) / std::sqrt(2), epsilon, "Intersection ratio 2 with bisection line search for sphere wrong!");
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(delta2, (std::sqrt(2_r) - 1_r) / std::sqrt(2_r), epsilon, "Intersection ratio 2 with bisection line search for sphere wrong!");
}
return 0;