Merge branch 'mesapd_polyhedra_for_coupling' into 'master'

OpenMesh improvements: more mesapd functionality, improved cmake

See merge request !541

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()

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

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

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 )
 

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

tests/mesa_pd/common/ContainsPointConvexPolyhedron.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/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);
+}