diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7205295fa7f36fd37d7837114cead443626e0293..0516669b5a1d54d2d4cc7970107ad29b4671194b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -956,6 +956,7 @@ if( (NOT DEFINED WALBERLA_BUILD_WITH_OPENMESH) OR WALBERLA_BUILD_WITH_OPENMESH )
add_definitions(-D_USE_MATH_DEFINES )
endif()
else()
+ message(" If OpenMesh required, set OPENMESH_DIR to the OpenMesh directory.")
set( WALBERLA_BUILD_WITH_OPENMESH OFF CACHE BOOL "Build with OpenMesh support" FORCE )
endif()
endif()
diff --git a/cmake/FindOpenMesh.cmake b/cmake/FindOpenMesh.cmake
index c9359b099e05242a2000eb54c8a152fbae947644..dfe98d2dd3758f50fd042bb3d265e79131447ff0 100644
--- a/cmake/FindOpenMesh.cmake
+++ b/cmake/FindOpenMesh.cmake
@@ -98,7 +98,7 @@ IF (NOT OPENMESH_FOUND)
"C:/libs/OpenMesh 3.0"
"C:/libs/OpenMesh 2.4.1"
"C:/libs/OpenMesh 2.4"
- "${OPENMESH_LIBRARY_DIR}"
+ "${OPENMESH_DIR}"
)
FIND_PATH (OPENMESH_INCLUDE_DIR OpenMesh/Core/Mesh/PolyMeshT.hh
diff --git a/src/mesa_pd/common/Contains.h b/src/mesa_pd/common/Contains.h
index 6f3bd1f1a4fa8b42da9fd69f7a1822ed634b7c63..a9e583440a6207f9d54f6ad811daea2eb62eaaa8 100644
--- a/src/mesa_pd/common/Contains.h
+++ b/src/mesa_pd/common/Contains.h
@@ -31,6 +31,7 @@
#include <mesa_pd/data/shape/Ellipsoid.h>
#include <mesa_pd/data/shape/HalfSpace.h>
#include <mesa_pd/data/shape/Sphere.h>
+#include <mesa_pd/data/shape/ConvexPolyhedron.h>
namespace walberla {
namespace mesa_pd {
@@ -75,6 +76,23 @@ bool isPointInsideCylindricalBoundary(const Vec3& point,
return distanceFromCylinderCenterLine.sqrLength() >= radius*radius;
}
+#ifdef WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE
+bool isPointInsideConvexPolyhedronBF(const Vec3& point, const mesh::TriangleMesh& mesh)
+{
+ return std::none_of(mesh.faces().begin(),
+ mesh.faces().end(),
+ [&](auto fh)
+ {
+ //check if point is on positive side of the face
+ const mesh::TriangleMesh::Normal &n = mesh.normal(fh); // Plane normal
+ const mesh::TriangleMesh::Point &pp = mesh.point(mesh.to_vertex_handle(mesh.halfedge_handle(fh))); // Point on plane
+
+ // normal * (p - pointOnPlane)
+ return (n[0] * (point[0] - pp[0]) + n[1] * (point[1] - pp[1]) + n[2] * (point[2] - pp[2]) >= real_t(0));
+ });
+}
+#endif
+
struct ContainsPointFunctor
{
template<typename ParticleAccessor_T, typename Shape_T>
@@ -123,6 +141,20 @@ struct ContainsPointFunctor
return isPointInsideCylindricalBoundary(point, ac.getPosition(particleIdx), cylindricalBoundary.getRadius(), cylindricalBoundary.getAxis());
}
+#ifdef WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE
+ template<typename ParticleAccessor_T>
+ bool operator()(const size_t particleIdx, const mesa_pd::data::ConvexPolyhedron& convexPolyhedron, 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");
+
+ auto pointBF = mesa_pd::transformPositionFromWFtoBF(particleIdx, ac, point);
+ if( pointBF.sqrLength() > convexPolyhedron.getBoundingSphereRadius() * convexPolyhedron.getBoundingSphereRadius() )
+ return false;
+
+ return isPointInsideConvexPolyhedronBF(pointBF, convexPolyhedron.getMesh());
+ }
+#endif
+
};
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index d2bf597b215594b57fa9f568cc9306135172c448..423dda6ef8a2811cd5225d8e34ed12752283c07a 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -231,6 +231,9 @@ if (WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE)
waLBerla_compile_test( NAME MESA_PD_VTK_ConvexPolyhedron FILES vtk/ConvexPolyhedronVTKOutput.cpp DEPENDS core mesa_pd mesh_common vtk )
waLBerla_execute_test( NAME MESA_PD_VTK_ConvexPolyhedron )
+ waLBerla_compile_test( NAME MESA_PD_Common_ContainsPointConvexPolyhedron FILES common/ContainsPointConvexPolyhedron.cpp DEPENDS core mesa_pd mesh_common )
+ waLBerla_execute_test( NAME MESA_PD_Common_ContainsPointConvexPolyhedron )
+
waLBerla_compile_test( NAME MESA_PD_Data_ConvexPolyhedron FILES data/ConvexPolyhedron.cpp DEPENDS core mesa_pd mesh_common )
waLBerla_execute_test( NAME MESA_PD_Data_ConvexPolyhedron )
diff --git a/tests/mesa_pd/common/ContainsPoint.cpp b/tests/mesa_pd/common/ContainsPoint.cpp
index be5b07996aca428761943d5c30ea0ad3adad6a5e..94a1a0ef806096754ec54c516ae44bfe2902ce46 100644
--- a/tests/mesa_pd/common/ContainsPoint.cpp
+++ b/tests/mesa_pd/common/ContainsPoint.cpp
@@ -41,7 +41,7 @@ using namespace walberla;
using mesa_pd::Vec3;
-/*!\brief Tests the contains pint functionality implemented in mesa_pd/common/Contains.h
+/*!\brief Tests the contains point functionality implemented in mesa_pd/common/Contains.h
*
* Currently the following shapes are tested:
* - sphere
diff --git a/tests/mesa_pd/common/ContainsPointConvexPolyhedron.cpp b/tests/mesa_pd/common/ContainsPointConvexPolyhedron.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8213d89a88eecef9e4784a6c12d671b0674eb264
--- /dev/null
+++ b/tests/mesa_pd/common/ContainsPointConvexPolyhedron.cpp
@@ -0,0 +1,122 @@
+//======================================================================================================================
+//
+// 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/kernel/SingleCast.h"
+
+#include <mesa_pd/data/shape/ConvexPolyhedron.h>
+#include <mesh_common/QHull.h>
+
+namespace contains_point_convex_polyhedron_test
+{
+using namespace walberla;
+using mesa_pd::Vec3;
+
+
+/*!\brief Tests the contains point functionality for ConvexPolyhedron shape, implemented in mesa_pd/common/Contains.h
+ * See tests/mesa_pd/common/ContainsPoint.cpp for other shapes
+ */
+
+//////////
+// 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;
+
+ /////////////////////////////////
+ // CONVEX POLYHEDRON ( = BOX ) //
+ /////////////////////////////////
+ {
+ Vec3 position(1_r, 0_r, 0_r);
+ Vec3 edgeLength{1_r};
+ Vec3 startVertex = - 0.5_r * edgeLength;
+
+ // manually construct triangle mesh of a cube
+ std::vector<Vec3> cornerVertices{startVertex,
+ startVertex + Vec3{edgeLength[0], 0, 0},
+ startVertex + Vec3{0, edgeLength[1], 0},
+ startVertex + Vec3{0, 0, edgeLength[2]},
+ startVertex + Vec3{edgeLength[0], edgeLength[1], 0},
+ startVertex + Vec3{edgeLength[0], 0, edgeLength[2]},
+ startVertex + Vec3{0, edgeLength[1], edgeLength[2]},
+ startVertex + Vec3{edgeLength[0], edgeLength[1], edgeLength[2]}};
+
+ mesh::QHull<mesh::TriangleMesh> qHull(cornerVertices);
+ qHull.run();
+
+ auto cpShape = shapeStorage->create<mesa_pd::data::ConvexPolyhedron>( qHull.mesh() );
+
+ mesa_pd::data::Particle&& p = *ps->create();
+ p.setPosition(position);
+ p.setShapeID(cpShape);
+ 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], "Convex Polyhedron 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], "Convex Polyhedron rotation check, case " << i << ": Wrong containment info for position " << testPositions[i] );
+ }
+ }
+
+ return 0;
+
+}
+
+} //namespace contains_point_convex_polyhedron_test
+
+int main( int argc, char **argv ){
+ contains_point_convex_polyhedron_test::main(argc, argv);
+}