Merge branch 'mr_extended_mesapd_shape_functionality' into 'master'

Added several new shape functionality for mesapd

See merge request walberla/walberla!526

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.
  */

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());
+   }
+
 };
 
 

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.
  */

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

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 )
 

tests/mesa_pd/common/ContainsPoint.cpp

+//======================================================================================================================
+//
+//  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);
+}

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;