diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt
index d834bd20ecc74452643579637bad66d3897e727c..58953c4bbedfbf38d0f883edbb882d2f2880228d 100644
--- a/apps/CMakeLists.txt
+++ b/apps/CMakeLists.txt
@@ -34,5 +34,4 @@ if ( WALBERLA_BUILD_WITH_PYTHON )
     add_subdirectory( pythonmodule )
     # no else with "EXLUDE_FROM_ALL" here, since if WALBERLA_BUILD_WITH_PYTHON_MODULE is not activated
     # waLBerla was build without -fPIC , so no linking into shared library is possible
-endif()
-
+endif()
\ No newline at end of file
diff --git a/apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp b/apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp
index 95678fab1449cf155d21063374f449c2291ee7b1..a9cca407c81d4ff61b905bb58abc754b880b51f5 100644
--- a/apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp
+++ b/apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp
@@ -24,53 +24,45 @@
 #include "blockforest/communication/UniformBufferedScheme.h"
 #include "blockforest/loadbalancing/StaticParMetis.h"
 
-#include "core/SharedFunctor.h"
 #include "core/Environment.h"
+#include "core/SharedFunctor.h"
 #include "core/logging/Logging.h"
 #include "core/math/IntegerFactorization.h"
+#include "core/timing/RemainingTimeLogger.h"
 
 #include "domain_decomposition/SharedSweep.h"
 
 #include "field/AddToStorage.h"
 #include "field/StabilityChecker.h"
 
+#include "geometry/InitBoundaryHandling.h"
+#include "geometry/mesh/TriangleMesh.h"
+#include "geometry/mesh/TriangleMeshIO.h"
+
+#include "lbm/BlockForestEvaluation.h"
+#include "lbm/PerformanceEvaluation.h"
+#include "lbm/PerformanceLogger.h"
 #include "lbm/boundary/factories/DefaultBoundaryHandling.h"
 #include "lbm/communication/PdfFieldPackInfo.h"
 #include "lbm/communication/SparsePdfFieldPackInfo.h"
 #include "lbm/field/AddToStorage.h"
+#include "lbm/lattice_model/CollisionModel.h"
 #include "lbm/lattice_model/D3Q19.h"
 #include "lbm/lattice_model/D3Q27.h"
-#include "lbm/lattice_model/CollisionModel.h"
 #include "lbm/lattice_model/ForceModel.h"
 #include "lbm/refinement/TimeStep.h"
 #include "lbm/sweeps/CellwiseSweep.h"
 #include "lbm/sweeps/SplitPureSweep.h"
 #include "lbm/vtk/VTKOutput.h"
-#include "lbm/BlockForestEvaluation.h"
-#include "lbm/PerformanceEvaluation.h"
-#include "lbm/PerformanceLogger.h"
 
-#include "geometry/mesh/TriangleMesh.h"
-#include "geometry/mesh/TriangleMeshIO.h"
-#include "geometry/InitBoundaryHandling.h"
-
-#include "mesh_common/TriangleMeshes.h"
-#include "mesh_common/MeshOperations.h"
-#include "mesh_common/DistanceComputations.h"
-#include "mesh_common/DistanceFunction.h"
-#include "mesh_common/MeshIO.h"
-#include "mesh_common/MatrixVectorOperations.h"
-#include "mesh_common/distance_octree/DistanceOctree.h"
-#include "mesh_common/vtk/VTKMeshWriter.h"
-#include "mesh_common/vtk/CommonDataSources.h"
+#include "mesh/blockforest/BlockExclusion.h"
 #include "mesh/blockforest/BlockForestInitialization.h"
 #include "mesh/blockforest/BlockWorkloadMemory.h"
-#include "mesh/blockforest/BlockExclusion.h"
 #include "mesh/blockforest/RefinementSelection.h"
-#include "mesh/boundary/BoundarySetup.h"
 #include "mesh/boundary/BoundaryInfo.h"
 #include "mesh/boundary/BoundaryLocation.h"
 #include "mesh/boundary/BoundaryLocationFunction.h"
+#include "mesh/boundary/BoundarySetup.h"
 #include "mesh/boundary/BoundaryUIDFaceDataSource.h"
 #include "mesh/boundary/ColorToBoundaryMapper.h"
 
@@ -78,11 +70,19 @@
 
 #include "timeloop/SweepTimeloop.h"
 
-#include "core/timing/RemainingTimeLogger.h"
-
 #include <cmath>
-#include <vector>
 #include <string>
+#include <vector>
+
+#include "mesh_common/DistanceComputations.h"
+#include "mesh_common/DistanceFunction.h"
+#include "mesh_common/MatrixVectorOperations.h"
+#include "mesh_common/MeshIO.h"
+#include "mesh_common/MeshOperations.h"
+#include "mesh_common/TriangleMeshes.h"
+#include "mesh_common/distance_octree/DistanceOctree.h"
+#include "mesh_common/vtk/CommonDataSources.h"
+#include "mesh_common/vtk/VTKMeshWriter.h"
 
 namespace walberla {
 
diff --git a/apps/showcases/CMakeLists.txt b/apps/showcases/CMakeLists.txt
index cda00720ff23502b5664f43fe43c140be2d026ee..cfdbf95914d9d7f94bd7070485a75b113002f1b0 100644
--- a/apps/showcases/CMakeLists.txt
+++ b/apps/showcases/CMakeLists.txt
@@ -2,6 +2,7 @@ add_subdirectory( BidisperseFluidizedBed )
 add_subdirectory( CombinedResolvedUnresolved )
 add_subdirectory( HeatConduction )
 add_subdirectory( Mixer )
+add_subdirectory( PegIntoSphereBed )
 if ( WALBERLA_BUILD_WITH_CODEGEN)
 add_subdirectory( PhaseFieldAllenCahn )
 endif()
diff --git a/apps/showcases/PegIntoSphereBed/CMakeLists.txt b/apps/showcases/PegIntoSphereBed/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3ddd2e024bbcf641da576604c9161f4eb79f528d
--- /dev/null
+++ b/apps/showcases/PegIntoSphereBed/CMakeLists.txt
@@ -0,0 +1,7 @@
+waLBerla_link_files_to_builddir( *.cfg )
+
+if (OPENMESH_FOUND AND WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE)
+    waLBerla_add_executable ( NAME MESA_PD_MESH_APP_PegIntoSphereBed
+            FILES PegIntoSphereBed.cpp
+            DEPENDS blockforest mesh_common mesa_pd core vtk )
+endif()
\ No newline at end of file
diff --git a/apps/showcases/PegIntoSphereBed/PegIntoSphereBed.cfg b/apps/showcases/PegIntoSphereBed/PegIntoSphereBed.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..a3e11a71c5007b46460bedeaaa7d4e29fe1f6a7d
--- /dev/null
+++ b/apps/showcases/PegIntoSphereBed/PegIntoSphereBed.cfg
@@ -0,0 +1,25 @@
+PegIntoSphereBed
+{
+   simulationCorner < 0, 0, 0 >;
+   simulationDomain < 2, 2, 3 >;
+   blocks < 1, 1, 1 >;
+   isPeriodic < 0, 0, 0 >;
+   initialRefinementLevel 1;
+   shift < 0.01, 0.01, 0.01 >;
+
+   sphereBedHeight 1.5;
+   sphereRadius 0.1;
+   sphereSpacing 0.2;
+   sphereDensity 5000;
+
+   pegBodyHeight 3;
+   pegPikeHeight 1;
+   pegRadius 0.5;
+   pegNumSideEdges 20;
+   pegPikeTipPosition < 1, 1, 1.5 >;
+   pegVelocity < 0, 0, -0.05 >;
+
+   dt                 0.0003;
+   simulationSteps    100000;
+   visSpacing         100;
+}
diff --git a/apps/showcases/PegIntoSphereBed/PegIntoSphereBed.cpp b/apps/showcases/PegIntoSphereBed/PegIntoSphereBed.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0a91f03c9021481f7f1fd983c031d27879c43e9d
--- /dev/null
+++ b/apps/showcases/PegIntoSphereBed/PegIntoSphereBed.cpp
@@ -0,0 +1,449 @@
+//======================================================================================================================
+//
+//  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   PegIntoSphereBed.cpp
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+/* This show case demonstrates a peg (consisting of a convex mesh) being rammed into
+ * a bed of spheres. It is purely based on MESA-PD and OpenMesh.
+ * The properties of the spheres and the peg, and additionally the height of the bed
+ * can be controlled by adjusting the appropriate parameters in the config file.
+ * The approximation of the peg's round body is controlled by the number of side
+ * edges. Its initial position is given with respect to its tip, the velocity with
+ * which it its afterwards moved is configurable as well.
+ *
+ *          peg radius
+ *           <------>
+ *            ______
+ *           |      |            ^
+ *           |      |            | peg
+ *           |      |            | body height
+ * oooooooooo|      |oooooooooo  v                   ^
+ * ooooooooooo\    /ooooooooooo  ^                   |
+ * oooooooooooo\  /oooooooooooo  | peg pike height   | sphere bed height
+ * ooooooooooooo\/ooooooooooooo  v                   |
+ * oooooooooooooooooooooooooooo                      |
+ * oooooooooooooooooooooooooooo                      v
+ * */
+
+#include "blockforest/Initialization.h"
+#include "blockforest/StructuredBlockForest.h"
+
+#include "mesa_pd/vtk/ConvexPolyhedron/data_sources/SurfaceVelocityVertexDataSource.h"
+
+#include "vtk/all.h"
+
+#include <OpenMesh/Core/IO/MeshIO.hh>
+#include <core/Environment.h>
+#include <core/grid_generator/SCIterator.h>
+#include <core/logging/Logging.h>
+#include <mesa_pd/collision_detection/GeneralContactDetection.h>
+#include <mesa_pd/common/ParticleFunctions.h>
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/Flags.h>
+#include <mesa_pd/data/LinkedCells.h>
+#include <mesa_pd/data/shape/HalfSpace.h>
+#include <mesa_pd/kernel/DoubleCast.h>
+#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
+#include <mesa_pd/kernel/SpringDashpot.h>
+#include <mesa_pd/kernel/SpringDashpotSpring.h>
+#include <mesa_pd/mpi/ContactFilter.h>
+#include <mesa_pd/mpi/ReduceContactHistory.h>
+#include <mesa_pd/mpi/ReduceProperty.h>
+#include <mesa_pd/mpi/notifications/ForceTorqueNotification.h>
+#include <mesa_pd/vtk/ParticleVtkOutput.h>
+#include <mesh_common/MeshOperations.h>
+#include <mesh_common/vtk/DistributedVTKMeshWriter.h>
+
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/data/shape/ConvexPolyhedron.h"
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/mpi/SyncNextNeighbors.h"
+#include "mesa_pd/vtk/ConvexPolyhedron/MeshParticleVTKOutput.h"
+#include "mesh_common/DistanceComputations.h"
+#include "mesh_common/MeshIO.h"
+#include "mesh_common/MeshOperations.h"
+#include "mesh_common/PolyMeshes.h"
+#include "mesh_common/TriangleMeshes.h"
+
+namespace walberla {
+
+void createPeg(mesh::TriangleMesh & mesh, mesa_pd::Vec3 & pegPikeTipPosition, real_t bodyHeight, real_t pikeHeight, real_t radius, uint_t numSideEdges) {
+   real_t alpha = real_t(2) * math::pi / real_t(numSideEdges); // 360° / numSideEdges -> approximation of cylinder and cone
+   real_t topCornerZ = pikeHeight + bodyHeight;
+   real_t bottomCornerZ = pikeHeight;
+
+   mesh::TriangleMesh::Point topCenterPoint(radius, radius, topCornerZ);
+   mesh::TriangleMesh::VertexHandle topCenterVertex = mesh.add_vertex(topCenterPoint);
+
+   mesh::TriangleMesh::Point bottomCenterPoint(radius, radius, real_t(0));
+   mesh::TriangleMesh::VertexHandle bottomCenterVertex = mesh.add_vertex(bottomCenterPoint);
+
+   mesh::TriangleMesh::VertexHandle firstTopVertex;
+   mesh::TriangleMesh::VertexHandle firstBottomVertex;
+
+   mesh::TriangleMesh::VertexHandle lastTopVertex;
+   mesh::TriangleMesh::VertexHandle lastBottomVertex;
+
+   for (uint_t e = 0; e < numSideEdges; ++e) {
+      real_t x_corner = radius + radius * std::sin(real_t(e) * alpha);
+      real_t y_corner = radius + radius * std::cos(real_t(e) * alpha);
+
+      mesh::TriangleMesh::Point topPoint(real_t(x_corner), real_t(y_corner), topCornerZ);
+      mesh::TriangleMesh::Point bottomPoint(real_t(x_corner), real_t(y_corner), bottomCornerZ);
+
+      mesh::TriangleMesh::VertexHandle newTopVertex = mesh.add_vertex(topPoint);
+      mesh::TriangleMesh::VertexHandle newBottomVertex = mesh.add_vertex(bottomPoint);
+
+      if (e > 0) {
+         // In the following, the order of vertices added to a face is very important!
+         // The half edge data structure needs to be valid, normals need to show into the correct direction.
+
+         // side faces (rectangle built up by two triangles)
+         mesh.add_face(newTopVertex, newBottomVertex, lastBottomVertex);
+         mesh.add_face(lastTopVertex, newTopVertex, lastBottomVertex);
+
+         // bottom - "pike" ("pizza slices")
+         mesh.add_face(newBottomVertex, bottomCenterVertex, lastBottomVertex);
+
+         // top face ("pizza slices")
+         mesh.add_face(topCenterVertex, newTopVertex, lastTopVertex);
+      } else {
+         firstTopVertex = newTopVertex;
+         firstBottomVertex = newBottomVertex;
+      }
+
+      lastTopVertex = newTopVertex;
+      lastBottomVertex = newBottomVertex;
+   }
+
+   // connect the first and the last sides
+   // side faces (rectangle built up by two triangles)
+   mesh.add_face(firstTopVertex, firstBottomVertex, lastBottomVertex);
+   mesh.add_face(lastTopVertex, firstTopVertex, lastBottomVertex);
+   // bottom - "pike"
+   mesh.add_face(firstBottomVertex, bottomCenterVertex, lastBottomVertex);
+   // top face
+   mesh.add_face(topCenterVertex, firstTopVertex, lastTopVertex);
+
+   // shift the mesh such that the centroid lies in 0,0,0 in its coordinate system (required)
+   // the pike tip is at radius,radius,0 up to now, after shifting the centroid by centroidShift to 0,0,0,
+   // the pike tip will be at -centroidShift.
+   auto centroidShift = mesh::toWalberla(mesh::computeCentroid(mesh));
+   mesh::translate(mesh, -centroidShift);
+
+   pegPikeTipPosition = -centroidShift + mesa_pd::Vec3(radius, radius, real_t(0));
+}
+
+mesa_pd::data::ParticleStorage::iterator createPlane( mesa_pd::data::ParticleStorage& ps,
+                                             mesa_pd::data::ShapeStorage& ss,
+                                             const Vector3<real_t>& pos,
+                                             const Vector3<real_t>& normal ) {
+   auto p0              = ps.create(true);
+   p0->getPositionRef() = pos;
+   p0->getShapeIDRef()  = ss.create<mesa_pd::data::HalfSpace>( normal );
+   p0->getOwnerRef()    = walberla::mpi::MPIManager::instance()->rank();
+   p0->getTypeRef()     = 0;
+   mesa_pd::data::particle_flags::set(p0->getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p0->getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+   mesa_pd::data::particle_flags::set(p0->getFlagsRef(), mesa_pd::data::particle_flags::NON_COMMUNICATING);
+   return p0;
+}
+
+class DEM
+{
+public:
+   DEM(const std::shared_ptr<mesa_pd::domain::BlockForestDomain>& domain, real_t dt, real_t mass)
+         : dt_(dt)
+         , domain_(domain)
+   {
+      sd_.setDampingT(0, 0, real_t(0));
+      sd_.setFriction(0, 0, real_t(0));
+      sd_.setParametersFromCOR(0, 0, real_t(0.9), dt*real_t(20), mass * real_t(0.5));
+
+      sds_.setParametersFromCOR(0, 0, real_t(0.9), dt*real_t(20), mass * real_t(0.5));
+      sds_.setCoefficientOfFriction(0,0,real_t(0.4));
+      sds_.setStiffnessT(0,0,real_t(0.9) * sds_.getStiffnessN(0,0));
+   }
+
+   inline
+   void operator()(const size_t idx1, const size_t idx2, mesa_pd::data::ParticleAccessorWithShape& ac)
+   {
+
+      ++contactsChecked_;
+      if (double_cast_(idx1, idx2, ac, gcd_, ac)) {
+         ++contactsDetected_;
+         if (contact_filter_(gcd_.getIdx1(), gcd_.getIdx2(), ac, gcd_.getContactPoint(),
+                             *domain_)) {
+            ++contactsTreated_;
+//            sd_(acd_.getIdx1(), acd_.getIdx2(), ac, acd_.getContactPoint(),
+//                 acd_.getContactNormal(), acd_.getPenetrationDepth());
+            sds_(gcd_.getIdx1(), gcd_.getIdx2(), ac, gcd_.getContactPoint(),
+                 gcd_.getContactNormal(), gcd_.getPenetrationDepth(), dt_);
+         }
+      }
+   }
+
+   inline void resetCounters() {
+      contactsChecked_ = 0;
+      contactsDetected_ = 0;
+      contactsTreated_ = 0;
+   }
+
+   inline int64_t getContactsChecked() const {
+      return contactsChecked_;
+   }
+
+   inline int64_t getContactsDetected() const {
+      return contactsDetected_;
+   }
+
+   inline int64_t getContactsTreated() const {
+      return contactsTreated_;
+   }
+
+private:
+   real_t dt_;
+   mesa_pd::kernel::DoubleCast double_cast_;
+   mesa_pd::mpi::ContactFilter contact_filter_;
+   std::shared_ptr<mesa_pd::domain::BlockForestDomain> domain_;
+   mesa_pd::kernel::SpringDashpot sd_ = mesa_pd::kernel::SpringDashpot(3);
+   mesa_pd::kernel::SpringDashpotSpring sds_ = mesa_pd::kernel::SpringDashpotSpring(3);
+   mesa_pd::collision_detection::GeneralContactDetection gcd_;
+   int64_t contactsChecked_ = 0;
+   int64_t contactsDetected_ = 0;
+   int64_t contactsTreated_ = 0;
+};
+
+void updatePegPosition(mesa_pd::data::ParticleAccessorWithShape & accessor, walberla::id_t pegUID, real_t dt) {
+   auto pegIdx = accessor.uidToIdx(pegUID);
+   WALBERLA_CHECK(pegIdx != accessor.getInvalidIdx());
+   auto newPegPosition = accessor.getPosition(pegIdx) + dt * accessor.getLinearVelocity(pegIdx);
+   accessor.setPosition(pegIdx, newPegPosition);
+}
+
+int main( int argc, char ** argv ) {
+   /// Setup
+   Environment env(argc, argv);
+
+   /// Config
+   auto cfg = env.config();
+   if (cfg == nullptr) WALBERLA_ABORT("No config specified!");
+   const Config::BlockHandle mainConf = cfg->getBlock( "PegIntoSphereBed" );
+
+   uint_t simulationSteps = mainConf.getParameter<uint_t>("simulationSteps", uint_t(1000));
+   uint_t visSpacing = mainConf.getParameter<uint_t>("visSpacing", uint_t(100));
+   real_t dt = mainConf.getParameter<real_t>("dt", real_t(0.0003));
+   Vector3<real_t> shift = mainConf.getParameter<Vector3<real_t>>("shift", Vector3<real_t>(real_t(0.01)));
+
+   real_t sphereBedHeight = mainConf.getParameter<real_t>("sphereBedHeight", real_t(1));
+   real_t sphereRadius = mainConf.getParameter<real_t>("sphereRadius", real_t(0.5));
+   real_t sphereSpacing = mainConf.getParameter<real_t>("sphereSpacing", real_t(0.5));
+   real_t sphereDensity = mainConf.getParameter<real_t>("sphereDensity", real_t(1000));
+
+   real_t pegBodyHeight = mainConf.getParameter<real_t>("pegBodyHeight", real_t(5));
+   real_t pegPikeHeight = mainConf.getParameter<real_t>("pegPikeHeight", real_t(2));
+   real_t pegRadius = mainConf.getParameter<real_t>("pegRadius", real_t(2));
+   uint_t pegNumSideEdges = mainConf.getParameter<uint_t>("pegNumSideEdges", uint_t(4));
+   Vector3<real_t> pegPikeTipPosition = mainConf.getParameter<Vector3<real_t>>("pegPikeTipPosition", Vector3<real_t>(real_t(1)));
+   Vector3<real_t> pegVelocity = mainConf.getParameter<Vector3<real_t>>("pegVelocity", Vector3<real_t>(0, 0, real_t(-0.05)));
+
+   std::string vtk_out = "vtk_out";
+
+   /// BlockForest
+   shared_ptr<BlockForest> forest = blockforest::createBlockForestFromConfig(mainConf);
+   auto domainSize = forest->getDomain();
+
+   /// MESAPD Domain
+   auto domain = std::make_shared<mesa_pd::domain::BlockForestDomain>(forest);
+
+   auto localDomain = forest->begin()->getAABB();
+   for (auto& blk : *forest) {
+      localDomain.merge(blk.getAABB());
+   }
+
+   /// MESAPD Data
+   auto ps = std::make_shared<mesa_pd::data::ParticleStorage>(4);
+   auto ss = std::make_shared<mesa_pd::data::ShapeStorage>();
+   mesa_pd::data::ParticleAccessorWithShape ac(ps, ss);
+   mesa_pd::data::LinkedCells lc(localDomain.getExtended(real_t(1)), real_t(2.1) * sphereRadius );
+   mesa_pd::mpi::SyncNextNeighbors SNN;
+
+   /// Mesh Peg
+   shared_ptr<mesh::TriangleMesh> pegMesh = make_shared<mesh::TriangleMesh>();
+   mesa_pd::Vec3 pegPikeTipOffset;
+   createPeg(*pegMesh, pegPikeTipOffset, pegBodyHeight, pegPikeHeight, pegRadius, pegNumSideEdges);
+   WALBERLA_LOG_INFO("Generated peg has " << pegMesh->n_vertices() << " vertices and " << pegMesh->n_faces() << " faces.");
+   /*if (!OpenMesh::IO::write_mesh(pegMesh, "peg_mesh.ply")) {
+      WALBERLA_ABORT("Error while writing peg mesh.");
+   } else {
+      WALBERLA_LOG_INFO("Wrote peg mesh to file.");
+   }*/
+   // We use simple distance calculation here due to a relatively low number of vertices. It it were much bigger, we would
+   // certainly be better off using a DistanceOctree for optimization.
+   auto pegTriDistance = make_shared<mesh::TriangleDistance<mesh::TriangleMesh>>(pegMesh);
+
+   /// MESAPD Particles
+   // Peg
+   auto pegShapeID = ss->create<mesa_pd::data::ConvexPolyhedron>(*pegMesh);
+   ss->shapes[pegShapeID]->updateMassAndInertia(real_t(1));
+   auto pegShape = dynamic_cast<mesa_pd::data::ConvexPolyhedron*>(ss->shapes[pegShapeID].get());
+
+   Vector3<real_t> pegPosition = pegPikeTipPosition - pegPikeTipOffset;
+   //WALBERLA_CHECK(domainSize.contains(pegPosition), "The pegs position " << pegPosition << " (defined at the mesh centroid) is outside the domain AABB!");
+   walberla::id_t pegUid = 0;
+
+   auto pegParticle = ps->create();
+
+   pegParticle->setShapeID(pegShapeID);
+   pegParticle->setPosition(pegPosition);
+   pegParticle->setOwner(walberla::MPIManager::instance()->rank());
+   pegParticle->setLinearVelocity(pegVelocity);
+   pegParticle->setInteractionRadius(pegShape->getBoundingSphereRadius());
+   mesa_pd::data::particle_flags::set(pegParticle->getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(pegParticle->getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+   mesa_pd::data::particle_flags::set(pegParticle->getFlagsRef(), mesa_pd::data::particle_flags::NON_COMMUNICATING);
+
+   pegUid = pegParticle->getUid();
+
+   // Spheres
+   auto sphereShapeId = ss->create<mesa_pd::data::Sphere>(sphereRadius);
+   ss->shapes[sphereShapeId]->updateMassAndInertia(real_t(sphereDensity));
+
+   auto extendedPegAABB = mesh::computeAABB(*pegMesh).getExtended(sphereRadius);
+   for (auto& iBlk : *forest) {
+      for (auto pt : grid_generator::SCGrid(iBlk.getAABB(),
+                                            Vector3<real_t>(sphereSpacing) * real_c(0.5) + shift,
+                                            sphereSpacing)) {
+         WALBERLA_CHECK(iBlk.getAABB().contains(pt));
+
+         if (pt[2] > sphereBedHeight) continue;
+         // check if a sphere at position pt would protrude into the peg
+         if (extendedPegAABB.contains(pt)) {
+            auto dSq = pegTriDistance->sqSignedDistance(mesh::toOpenMesh(pt - pegPosition));
+            if (dSq < sphereRadius*sphereRadius) {
+               continue;
+            }
+         }
+
+         auto sphereParticle = ps->create();
+
+         sphereParticle->setShapeID(sphereShapeId);
+         sphereParticle->setPosition(pt);
+         sphereParticle->setOwner(walberla::MPIManager::instance()->rank());
+         sphereParticle->setInteractionRadius(sphereRadius);
+      }
+   }
+   int64_t numParticles = int64_c(ps->size());
+   walberla::mpi::reduceInplace(numParticles, walberla::mpi::SUM);
+   WALBERLA_LOG_INFO_ON_ROOT("#particles created: " << numParticles);
+
+   // Confining Planes
+   auto planeShift = (sphereSpacing - sphereRadius - sphereRadius) * real_t(0.5);
+   auto confiningDomain = domainSize.getExtended(planeShift);
+   if (!forest->isPeriodic(0)) {
+      createPlane(*ps, *ss, confiningDomain.minCorner()+ shift, Vector3<real_t>(+1,0,0));
+      createPlane(*ps, *ss, confiningDomain.maxCorner()+ shift, Vector3<real_t>(-1,0,0));
+   }
+   if (!forest->isPeriodic(1)) {
+      createPlane(*ps, *ss, confiningDomain.minCorner()+ shift, Vector3<real_t>(0,+1,0));
+      createPlane(*ps, *ss, confiningDomain.maxCorner()+ shift, Vector3<real_t>(0,-1,0));
+   }
+   if (!forest->isPeriodic(2)) {
+      createPlane(*ps, *ss, confiningDomain.minCorner()+ shift, Vector3<real_t>(0,0,+1));
+      createPlane(*ps, *ss, confiningDomain.maxCorner()+ shift, Vector3<real_t>(0,0,-1));
+   }
+
+   SNN(*ps, *domain);
+
+   /// VTK Output
+   // domain output
+   auto vtkDomainOutput = vtk::createVTKOutput_DomainDecomposition(forest, "domain_decomposition",
+                                                                   uint_t(1), vtk_out, "simulation_step");
+   vtkDomainOutput->write();
+   // mesapd mesh output
+   mesa_pd::MeshParticleVTKOutput<mesh::PolyMesh> meshParticleVTK(ps, ss, "mesh", visSpacing);
+   meshParticleVTK.addFaceOutput<mesa_pd::data::SelectParticleUid>("Uid");
+   auto surfaceVelocityDataSource = make_shared<mesa_pd::SurfaceVelocityVertexDataSource<mesh::PolyMesh,
+         mesa_pd::data::ParticleAccessorWithShape>>("SurfaceVelocity", ac);
+   meshParticleVTK.addVertexDataSource(surfaceVelocityDataSource);
+   // mesapd particle output
+   auto particleVtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
+   particleVtkOutput->addOutput<mesa_pd::data::SelectParticleInteractionRadius>("interactionRadius");
+   particleVtkOutput->setParticleSelector([sphereShapeId](const mesa_pd::data::ParticleStorage::iterator& pIt){
+      return pIt->getShapeID() == sphereShapeId;
+   });
+   auto particleVtkWriter = walberla::vtk::createVTKOutput_PointData(particleVtkOutput, "particles", visSpacing,
+                                                                     vtk_out, "simulation_step");
+
+
+   /// MESAPD kernels
+   mesa_pd::kernel::ExplicitEuler explicitEuler(dt);
+   DEM dem(domain, dt, real_t(1));
+   mesa_pd::kernel::InsertParticleIntoLinkedCells ipilc;
+   mesa_pd::mpi::ReduceProperty RP;
+   mesa_pd::mpi::ReduceContactHistory RCH;
+
+   Vector3<real_t> globalAcceleration(real_t(0), real_t(0), real_t(-6));
+   auto addGravitationalForce = [&globalAcceleration](const size_t idx, mesa_pd::data::ParticleAccessorWithShape& ac_) {
+      auto mass = real_t(1) / ac_.getInvMass(idx);
+      auto force = mass * globalAcceleration;
+      mesa_pd::addForceAtomic(idx, ac_, force);
+   };
+
+   for (uint_t i = 0; i < simulationSteps; ++i) {
+      if(i % visSpacing == 0){
+         WALBERLA_LOG_INFO_ON_ROOT( "Timestep " << i << " / " << simulationSteps );
+      }
+
+      // VTK
+      meshParticleVTK();
+      particleVtkWriter->write();
+
+      // Prepare Data Structures
+      lc.clear();
+      ps->forEachParticle(true, mesa_pd::kernel::SelectAll(), ac, ipilc, ac, lc);
+
+      // Collision Resolution
+      dem.resetCounters();
+      lc.forEachParticlePairHalf(true, mesa_pd::kernel::SelectAll(), ac, dem, ac);
+      RP.operator()<mesa_pd::ForceTorqueNotification>(*ps);
+      RCH(*ps);
+
+      // Force Application
+      ps->forEachParticle(true, mesa_pd::kernel::SelectLocal(), ac, addGravitationalForce, ac);
+
+      // Integration
+      ps->forEachParticle(true, mesa_pd::kernel::SelectLocal(), ac, explicitEuler, ac);
+
+      updatePegPosition(ac, pegUid, dt);
+
+      SNN(*ps, *domain);
+   }
+
+   return EXIT_SUCCESS;
+}
+
+}
+
+int main( int argc, char* argv[] ) {
+   return walberla::main( argc, argv );
+}
\ No newline at end of file
diff --git a/apps/showcases/PegIntoSphereBed/showcase.png b/apps/showcases/PegIntoSphereBed/showcase.png
new file mode 100644
index 0000000000000000000000000000000000000000..48f19fe70c396b5c9f7d0e7b8a82edfc8bb8b3fb
Binary files /dev/null and b/apps/showcases/PegIntoSphereBed/showcase.png differ
diff --git a/cmake/waLBerlaFunctions.cmake b/cmake/waLBerlaFunctions.cmake
index 033d50814fc3906cd86baff5ad4df324957aac35..a5822aa849337b6550dc5d546899d85bd2978d9e 100644
--- a/cmake/waLBerlaFunctions.cmake
+++ b/cmake/waLBerlaFunctions.cmake
@@ -25,15 +25,18 @@ set ( WALBERLA_GLOB_FILES *.cpp
 #   [optional]           This is done using the ${arg}_FOUND variable.
 #                        Example: waLBerla_add_module( DEPENDS someModule BUILD_ONLY_IF_FOUND pe)
 #                                 This module is only built if PE_FOUND is true.
+#   OPTIONAL_DEPENDS     Lists modules, that this module might depend on. For example a module could depend on mesh_common if OpenMesh is
+#   [optional]           available.
 #
 #######################################################################################################################
 
 function ( waLBerla_add_module )
     set( options )
     set( oneValueArgs )
-    set( multiValueArgs DEPENDS EXCLUDE FILES BUILD_ONLY_IF_FOUND )
+    set( multiValueArgs DEPENDS EXCLUDE FILES BUILD_ONLY_IF_FOUND OPTIONAL_DEPENDS )
     cmake_parse_arguments( ARG "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN} )
 
+    set( ALL_DEPENDENCIES ${ARG_DEPENDS} ${ARG_OPTIONAL_DEPENDS})
     # Module name is the directory relative to WALBERLA_MODULE_DIRS
     get_current_module_name ( moduleName )
     get_module_library_name ( moduleLibraryName ${moduleName} )
@@ -76,7 +79,7 @@ function ( waLBerla_add_module )
     endif ( )
 
     # Dependency Check
-    check_dependencies( missingDeps additionalDeps FILES ${sourceFiles} EXPECTED_DEPS ${ARG_DEPENDS} ${moduleName} )
+    check_dependencies( missingDeps additionalDeps FILES ${sourceFiles} EXPECTED_DEPS ${ALL_DEPENDENCIES} ${moduleName} )
     if ( missingDeps )
         message ( WARNING "The module ${moduleName} depends on ${missingDeps} which are not listed as dependencies!" )
     endif()
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index bb3cb83c6279ccbd2e356747199045996afcebf6..3a36f7e0c5993f3ed44d1c21590a99a2ebc133aa 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -16,7 +16,7 @@ if __name__ == '__main__':
 
     mpd = Module(args.path)
     ps = mpd.add(data.ParticleStorage())
-    ps.set_shapes('Sphere', 'HalfSpace', 'CylindricalBoundary', 'Box', 'Ellipsoid')
+    ps.set_shapes('Sphere', 'HalfSpace', 'CylindricalBoundary', 'Box', 'Ellipsoid', 'ConvexPolyhedron')
     ps.add_property("position", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
     ps.add_property("linearVelocity", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
     ps.add_property("invMass", "walberla::real_t", defValue="real_t(1)", syncMode="ON_GHOST_CREATION")
diff --git a/src/mesa_pd/CMakeLists.txt b/src/mesa_pd/CMakeLists.txt
index d953b4269777c9357add6fcf2237f591437d26f6..6fdb08da0a7bebc3c83fc71944fa74af7cdfc613 100644
--- a/src/mesa_pd/CMakeLists.txt
+++ b/src/mesa_pd/CMakeLists.txt
@@ -5,4 +5,11 @@
 #
 ###################################################################################################
 
-waLBerla_add_module( DEPENDS blockforest core stencil vtk )
+waLBerla_add_module( DEPENDS blockforest core stencil vtk OPTIONAL_DEPENDS mesh_common )
+if(OPENMESH_CORE_FOUND)
+    set( WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE ON CACHE INTERNAL "")
+    message( STATUS "MESA-PD: ConvexPolyhedron shape is available (OpenMesh dependency satisfied)" )
+else()
+    set( WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE OFF CACHE INTERNAL "")
+    message( STATUS "MESA-PD: ConvexPolyhedron shape is unavailable (OpenMesh not found)" )
+endif()
\ No newline at end of file
diff --git a/src/mesa_pd/data/ShapeStorage.h b/src/mesa_pd/data/ShapeStorage.h
index a8baee5c713ddf991286a6e6ef6bc21055d4affb..7fdb72a6a10580f6edaa590936a31c827ec62745 100644
--- a/src/mesa_pd/data/ShapeStorage.h
+++ b/src/mesa_pd/data/ShapeStorage.h
@@ -34,6 +34,7 @@
 #include <mesa_pd/data/shape/CylindricalBoundary.h>
 #include <mesa_pd/data/shape/Box.h>
 #include <mesa_pd/data/shape/Ellipsoid.h>
+#include <mesa_pd/data/shape/ConvexPolyhedron.h>
 
 #include <core/Abort.h>
 #include <core/debug/Debug.h>
@@ -64,12 +65,17 @@ static_assert( Sphere::SHAPE_TYPE != HalfSpace::SHAPE_TYPE, "Shape types have to
 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!" );
 static_assert( Sphere::SHAPE_TYPE != Ellipsoid::SHAPE_TYPE, "Shape types have to be different!" );
+static_assert( Sphere::SHAPE_TYPE != ConvexPolyhedron::SHAPE_TYPE, "Shape types have to be different!" );
 static_assert( HalfSpace::SHAPE_TYPE != CylindricalBoundary::SHAPE_TYPE, "Shape types have to be different!" );
 static_assert( HalfSpace::SHAPE_TYPE != Box::SHAPE_TYPE, "Shape types have to be different!" );
 static_assert( HalfSpace::SHAPE_TYPE != Ellipsoid::SHAPE_TYPE, "Shape types have to be different!" );
+static_assert( HalfSpace::SHAPE_TYPE != ConvexPolyhedron::SHAPE_TYPE, "Shape types have to be different!" );
 static_assert( CylindricalBoundary::SHAPE_TYPE != Box::SHAPE_TYPE, "Shape types have to be different!" );
 static_assert( CylindricalBoundary::SHAPE_TYPE != Ellipsoid::SHAPE_TYPE, "Shape types have to be different!" );
+static_assert( CylindricalBoundary::SHAPE_TYPE != ConvexPolyhedron::SHAPE_TYPE, "Shape types have to be different!" );
 static_assert( Box::SHAPE_TYPE != Ellipsoid::SHAPE_TYPE, "Shape types have to be different!" );
+static_assert( Box::SHAPE_TYPE != ConvexPolyhedron::SHAPE_TYPE, "Shape types have to be different!" );
+static_assert( Ellipsoid::SHAPE_TYPE != ConvexPolyhedron::SHAPE_TYPE, "Shape types have to be different!" );
 
 template <typename ShapeT, typename... Args>
 size_t ShapeStorage::create(Args&&... args)
@@ -90,6 +96,7 @@ ReturnType ShapeStorage::singleDispatch( ParticleStorage& ps, size_t idx, func&
       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()));
       case Ellipsoid::SHAPE_TYPE : return f(ps, idx, *static_cast<Ellipsoid*>(shapes[ps.getShapeID(idx)].get()));
+      case ConvexPolyhedron::SHAPE_TYPE : return f(ps, idx, *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idx)].get()));
       default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idx)]->getShapeType() << ") could not be determined!");
    }
 }
@@ -130,6 +137,11 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
                                                    *static_cast<Sphere*>(shapes[ps.getShapeID(idx)].get()),
                                                    idy,
                                                    *static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
+            case ConvexPolyhedron::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<Sphere*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idy)].get()));
             default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
          }
       case HalfSpace::SHAPE_TYPE :
@@ -160,6 +172,11 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
                                                    *static_cast<HalfSpace*>(shapes[ps.getShapeID(idx)].get()),
                                                    idy,
                                                    *static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
+            case ConvexPolyhedron::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<HalfSpace*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idy)].get()));
             default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
          }
       case CylindricalBoundary::SHAPE_TYPE :
@@ -190,6 +207,11 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
                                                    *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idx)].get()),
                                                    idy,
                                                    *static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
+            case ConvexPolyhedron::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idy)].get()));
             default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
          }
       case Box::SHAPE_TYPE :
@@ -220,6 +242,11 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
                                                    *static_cast<Box*>(shapes[ps.getShapeID(idx)].get()),
                                                    idy,
                                                    *static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
+            case ConvexPolyhedron::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<Box*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idy)].get()));
             default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
          }
       case Ellipsoid::SHAPE_TYPE :
@@ -250,6 +277,46 @@ ReturnType ShapeStorage::doubleDispatch( ParticleStorage& ps, size_t idx, size_t
                                                    *static_cast<Ellipsoid*>(shapes[ps.getShapeID(idx)].get()),
                                                    idy,
                                                    *static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
+            case ConvexPolyhedron::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<Ellipsoid*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idy)].get()));
+            default : WALBERLA_ABORT("Shape type (" << shapes[ps.getShapeID(idy)]->getShapeType() << ") could not be determined!");
+         }
+      case ConvexPolyhedron::SHAPE_TYPE :
+         switch (shapes[ps.getShapeID(idy)]->getShapeType())
+         {
+            case Sphere::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<Sphere*>(shapes[ps.getShapeID(idy)].get()));
+            case HalfSpace::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<HalfSpace*>(shapes[ps.getShapeID(idy)].get()));
+            case CylindricalBoundary::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<CylindricalBoundary*>(shapes[ps.getShapeID(idy)].get()));
+            case Box::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<Box*>(shapes[ps.getShapeID(idy)].get()));
+            case Ellipsoid::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<Ellipsoid*>(shapes[ps.getShapeID(idy)].get()));
+            case ConvexPolyhedron::SHAPE_TYPE : return f(ps,
+                                                   idx,
+                                                   *static_cast<ConvexPolyhedron*>(shapes[ps.getShapeID(idx)].get()),
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(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/ConvexPolyhedron.h b/src/mesa_pd/data/shape/ConvexPolyhedron.h
new file mode 100644
index 0000000000000000000000000000000000000000..cad5c65863e3ad4f2fc7a1bace071454c787aecd
--- /dev/null
+++ b/src/mesa_pd/data/shape/ConvexPolyhedron.h
@@ -0,0 +1,356 @@
+//======================================================================================================================
+//
+//  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 ConvexPolyhedron.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/Abort.h"
+
+#include <waLBerlaDefinitions.h>
+
+#include <mesa_pd/data/shape/BaseShape.h>
+
+
+#ifdef WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE
+
+/**
+ * Full ConvexPolyhedron shape supporting all features, as the mesh_common module is available.
+ */
+
+#include "mesh_common/TriangleMeshes.h"
+#include "mesh_common/MeshOperations.h"
+
+#include <core/math/Constants.h>
+
+#include "core/mpi/BufferDataTypeExtensions.h"
+#include "mesh_common/OpenMeshBufferTypeExtensions.h"
+
+namespace walberla {
+namespace mesa_pd {
+namespace data {
+
+using namespace walberla::mesa_pd;
+
+/**
+ * \brief Convex Polyhedron shape containing a mesh.
+ *
+ * Use ConvexPolyhedron::updateMeshQuantities() to initialize the particle properties (volume, mass, inertia).
+ * This is done automatically upon initialization (and after shape unpacking), but if the mesh is manipulated
+ * while the particle already exists, it has to be called manually.
+ *
+ * \attention The origin of the meshes coordinate system is required to be at the position of its centroid! As such, the particle and mesh COS also overlap.
+ */
+class ConvexPolyhedron : public walberla::mesa_pd::data::BaseShape {
+public:
+   explicit ConvexPolyhedron(const mesh::TriangleMesh& mesh = mesh::TriangleMesh())
+         : BaseShape(ConvexPolyhedron::SHAPE_TYPE), mesh_(mesh)
+   {
+      if (mesh_.n_vertices() > 0) {
+         WALBERLA_CHECK_GREATER_EQUAL(mesh_.n_vertices(), 4);
+         updateMeshQuantities();
+      }
+   }
+
+   constexpr static int IS_AVAILABLE = true;
+
+   const mesh::TriangleMesh& getMesh() const { return mesh_; }
+   
+   mesh::TriangleMesh::VertexHandle supportVertex( const mesh::TriangleMesh::Normal & d,
+                                                   const mesh::TriangleMesh::VertexHandle startVertex ) const;
+
+   void updateMeshQuantities();
+
+   real_t getBoundingSphereRadius() const;
+   real_t getVolume() const override;
+   void updateMassAndInertia(const real_t density) override;
+
+   Vec3 support( const Vec3& d ) const override;
+
+   void pack(walberla::mpi::SendBuffer& buf) override;
+   void unpack(walberla::mpi::RecvBuffer& buf) override;
+
+   constexpr static int SHAPE_TYPE = 5; ///< Unique shape type identifier for convex polyhedrons.\ingroup mesa_pd_shape
+
+   bool meshQuantitiesAvailable = false;
+   mesh::TriangleMesh::VertexHandle octandVertices_[8];
+   
+private:
+   mesh::TriangleMesh mesh_;
+
+   real_t unitMass_; // mass for density 1, equals volume
+   Mat3 unitInertiaBF_; // inertia for density 1 relative to centroid
+   real_t boundingSphereRadius_;
+};
+
+/**
+ * \brief Get the bounding sphere radius around the centroid.
+ * \return Interaction radius.
+ */
+inline real_t ConvexPolyhedron::getBoundingSphereRadius() const {
+   return boundingSphereRadius_;
+}
+
+inline real_t ConvexPolyhedron::getVolume() const {
+   WALBERLA_CHECK(meshQuantitiesAvailable,
+         "unitMass and unitInertia first have to be initialized using updateMeshQuantities!");
+
+   return unitMass_; // V = m*d = m*1 = m
+}
+
+inline void ConvexPolyhedron::updateMassAndInertia(const real_t density) {
+   WALBERLA_CHECK(meshQuantitiesAvailable,
+                  "unitMass and unitInertia first have to be initialized using updateMeshQuantities!");
+
+   const real_t m = unitMass_ * density;
+   const Mat3 I = unitInertiaBF_ * density;
+
+   mass_ = m;
+   invMass_ = real_t(1.) / m;
+
+   inertiaBF_ = I;
+   invInertiaBF_ = I.getInverse();
+}
+
+inline void ConvexPolyhedron::updateMeshQuantities() {
+   WALBERLA_CHECK_GREATER(mesh_.n_vertices(), 0, "Cannot compute mesh quantities for an empty mesh!");
+
+   octandVertices_[0] = supportVertex(mesh::TriangleMesh::Normal(real_t(1), real_t(1), real_t(1)), *mesh_.vertices_begin());
+   octandVertices_[1] = supportVertex(mesh::TriangleMesh::Normal(real_t(1), real_t(1), real_t(-1)), *mesh_.vertices_begin());
+   octandVertices_[2] = supportVertex(mesh::TriangleMesh::Normal(real_t(1), real_t(-1), real_t(1)), *mesh_.vertices_begin());
+   octandVertices_[3] = supportVertex(mesh::TriangleMesh::Normal(real_t(1), real_t(-1), real_t(-1)), *mesh_.vertices_begin());
+   octandVertices_[4] = supportVertex(mesh::TriangleMesh::Normal(real_t(-1), real_t(1), real_t(1)), *mesh_.vertices_begin());
+   octandVertices_[5] = supportVertex(mesh::TriangleMesh::Normal(real_t(-1), real_t(1), real_t(-1)), *mesh_.vertices_begin());
+   octandVertices_[6] = supportVertex(mesh::TriangleMesh::Normal(real_t(-1), real_t(-1), real_t(1)), *mesh_.vertices_begin());
+   octandVertices_[7] = supportVertex(mesh::TriangleMesh::Normal(real_t(-1), real_t(-1), real_t(-1)), *mesh_.vertices_begin());
+
+   Vec3 centroid;
+   mesh::computeMassProperties(mesh_, real_t(1), centroid, unitInertiaBF_, unitMass_);
+   WALBERLA_CHECK_FLOAT_EQUAL(centroid, Vector3<real_t>(0,0,0), "The mesh has to have its centroid at the origin of its coordinate system! Use mesh::computeCentroid and mesh::translate.");
+
+   real_t maxSqRadius(0);
+   for(auto vh : mesh_.vertices()) {
+      auto v = mesh::toWalberla(mesh_.point(vh));
+      auto centroidToVSqr = v.sqrLength();
+
+      if (centroidToVSqr > maxSqRadius) {
+         maxSqRadius = centroidToVSqr;
+      }
+   }
+   boundingSphereRadius_ = std::sqrt(maxSqRadius);
+
+   meshQuantitiesAvailable = true;
+}
+
+inline Vec3 ConvexPolyhedron::support( const Vec3& d ) const {
+   WALBERLA_CHECK_UNEQUAL(mesh_.n_vertices(), 0, "Cannot compute support for a mesh with 0 vertices.");
+   WALBERLA_CHECK(meshQuantitiesAvailable, "Octand vertices of this mesh first have to be initialized using updateMeshQuantities!");
+   // taken from pe implementation
+
+   if (math::equal(d.length(), real_t(0))) return Vec3(0,0,0);
+
+   // d is already in the mesh coordinate system, as the origins of the particle and the mesh COS overlap
+   mesh::TriangleMesh::Normal d_loc = mesh::toOpenMesh(d);
+
+   mesh::TriangleMesh::VertexHandle startVertex;
+   if(d_loc[0] >= real_t( 0 ))
+   {
+      if(d_loc[1] >= real_t( 0 ))
+      {
+         startVertex = d_loc[2] >= real_t( 0 ) ? octandVertices_[0] : octandVertices_[1];
+      }
+      else // d_loc[1] < 0
+      {
+         startVertex = d_loc[2] >= real_t( 0 ) ? octandVertices_[2] : octandVertices_[3];
+      }
+   }
+   else // d_loc[0] < 0
+   {
+      if(d_loc[1] >= real_t( 0 ))
+      {
+         startVertex = d_loc[2] >= real_t( 0 ) ? octandVertices_[4] : octandVertices_[5];
+      }
+      else // d_loc[1] < 0
+      {
+         startVertex = d_loc[2] >= real_t( 0 ) ? octandVertices_[6] : octandVertices_[7];
+      }
+   }
+
+   mesh::TriangleMesh::VertexHandle vh = supportVertex( d_loc, startVertex );
+
+   // the resulting vertex has to be shifted back into the mesh particle coordinate system
+   auto relativeSupport = mesh::toWalberla(mesh_.point( vh ));
+
+   //WALBERLA_LOG_INFO("Conv poly support: " << relativeSupport << ", d: " << d);
+
+   return relativeSupport;
+}
+
+inline mesh::TriangleMesh::VertexHandle ConvexPolyhedron::supportVertex( const mesh::TriangleMesh::Normal & d,
+      const mesh::TriangleMesh::VertexHandle startVertex ) const {
+   // taken from pe implementation
+
+   mesh::TriangleMesh::VertexHandle maxScalarProductVertex = startVertex;
+   real_t maxScalarProduct = mesh_.point(maxScalarProductVertex) | d;
+
+   bool isExtremum = false;
+   while( !isExtremum ) {
+      isExtremum = true;
+      for(auto vh : mesh_.vv_range( maxScalarProductVertex ))
+      {
+         const real_t sp = mesh_.point(vh) | d;
+         if(sp > maxScalarProduct)
+         {
+            isExtremum = false;
+            maxScalarProductVertex = vh;
+            maxScalarProduct = sp;
+            break;
+         }
+      }
+   }
+
+   return maxScalarProductVertex;
+}
+
+inline void ConvexPolyhedron::pack(walberla::mpi::SendBuffer& buf){
+   // partially taken from pe implementation
+
+   BaseShape::pack(buf);
+
+   buf << mesh_.n_vertices(); // number of vertices
+
+   int dbgIndex = 0;
+   WALBERLA_UNUSED(dbgIndex);
+   for(auto vh : mesh_.vertices()) {
+      WALBERLA_ASSERT_EQUAL( vh.idx(), dbgIndex++ ); // assume vertices are compactly stored
+
+      buf << mesh_.point(vh);
+   }
+
+   buf << mesh_.n_faces();
+
+   for(auto fh: mesh_.faces()) {
+      for (auto vhIt = mesh_.cfv_ccwbegin(fh); vhIt != mesh_.cfv_ccwend(fh); ++vhIt) {
+         WALBERLA_ASSERT_GREATER_EQUAL(vhIt->idx(), 0);
+         WALBERLA_ASSERT_LESS(vhIt->idx(), mesh_.n_vertices());
+
+         buf << vhIt->idx();
+      }
+   }
+}
+
+inline void ConvexPolyhedron::unpack(walberla::mpi::RecvBuffer& buf){
+   // partially taken from pe implementation
+
+   BaseShape::unpack(buf);
+
+   WALBERLA_CHECK_EQUAL(mesh_.n_vertices(), 0, "Mesh needs to be empty!");
+   WALBERLA_CHECK_EQUAL(mesh_.n_faces(), 0, "Mesh needs to be empty!");
+
+   size_t numVertices;
+   buf >> numVertices;
+   std::vector<mesh::TriangleMesh::VertexHandle> vertexHandles(numVertices);
+   for(size_t i = 0; i < numVertices; ++i) {
+      mesh::TriangleMesh::Point p;
+      buf >> p;
+      vertexHandles[size_t(i)] = mesh_.add_vertex( p );
+   }
+
+   size_t numFaces;
+   buf >> numFaces;
+   for(size_t i = 0; i < numFaces; ++i){
+      int v0;
+      int v1;
+      int v2;
+      buf >> v0 >> v1 >> v2;
+      WALBERLA_ASSERT_GREATER_EQUAL( v0, 0 );
+      WALBERLA_ASSERT_GREATER_EQUAL( v1, 0 );
+      WALBERLA_ASSERT_GREATER_EQUAL( v2, 0 );
+      WALBERLA_ASSERT_LESS( v0, numVertices );
+      WALBERLA_ASSERT_LESS( v1, numVertices );
+      WALBERLA_ASSERT_LESS( v2, numVertices );
+
+      mesh_.add_face( vertexHandles[size_t(v0)], vertexHandles[size_t(v1)], vertexHandles[size_t(v2)] );
+   }
+
+   updateMeshQuantities();
+}
+
+}
+}
+}
+
+#else
+
+/**
+ * Replacement for the "full" ConvexPolyhedron, which will throw compile time errors if any of its features are
+ * actually used.
+ */
+
+namespace walberla {
+
+namespace mesh {
+// forward declaration failing if ConvexPolyhedron is actually used
+class TriangleMesh;
+}
+
+namespace mesa_pd {
+namespace data {
+
+using namespace walberla::mesa_pd;
+
+class ConvexPolyhedron : public walberla::mesa_pd::data::BaseShape {
+public:
+   explicit ConvexPolyhedron(const mesh::TriangleMesh&) : BaseShape(ConvexPolyhedron::SHAPE_TYPE) {
+      WALBERLA_ABORT("Shape ConvexPolyhedron is not available! Ensure waLBerla is configured with OpenMesh support.");
+   }
+
+   ConvexPolyhedron() : BaseShape(ConvexPolyhedron::SHAPE_TYPE) {
+      WALBERLA_ABORT("Shape ConvexPolyhedron is not available! Ensure waLBerla is configured with OpenMesh support.");
+   }
+
+   constexpr static int IS_AVAILABLE = false;
+   constexpr static int SHAPE_TYPE = 5; ///< Unique shape type identifier for convex polyhedrons.\ingroup mesa_pd_shape
+
+   void updateMassAndInertia(const real_t /*density*/) override {
+      WALBERLA_ABORT("Shape ConvexPolyhedron is not available! Ensure waLBerla is configured with OpenMesh support.");
+   };
+
+   real_t getVolume() const override {
+      WALBERLA_ABORT("Shape ConvexPolyhedron is not available! Ensure waLBerla is configured with OpenMesh support.");
+   }
+
+   Vec3 support( const Vec3& /*d*/ ) const override {
+      WALBERLA_ABORT("Shape ConvexPolyhedron is not available! Ensure waLBerla is configured with OpenMesh support.");
+   }
+
+   void pack(walberla::mpi::SendBuffer& /*buf*/) override {
+      WALBERLA_ABORT("Shape ConvexPolyhedron is not available! Ensure waLBerla is configured with OpenMesh support.");
+   };
+
+   void unpack(walberla::mpi::RecvBuffer& /*buf*/) override {
+      WALBERLA_ABORT("Shape ConvexPolyhedron is not available! Ensure waLBerla is configured with OpenMesh support.");
+   };
+};
+
+}
+}
+}
+
+
+#endif
\ No newline at end of file
diff --git a/src/mesa_pd/kernel/DoubleCast.h b/src/mesa_pd/kernel/DoubleCast.h
index 46e917cc7adf7f2360763256b6c53c51b4a762ac..6f7851bfe4f86670954880daf5d87090de801f82 100644
--- a/src/mesa_pd/kernel/DoubleCast.h
+++ b/src/mesa_pd/kernel/DoubleCast.h
@@ -34,6 +34,7 @@
 #include <mesa_pd/data/shape/CylindricalBoundary.h>
 #include <mesa_pd/data/shape/Box.h>
 #include <mesa_pd/data/shape/Ellipsoid.h>
+#include <mesa_pd/data/shape/ConvexPolyhedron.h>
 
 #include <core/Abort.h>
 #include <core/debug/Debug.h>
@@ -94,6 +95,11 @@ auto DoubleCast::operator()( size_t idx, size_t idy, Accessor& ac, func& f, Args
                                                    *static_cast<Sphere*>(ac.getShape(idx)),
                                                    *static_cast<Ellipsoid*>(ac.getShape(idy)),
                                                    std::forward<Args>(args)...);
+            case ConvexPolyhedron::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<Sphere*>(ac.getShape(idx)),
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idy)),
+                                                   std::forward<Args>(args)...);
             default : WALBERLA_ABORT("Shape type (" << ac.getShape(idy)->getShapeType() << ") could not be determined!");
          }
       case HalfSpace::SHAPE_TYPE :
@@ -124,6 +130,11 @@ auto DoubleCast::operator()( size_t idx, size_t idy, Accessor& ac, func& f, Args
                                                    *static_cast<HalfSpace*>(ac.getShape(idx)),
                                                    *static_cast<Ellipsoid*>(ac.getShape(idy)),
                                                    std::forward<Args>(args)...);
+            case ConvexPolyhedron::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<HalfSpace*>(ac.getShape(idx)),
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idy)),
+                                                   std::forward<Args>(args)...);
             default : WALBERLA_ABORT("Shape type (" << ac.getShape(idy)->getShapeType() << ") could not be determined!");
          }
       case CylindricalBoundary::SHAPE_TYPE :
@@ -154,6 +165,11 @@ auto DoubleCast::operator()( size_t idx, size_t idy, Accessor& ac, func& f, Args
                                                    *static_cast<CylindricalBoundary*>(ac.getShape(idx)),
                                                    *static_cast<Ellipsoid*>(ac.getShape(idy)),
                                                    std::forward<Args>(args)...);
+            case ConvexPolyhedron::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<CylindricalBoundary*>(ac.getShape(idx)),
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idy)),
+                                                   std::forward<Args>(args)...);
             default : WALBERLA_ABORT("Shape type (" << ac.getShape(idy)->getShapeType() << ") could not be determined!");
          }
       case Box::SHAPE_TYPE :
@@ -184,6 +200,11 @@ auto DoubleCast::operator()( size_t idx, size_t idy, Accessor& ac, func& f, Args
                                                    *static_cast<Box*>(ac.getShape(idx)),
                                                    *static_cast<Ellipsoid*>(ac.getShape(idy)),
                                                    std::forward<Args>(args)...);
+            case ConvexPolyhedron::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<Box*>(ac.getShape(idx)),
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idy)),
+                                                   std::forward<Args>(args)...);
             default : WALBERLA_ABORT("Shape type (" << ac.getShape(idy)->getShapeType() << ") could not be determined!");
          }
       case Ellipsoid::SHAPE_TYPE :
@@ -214,6 +235,46 @@ auto DoubleCast::operator()( size_t idx, size_t idy, Accessor& ac, func& f, Args
                                                    *static_cast<Ellipsoid*>(ac.getShape(idx)),
                                                    *static_cast<Ellipsoid*>(ac.getShape(idy)),
                                                    std::forward<Args>(args)...);
+            case ConvexPolyhedron::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<Ellipsoid*>(ac.getShape(idx)),
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idy)),
+                                                   std::forward<Args>(args)...);
+            default : WALBERLA_ABORT("Shape type (" << ac.getShape(idy)->getShapeType() << ") could not be determined!");
+         }
+      case ConvexPolyhedron::SHAPE_TYPE :
+         switch (ac.getShape(idy)->getShapeType())
+         {
+            case Sphere::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idx)),
+                                                   *static_cast<Sphere*>(ac.getShape(idy)),
+                                                   std::forward<Args>(args)...);
+            case HalfSpace::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idx)),
+                                                   *static_cast<HalfSpace*>(ac.getShape(idy)),
+                                                   std::forward<Args>(args)...);
+            case CylindricalBoundary::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idx)),
+                                                   *static_cast<CylindricalBoundary*>(ac.getShape(idy)),
+                                                   std::forward<Args>(args)...);
+            case Box::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idx)),
+                                                   *static_cast<Box*>(ac.getShape(idy)),
+                                                   std::forward<Args>(args)...);
+            case Ellipsoid::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idx)),
+                                                   *static_cast<Ellipsoid*>(ac.getShape(idy)),
+                                                   std::forward<Args>(args)...);
+            case ConvexPolyhedron::SHAPE_TYPE : return f(idx,
+                                                   idy,
+                                                   *static_cast<ConvexPolyhedron*>(ac.getShape(idx)),
+                                                   *static_cast<ConvexPolyhedron*>(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 12f85b93009ddba097dd3b8f5cf05b0a01ee0354..b8ae7571732ba1247c7ccbca3acf3048103f5a7c 100644
--- a/src/mesa_pd/kernel/SingleCast.h
+++ b/src/mesa_pd/kernel/SingleCast.h
@@ -34,6 +34,7 @@
 #include <mesa_pd/data/shape/CylindricalBoundary.h>
 #include <mesa_pd/data/shape/Box.h>
 #include <mesa_pd/data/shape/Ellipsoid.h>
+#include <mesa_pd/data/shape/ConvexPolyhedron.h>
 
 #include <core/Abort.h>
 #include <core/debug/Debug.h>
@@ -70,6 +71,7 @@ auto SingleCast::operator()( size_t idx, Accessor& ac, func& f, 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)...);
       case Ellipsoid::SHAPE_TYPE : return f(idx, *static_cast<Ellipsoid*>(ac.getShape(idx)), std::forward<Args>(args)...);
+      case ConvexPolyhedron::SHAPE_TYPE : return f(idx, *static_cast<ConvexPolyhedron*>(ac.getShape(idx)), std::forward<Args>(args)...);
       default : WALBERLA_ABORT("Shape type (" << ac.getShape(idx)->getShapeType() << ") could not be determined!");
    }
 }
diff --git a/src/mesa_pd/mpi/ShapePackUnpack.h b/src/mesa_pd/mpi/ShapePackUnpack.h
index 0b5ab10e22738c55482d50d36e08a62367014aa3..dd90703559ec4da6660a8fe5e1706bf01242d5ed 100644
--- a/src/mesa_pd/mpi/ShapePackUnpack.h
+++ b/src/mesa_pd/mpi/ShapePackUnpack.h
@@ -32,6 +32,7 @@
 #include <mesa_pd/data/shape/CylindricalBoundary.h>
 #include <mesa_pd/data/shape/Box.h>
 #include <mesa_pd/data/shape/Ellipsoid.h>
+#include <mesa_pd/data/shape/ConvexPolyhedron.h>
 
 #include <core/mpi/RecvBuffer.h>
 #include <core/mpi/SendBuffer.h>
@@ -83,6 +84,10 @@ namespace mpi {
             bs = std::make_unique<mesa_pd::data::Ellipsoid>();
             bs->unpack(buf);
             break;
+         case ConvexPolyhedron::SHAPE_TYPE :
+            bs = std::make_unique<mesa_pd::data::ConvexPolyhedron>();
+            bs->unpack(buf);
+            break;
          default : WALBERLA_ABORT("Shape type (" << shapeType << ") could not be determined!");
       }
       return buf;
diff --git a/src/mesa_pd/vtk/ConvexPolyhedron/MeshParticleVTKOutput.h b/src/mesa_pd/vtk/ConvexPolyhedron/MeshParticleVTKOutput.h
new file mode 100644
index 0000000000000000000000000000000000000000..f39b7c5bc32a7125120b472de3832603e498409d
--- /dev/null
+++ b/src/mesa_pd/vtk/ConvexPolyhedron/MeshParticleVTKOutput.h
@@ -0,0 +1,201 @@
+//======================================================================================================================
+//
+//  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 MeshVTKOutput.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <OpenMesh/Core/Utils/PropertyManager.hh>
+#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/data/ShapeStorage.h>
+#include <mesa_pd/data/shape/ConvexPolyhedron.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/Types.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/DataSourceAdapters.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/FaceDataSource.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/OutputSelectorFaceDataSource.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/OutputSelectorVertexDataSource.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/VertexDataSource.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/tesselation/ConvexPolyhedronTesselation.h>
+#include <mesh_common/vtk/DistributedVTKMeshWriter.h>
+#include <utility>
+
+namespace walberla {
+namespace mesa_pd {
+
+template<typename MeshType>
+class MeshParticleVTKOutput {
+   static_assert(MeshType::IsPolyMesh == 1, "We need polygonal meshes here!");
+
+public:
+   using ParticleSelectorFunc = std::function<bool (const walberla::mesa_pd::data::ParticleStorage::iterator& pIt)>;
+
+   typedef typename mesh::DistributedVTKMeshWriter<MeshType>::Vertices Vertices;
+
+   MeshParticleVTKOutput( shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps,
+                          shared_ptr<walberla::mesa_pd::data::ShapeStorage> ss, const std::string & identifier,
+                          const uint_t writeFrequency, const std::string & baseFolder = "vtk_out")
+                          : ps_(std::move(ps)), ss_(std::move(ss)), mesh_(make_shared<MeshType>()),
+                          faceToParticleIdxManager_(*mesh_, "particle"), vertexToParticleIdxManager_(*mesh_, "particle"),
+                          meshWriter_(mesh_, identifier, writeFrequency, baseFolder) {
+
+   }
+
+   void setParticleSelector( const ParticleSelectorFunc& func) {particleSelector_ = func;}
+   void assembleMesh();
+
+   template <typename Selector>
+   void addVertexOutput(const std::string& name);
+   template<typename Selector>
+   void addFaceOutput(const std::string &name);
+
+   template <typename DataSourceType>
+   void addVertexDataSource(const shared_ptr<DataSourceType> & dataSource);
+   template <typename DataSourceType>
+   void addFaceDataSource(const shared_ptr<DataSourceType> & dataSource);
+
+   void operator()() {
+      assembleMesh();
+      // the mesh writer writes the mesh to vtk files and adds properties as defined by the data sources
+      meshWriter_();
+      mesh_->clean(); // the output mesh is no longer needed, thus discard its contents
+   }
+
+private:
+   const shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps_;
+   const shared_ptr<walberla::mesa_pd::data::ShapeStorage> ss_;
+   shared_ptr<MeshType> mesh_; ///< the output mesh
+
+   // these "managers" (which are just maps basically) map the faces and vertices of the output mesh to particles
+   // such that we can assign those faces and vertices properties of the particle (i.e. ID or velocity)
+   ParticleIdxFacePropertyManager<MeshType> faceToParticleIdxManager_;
+   ParticleIdxVertexPropertyManager<MeshType> vertexToParticleIdxManager_;
+
+   mesh::DistributedVTKMeshWriter<MeshType> meshWriter_;
+
+   ParticleSelectorFunc particleSelector_ = [](const walberla::mesa_pd::data::ParticleStorage::iterator& /*pIt*/){
+      return true;
+   };
+};
+
+/**
+ * \brief Adds a output selector for the vertices of the mesh
+ * Similar to MESAPD's vtk output, one can add a property selector to select a property for a vertex
+ * corresponding to the associated particle.
+ *
+ * \tparam MeshType
+ * \tparam Selector Type of the selector
+ * \param name Name of the output
+ */
+template<typename MeshType>
+template<typename Selector>
+void MeshParticleVTKOutput<MeshType>::addVertexOutput(const std::string &name) {
+   typedef OutputSelectorVertexDataSource<MeshType, Selector, typename Selector::return_type> DataSourceType;
+   auto ds = make_shared<DataSourceType>(name, Selector());
+   addVertexDataSource(ds);
+}
+
+/**
+ * \brief Adds a output selector for the faces of the mesh
+ * Similar to MESAPD's vtk output, one can add a property selector to select a property for a face
+ * corresponding to the associated particle.
+ *
+ * \tparam MeshType
+ * \tparam Selector Type of the selector
+ * \param name Name of the output
+ */
+template<typename MeshType>
+template<typename Selector>
+void MeshParticleVTKOutput<MeshType>::addFaceOutput(const std::string &name) {
+   typedef OutputSelectorFaceDataSource<MeshType, Selector, typename Selector::return_type> DataSourceType;
+   auto ds = make_shared<DataSourceType>(name, Selector());
+   addFaceDataSource(ds);
+}
+
+/**
+ * \brief Add a vertex data source assigning a piece of data to a vertex.
+ * \tparam MeshType
+ * \tparam DataSourceType Type of data source (has to be derived from VertexDataSource).
+ * \param dataSource Data source responsible for picking data for a vertex.
+ */
+template<typename MeshType>
+template<typename DataSourceType>
+void MeshParticleVTKOutput<MeshType>::addVertexDataSource(const shared_ptr<DataSourceType> & dataSource) {
+   typedef internal::VertexDataSourceAdapter<MeshType, typename DataSourceType::ComponentType> AdapterType;
+   meshWriter_.addDataSource(make_shared<AdapterType>(dataSource, vertexToParticleIdxManager_, ps_));
+}
+
+/**
+ * \brief Add a face data source assigning a piece of data to a face.
+ * \tparam MeshType
+ * \tparam DataSourceType Type of data source (has to be derived from FaceDataSource).
+ * \param dataSource Data source responsible for picking data for a face.
+ */
+template<typename MeshType>
+template<typename DataSourceType>
+void MeshParticleVTKOutput<MeshType>::addFaceDataSource(const shared_ptr<DataSourceType> & dataSource) {
+   typedef internal::FaceDataSourceAdapter<MeshType, typename DataSourceType::ComponentType> AdapterType;
+   meshWriter_.addDataSource(make_shared<AdapterType>(dataSource, faceToParticleIdxManager_, ps_));
+}
+
+
+/**
+ * \brief Creates the output mesh and writes it to mesh_.
+ * \tparam MeshType
+ */
+template<typename MeshType>
+void MeshParticleVTKOutput<MeshType>::assembleMesh() {
+   // those will save the newly created vertices and faces for each mesh
+   // to make it possible to map a vertex/face to the corresponding particle
+   std::vector<typename MeshType::VertexHandle> newVertices;
+   std::vector<typename MeshType::FaceHandle> newFaces;
+
+   // ensure the mesh is empty, as this will contain the new output
+   mesh_->clean();
+
+   // then iterate over every particle and tesselate it to include it in the output mesh
+   for (auto pIt = ps_->begin(); pIt != ps_->end(); ++pIt) {
+      if (!particleSelector_(pIt)) continue;
+
+      auto& shape = ss_->shapes[pIt->getShapeID()];
+
+      newVertices.clear();
+      newFaces.clear();
+
+      if (shape->getShapeType() == walberla::mesa_pd::data::ConvexPolyhedron::SHAPE_TYPE) {
+         const auto& convexPolyhedron = *static_cast<walberla::mesa_pd::data::ConvexPolyhedron*>(shape.get());
+         const auto& particle = *pIt;
+
+         // tesselate: add the shape at the particle's position into the output mesh
+         tesselate(convexPolyhedron, particle, mesh_, newVertices, newFaces);
+      }
+
+      // save particle idx to managers
+      for (const auto & vertex: newVertices) {
+         vertexToParticleIdxManager_[vertex] = pIt->getIdx();
+      }
+      for (const auto & face: newFaces) {
+         faceToParticleIdxManager_[face] = pIt->getIdx();
+      }
+   }
+
+   //WALBERLA_LOG_INFO("MESA-PD VTK output mesh contains " << mesh_->n_vertices() << " vertices.")
+}
+
+
+}
+}
\ No newline at end of file
diff --git a/src/mesa_pd/vtk/ConvexPolyhedron/Types.h b/src/mesa_pd/vtk/ConvexPolyhedron/Types.h
new file mode 100644
index 0000000000000000000000000000000000000000..46aa1bbcd8ce1a5a9309406621b70c49ba2b0d3f
--- /dev/null
+++ b/src/mesa_pd/vtk/ConvexPolyhedron/Types.h
@@ -0,0 +1,16 @@
+#include <OpenMesh/Core/Utils/Property.hh>
+#include <OpenMesh/Core/Utils/PropertyManager.hh>
+
+namespace walberla {
+namespace mesa_pd {
+
+typedef typename OpenMesh::FPropHandleT<size_t> ParticleIdxFacePropertyHandle;
+typedef typename OpenMesh::VPropHandleT<size_t> ParticleIdxVertexPropertyHandle;
+template <typename MeshType>
+using ParticleIdxFacePropertyManager = OpenMesh::PropertyManager<ParticleIdxFacePropertyHandle, MeshType>;
+template <typename MeshType>
+using ParticleIdxVertexPropertyManager = typename OpenMesh::PropertyManager<ParticleIdxVertexPropertyHandle, MeshType>;
+
+
+}
+}
\ No newline at end of file
diff --git a/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/DataSourceAdapters.h b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/DataSourceAdapters.h
new file mode 100644
index 0000000000000000000000000000000000000000..8eb8a03f4586c3455e2b2fd3e31216ba7e1ca782
--- /dev/null
+++ b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/DataSourceAdapters.h
@@ -0,0 +1,92 @@
+//======================================================================================================================
+//
+//  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 DataSourceAdapters.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <core/DataTypes.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/FaceDataSource.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/VertexDataSource.h>
+#include <mesh_common/vtk/DistributedVTKMeshWriter.h>
+
+namespace walberla {
+namespace mesa_pd {
+namespace internal {
+
+/**
+ * \brief Adapts a vertex data source for the MESAPD mesh output to the generic vertex data source class.
+ * \tparam MeshType
+ * \tparam T output type
+ */
+template<typename MeshType, typename T >
+class VertexDataSourceAdapter : public mesh::DistributedVTKMeshWriter<MeshType>::template VertexDataSource<T> {
+public:
+   typedef typename mesh::DistributedVTKMeshWriter<MeshType>::template VertexDataSource<T>::Vertices Vertices;
+
+   VertexDataSourceAdapter( const shared_ptr<VertexDataSource<MeshType, T>> & vertexDataSource,
+         const ParticleIdxVertexPropertyManager<MeshType> & vertexToParticleIdxManager,
+         const shared_ptr<walberla::mesa_pd::data::ParticleStorage> & ps)
+         : mesh::DistributedVTKMeshWriter<MeshType>::template VertexDataSource<T>( vertexDataSource->name() ),
+               vertexDataSource_(vertexDataSource), vertexToParticleIdxManager_(vertexToParticleIdxManager),
+               ps_(ps) { }
+
+   virtual void getData(const MeshType & mesh, const Vertices & vertices, std::vector<T> & data) {
+      return vertexDataSource_->getData( mesh, vertices, data, vertexToParticleIdxManager_, ps_ );
+   };
+
+   virtual uint_t numComponents() { return vertexDataSource_->numComponents(); }
+
+protected:
+   shared_ptr<VertexDataSource<MeshType, T>> vertexDataSource_;
+   const ParticleIdxVertexPropertyManager<MeshType> & vertexToParticleIdxManager_;
+   const shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps_;
+};
+
+/**
+ * \brief Adapts a face data source for the MESAPD mesh output to the generic face data source class.
+ * \tparam MeshType
+ * \tparam T output type
+ */
+template<typename MeshType, typename T >
+class FaceDataSourceAdapter : public mesh::DistributedVTKMeshWriter<MeshType>::template FaceDataSource<T> {
+public:
+   typedef typename mesh::DistributedVTKMeshWriter<MeshType>::template FaceDataSource<T>::Faces Faces;
+
+   FaceDataSourceAdapter(const shared_ptr<FaceDataSource<MeshType, T>> & faceDataSource,
+         const ParticleIdxFacePropertyManager<MeshType> & faceToParticleIdxManager,
+         shared_ptr<walberla::mesa_pd::data::ParticleStorage>  ps)
+   : mesh::DistributedVTKMeshWriter<MeshType>::template FaceDataSource<T>(faceDataSource->name()),
+         faceDataSource_(faceDataSource), faceToParticleIdxManager_(faceToParticleIdxManager),
+         ps_(std::move(ps)) { }
+
+   virtual void getData(const MeshType & mesh, const Faces & faces, std::vector<T> & data) {
+      return faceDataSource_->getData( mesh, faces, data, faceToParticleIdxManager_, ps_ );
+   };
+
+   virtual uint_t numComponents() { return faceDataSource_->numComponents(); }
+
+protected:
+   shared_ptr<FaceDataSource<MeshType, T>> faceDataSource_;
+   const ParticleIdxFacePropertyManager<MeshType> & faceToParticleIdxManager_;
+   const shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps_;
+};
+
+}
+}
+}
\ No newline at end of file
diff --git a/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/FaceDataSource.h b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/FaceDataSource.h
new file mode 100644
index 0000000000000000000000000000000000000000..94479be44d587c2c8465112058c167a15d13865b
--- /dev/null
+++ b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/FaceDataSource.h
@@ -0,0 +1,50 @@
+//======================================================================================================================
+//
+//  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 ParticleFaceDataSource.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/Types.h>
+#include <mesh_common/vtk/DistributedVTKMeshWriter.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+template<typename MeshType, typename T>
+class FaceDataSource {
+public:
+   typedef typename mesh::DistributedVTKMeshWriter<MeshType>::Faces Faces;
+
+   explicit FaceDataSource(const std::string & name) : name_(name) {}
+
+   const std::string & name() { return name_; }
+   virtual uint_t numComponents() = 0;
+   virtual void getData(const MeshType &, const Faces &faces, std::vector<T> &data,
+                        const ParticleIdxFacePropertyManager<MeshType> & faceToParticleIdxManager,
+                        shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps) = 0;
+
+   virtual ~FaceDataSource() {}
+
+protected:
+   std::string name_;
+};
+
+}
+}
\ No newline at end of file
diff --git a/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/OutputSelectorFaceDataSource.h b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/OutputSelectorFaceDataSource.h
new file mode 100644
index 0000000000000000000000000000000000000000..f1442527d9ee3fcb08515aedf7d402ff035b47ed
--- /dev/null
+++ b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/OutputSelectorFaceDataSource.h
@@ -0,0 +1,105 @@
+//======================================================================================================================
+//
+//  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 ParticleOutputSelectorFaceDataSource.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <core/DataTypes.h>
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/FaceDataSource.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+/**
+ * \brief Generic data source to pick data for a face using a particle selector.
+ * \attention The underlying mesh data sources don't support Vec3's etc. natively, thus specializations need to be given.
+ * \tparam MeshType
+ * \tparam Selector Type of the selector.
+ * \tparam Type Type of the data.
+ */
+template<typename MeshType, typename Selector, typename Type = typename Selector::return_type>
+class OutputSelectorFaceDataSource : public FaceDataSource<MeshType, Type> {
+public:
+   typedef FaceDataSource<MeshType, Type> Base;
+   typedef typename Base::Faces Faces;
+
+   typedef Type ComponentType;
+
+   OutputSelectorFaceDataSource(const std::string& name, Selector selector) : Base(name), selector_(selector) { }
+
+   virtual uint_t numComponents() {
+      return uint_t(1);
+   }
+
+   using Base::getData;
+   virtual void getData( const MeshType &, const Faces & faces, std::vector<Type> & data,
+                         const ParticleIdxFacePropertyManager<MeshType> & faceToParticleIdxManager,
+                         shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps) {
+      for (const auto & face: faces) {
+         size_t particleIdx = faceToParticleIdxManager[face];
+         auto p = (*ps)[particleIdx];
+         data.push_back(selector_(p));
+      }
+   };
+
+private:
+   Selector selector_;
+};
+
+/**
+ * \brief Data Source for particle selectors specialized for Vec3
+ * \tparam MeshType
+ * \tparam Selector
+ * \tparam Type
+ */
+template<typename MeshType, typename Selector, typename Type>
+class OutputSelectorFaceDataSource<MeshType, Selector, Vector3<Type>> : public FaceDataSource<MeshType, Type> {
+public:
+   typedef FaceDataSource<MeshType, Type> Base;
+   typedef typename Base::Faces Faces;
+
+   typedef Type ComponentType;
+
+   OutputSelectorFaceDataSource(const std::string& name, Selector selector) : Base(name), selector_(selector) { }
+
+   virtual uint_t numComponents() {
+      return uint_t(3);
+   }
+
+   using Base::getData;
+   virtual void getData( const MeshType &, const Faces & faces, std::vector<Type> & data,
+                         const ParticleIdxFacePropertyManager<MeshType> & faceToParticleIdxManager,
+                         shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps) {
+      for (const auto & face: faces) {
+         size_t particleIdx = faceToParticleIdxManager[face];
+         auto p = (*ps)[particleIdx];
+         const Vector3<Type> d = selector_(p);
+         data.push_back(d[0]);
+         data.push_back(d[1]);
+         data.push_back(d[2]);
+      }
+   };
+
+private:
+   Selector selector_;
+};
+
+}
+}
\ No newline at end of file
diff --git a/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/OutputSelectorVertexDataSource.h b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/OutputSelectorVertexDataSource.h
new file mode 100644
index 0000000000000000000000000000000000000000..1f8a7fd05e41ed5b78a399571264a53cc3a16697
--- /dev/null
+++ b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/OutputSelectorVertexDataSource.h
@@ -0,0 +1,97 @@
+//======================================================================================================================
+//
+//  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 ParticleOutputSelectorVertexDataSource.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <core/DataTypes.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/VertexDataSource.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+template<typename MeshType, typename Selector, typename Type = typename Selector::return_type>
+class OutputSelectorVertexDataSource : public VertexDataSource<MeshType, Type> {
+public:
+   typedef VertexDataSource<MeshType, Type> Base;
+   typedef typename Base::Vertices Vertices;
+
+   typedef Type ComponentType;
+
+   OutputSelectorVertexDataSource(const std::string& name, Selector selector) : Base(name), selector_(selector) { }
+
+   virtual uint_t numComponents() {
+      return uint_t(1);
+   }
+
+   using Base::getData;
+   virtual void getData( const MeshType &, const Vertices & vertices, std::vector<Type> & data,
+                         const ParticleIdxVertexPropertyManager<MeshType> & vertexToParticleIdxManager,
+                         shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps) {
+      for (const auto & vertex: vertices) {
+         size_t particleIdx = vertexToParticleIdxManager[vertex];
+         auto p = (*ps)[particleIdx];
+         data.push_back(selector_(p));
+      }
+   };
+
+private:
+   Selector selector_;
+};
+
+/**
+ * \brief Data Source specialized for Vec3
+ * \tparam MeshType
+ * \tparam Selector
+ * \tparam Type
+ */
+template<typename MeshType, typename Selector, typename Type>
+class OutputSelectorVertexDataSource<MeshType, Selector, Vector3<Type>> : public VertexDataSource<MeshType, Type> {
+public:
+   typedef VertexDataSource<MeshType, Type> Base;
+   typedef typename Base::Vertices Vertices;
+
+   typedef Type ComponentType;
+
+   OutputSelectorVertexDataSource(const std::string& name, Selector selector) : Base(name), selector_(selector) { }
+
+   virtual uint_t numComponents() {
+      return uint_t(3);
+   }
+
+   using Base::getData;
+   virtual void getData( const MeshType &, const Vertices & vertices, std::vector<Type> & data,
+                         const ParticleIdxVertexPropertyManager<MeshType> & vertexToParticleIdxManager,
+                         shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps) {
+      for (const auto & vertex: vertices) {
+         size_t particleIdx = vertexToParticleIdxManager[vertex];
+         auto p = (*ps)[particleIdx];
+         const Vector3<Type> d = selector_(p);
+         data.push_back(d[0]);
+         data.push_back(d[1]);
+         data.push_back(d[2]);
+      }
+   };
+
+private:
+   Selector selector_;
+};
+
+}
+}
\ No newline at end of file
diff --git a/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/SurfaceVelocityVertexDataSource.h b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/SurfaceVelocityVertexDataSource.h
new file mode 100644
index 0000000000000000000000000000000000000000..03eec9d9ed4d5d7bbf17b24c86b55e44389afb7f
--- /dev/null
+++ b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/SurfaceVelocityVertexDataSource.h
@@ -0,0 +1,63 @@
+//======================================================================================================================
+//
+//  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 ParticleOutputSelectorVertexDataSource.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <core/DataTypes.h>
+#include <mesa_pd/common/ParticleFunctions.h>
+#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/VertexDataSource.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+template<typename MeshType, typename Accessor, typename Type = real_t>
+class SurfaceVelocityVertexDataSource : public VertexDataSource<MeshType, Type> {
+public:
+   typedef VertexDataSource<MeshType, Type> Base;
+   typedef typename Base::Vertices Vertices;
+
+   typedef Type ComponentType;
+
+   SurfaceVelocityVertexDataSource(const std::string& name, const Accessor & ac) : Base(name), ac_(ac) { }
+
+   virtual uint_t numComponents() {
+      return uint_t(3);
+   }
+
+   using Base::getData;
+   virtual void getData( const MeshType & mesh, const Vertices & vertices, std::vector<Type> & data,
+                         const ParticleIdxVertexPropertyManager<MeshType> & vertexToParticleIdxManager,
+                         shared_ptr<walberla::mesa_pd::data::ParticleStorage>) {
+      for (const auto & vertex: vertices) {
+         size_t particleIdx = vertexToParticleIdxManager[vertex];
+         auto vertexPosition = mesh::toWalberlaNumericCast<real_t>(mesh.point(vertex));
+         const Vector3<Type> d = walberla::mesa_pd::getVelocityAtWFPoint(particleIdx, ac_, vertexPosition);
+         data.push_back(d[0]);
+         data.push_back(d[1]);
+         data.push_back(d[2]);
+      }
+   };
+private:
+   const Accessor & ac_;
+};
+
+}
+}
\ No newline at end of file
diff --git a/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/VertexDataSource.h b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/VertexDataSource.h
new file mode 100644
index 0000000000000000000000000000000000000000..bef238fa0de699c4da73597d6ddffe460bcf7add
--- /dev/null
+++ b/src/mesa_pd/vtk/ConvexPolyhedron/data_sources/VertexDataSource.h
@@ -0,0 +1,49 @@
+//======================================================================================================================
+//
+//  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 ParticleVertexDataSource.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/vtk/ConvexPolyhedron/Types.h>
+#include <mesh_common/vtk/DistributedVTKMeshWriter.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+template<typename MeshType, typename T>
+class VertexDataSource {
+public:
+   typedef typename mesh::DistributedVTKMeshWriter<MeshType>::Vertices Vertices;
+
+   explicit VertexDataSource(const std::string & name) : name_(name) {}
+
+   const std::string & name() { return name_; }
+   virtual uint_t numComponents() = 0;
+   virtual void getData(const MeshType &, const Vertices &vertices, std::vector<T> &data,
+         const ParticleIdxVertexPropertyManager<MeshType> & vertexToParticleIdxManager,
+         shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps) = 0;
+
+   virtual ~VertexDataSource() {}
+
+protected:
+   std::string name_;
+};
+
+}
+}
\ No newline at end of file
diff --git a/src/mesa_pd/vtk/ConvexPolyhedron/tesselation/ConvexPolyhedronTesselation.h b/src/mesa_pd/vtk/ConvexPolyhedron/tesselation/ConvexPolyhedronTesselation.h
new file mode 100644
index 0000000000000000000000000000000000000000..87bed413da182d220fed5ea869b29328e046d8c6
--- /dev/null
+++ b/src/mesa_pd/vtk/ConvexPolyhedron/tesselation/ConvexPolyhedronTesselation.h
@@ -0,0 +1,81 @@
+//======================================================================================================================
+//
+//  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 ConvexPolyhedronTesselation.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/data/shape/ConvexPolyhedron.h>
+#include <mesh_common/TriangleMeshes.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+template<typename MeshType>
+void tesselate(const walberla::mesa_pd::data::ConvexPolyhedron& shape,
+      const walberla::mesa_pd::data::ParticleStorage::Particle& particle,
+      shared_ptr<MeshType> mesh, std::vector<typename MeshType::VertexHandle>& newVertices,
+      std::vector<typename MeshType::FaceHandle>& newFaces) {
+   // partially taken from pe implementation
+   typedef typename MeshType::Scalar Scalar;
+   typedef typename MeshType::Point Point;
+   typedef typename MeshType::VertexHandle VertexHandle;
+
+   const auto& polyhedronMesh = shape.getMesh();
+   const auto& position = particle.getPosition();
+
+   newVertices.reserve(polyhedronMesh.n_vertices());
+   for (auto vh: polyhedronMesh.vertices()) {
+      // p0: point in mesh coordinate system (and thus, in particle COS, as they have to overlap at the origin)
+      // p1: point in world coordinate system
+      auto p0 = mesh::toWalberla(polyhedronMesh.point(vh));
+      Vector3<real_t> p1 = position + particle.getRotation().getMatrix()*p0;
+
+      Point p2(numeric_cast<Scalar>(p1[0]),
+            numeric_cast<Scalar>(p1[1]),
+            numeric_cast<Scalar>(p1[2]));
+
+      newVertices.push_back(mesh->add_vertex(p2));
+   }
+
+   newFaces.reserve( polyhedronMesh.n_faces() );
+   for (auto fh: polyhedronMesh.faces()) {
+      auto it = polyhedronMesh.cfv_ccwbegin(fh);
+      VertexHandle v0 = newVertices[ size_t( (it++)->idx() ) ];
+      VertexHandle v1 = newVertices[ size_t( (it++)->idx() ) ];
+      VertexHandle v2 = newVertices[ size_t( (it++)->idx() ) ];
+      WALBERLA_ASSERT_EQUAL(it, polyhedronMesh.cfv_ccwend(fh) );
+      newFaces.push_back( mesh->add_face(v0, v1, v2) );
+   }
+
+   if (mesh->has_face_normals()) {
+      for (auto fh: newVertices) {
+         mesh->update_normal(fh);
+      }
+   }
+
+   if (mesh->has_vertex_normals()) {
+      for (auto vh: newVertices) {
+         mesh->update_normal(vh);
+      }
+   }
+}
+
+}
+}
\ No newline at end of file
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 9883532626fc27dc04a61483b919caf2494c8e28..e53281ee904ac96c4570f5ae73ba4df065300fa2 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -5,7 +5,7 @@
 ###################################################################################################
 
 waLBerla_add_module( DEPENDS blockforest boundary core domain_decomposition
-                             python_coupling field geometry pe stencil mesa_pd mesh_common
+                             python_coupling field geometry pe stencil mesh_common
                      BUILD_ONLY_IF_FOUND OpenMesh)
 
 ###################################################################################################
diff --git a/src/mesh_common/DistanceComputations.h b/src/mesh_common/DistanceComputations.h
index 3a1199b480de58e5b321730060be3e32495eeff6..997546f1e7836fc81e56092e5aef20ac86b44d29 100644
--- a/src/mesh_common/DistanceComputations.h
+++ b/src/mesh_common/DistanceComputations.h
@@ -200,8 +200,8 @@ protected:
 template< typename MeshType >
 void TriangleDistance<MeshType>::computeNormals()
 {
-   // J. Bærentzen and H. Aanæs. Signed distance computation using the angle weighted pseudonormal.
-   // Visualization and Computer Graphics, IEEE Transactions on, 11(3):243–253, 2005.
+   // J. B�rentzen and H. Aan�s. Signed distance computation using the angle weighted pseudonormal.
+   // Visualization and Computer Graphics, IEEE Transactions on, 11(3):243�253, 2005.
 
    static_assert( MeshType::IsTriMesh == 1, "computeNormals only works with triangular meshes!" );
 
diff --git a/src/mesh_common/MeshIO.h b/src/mesh_common/MeshIO.h
index 63d8f531f55ff2b6a1f91f2b93dd0ca5f56ca987..0a6341961b68a2dfb3007f042beac83a90f2b995 100644
--- a/src/mesh_common/MeshIO.h
+++ b/src/mesh_common/MeshIO.h
@@ -16,6 +16,7 @@
 //! \file MeshIO.h
 //! \ingroup mesh
 //! \author Christian Godenschwager <christian.godenschwager@fau.de>
+//! \author Lukas Werner
 //
 //======================================================================================================================
 
@@ -42,6 +43,30 @@
 namespace walberla {
 namespace mesh {
 
+/**
+ * \brief Reads a mesh from a generic input stream.
+ *
+ * \tparam MeshType     The type of the OpenMesh
+ *
+ * \param inputStream   The input stream from which the mesh should be read
+ * \param mesh          The mesh data structure to be written to
+ * \param extension     The mesh file's extension
+ *
+ * \return Whether the read operation was successful.
+ */
+template< typename MeshType >
+bool readFromStream( std::istream & inputStream, MeshType & mesh, const std::string & extension,
+                     bool binaryFile = false )
+{
+   OpenMesh::IO::Options options;
+   if( mesh.has_face_colors() )
+      options += OpenMesh::IO::Options::FaceColor;
+   if ( binaryFile )
+      options += OpenMesh::IO::Options::Binary;
+   if( mesh.has_vertex_colors() )
+      options += OpenMesh::IO::Options::VertexColor;
+   return OpenMesh::IO::read_mesh( mesh, inputStream, extension, options );
+}
 
 /**
 * \brief Loads an OpenMesh in parallel
@@ -86,16 +111,41 @@ void readAndBroadcast( const std::string & filename, MeshType & mesh, bool binar
    if ( binaryFile )
       iss = std::istringstream( str, std::ifstream::in | std::ifstream::binary );
 
-   OpenMesh::IO::Options options;
-   if( mesh.has_face_colors() )
-      options += OpenMesh::IO::Options::FaceColor;
-   if ( binaryFile )
-      options += OpenMesh::IO::Options::Binary;
-   if( mesh.has_vertex_colors() )
-      options += OpenMesh::IO::Options::VertexColor;
-   if( !OpenMesh::IO::read_mesh( mesh, iss, extension, options ) )
+   if (!readFromStream<MeshType>(iss, mesh, extension, binaryFile)) {
+      WALBERLA_ABORT( "Error while reading file \"" << filename << "\"!" );
+   }
+}
+
+/**
+ * \brief Read a mesh from a file.
+ * \attention If reading in parallel, readAndBroadcast should be preferred!
+ *
+ * \tparam MeshType     The type of the OpenMesh
+ *
+ * \param filename      Filename of the mesh to be loaded
+ * \param mesh          The mesh data structure to be written to
+ */
+template< typename MeshType >
+void readFromFile( const std::string & filename, MeshType & mesh, bool binaryFile = false )
+{
+   if (!filesystem::exists(filename)) {
+      WALBERLA_ABORT( "The mesh file \"" << filename << "\" does not exist!" );
+   }
+
+   std::string extension = filesystem::path(filename).extension().string();
+
+   std::ios_base::openmode openMode = std::ifstream::in;
+   if (binaryFile) {
+      openMode |= std::ifstream::binary;
+   }
+   std::ifstream inputFileStream = std::ifstream(filename, openMode);
+
+   if (!readFromStream<MeshType>(inputFileStream, mesh, extension, binaryFile)) {
       WALBERLA_ABORT( "Error while reading file \"" << filename << "\"!" );
+   }
+   inputFileStream.close();
 }
 
+
 } // namespace mesh
 } // namespace walberla
diff --git a/src/mesh_common/MeshOperations.h b/src/mesh_common/MeshOperations.h
index d659d0d48a3b9a8b3a9a454861b8be2855f49381..b870817508fe2e00186fc7b0f5dfea7e944e9b70 100644
--- a/src/mesh_common/MeshOperations.h
+++ b/src/mesh_common/MeshOperations.h
@@ -253,6 +253,7 @@ typename MeshType::Point computeCentroid( const MeshType & mesh )
 template< typename MeshType >
 Matrix3<typename MeshType::Scalar> computeInertiaTensor( const MeshType & mesh )
 {
+   // Inertia tensor is calculated relative to the origin of the meshes coordinate system!
    static_assert( MeshType::IsTriMesh == 1, "computeInertiaTensor only works with triangular meshes!" );
 
    typedef typename MeshType::Point  Point;
@@ -315,6 +316,100 @@ typename MeshType::Point computeCentroid( const MeshType & mesh, const typename
    return centroid;
 }
 
+/**
+ * \brief Computes all mass properties (mass, centroid, inertia tensor at once).
+ *
+ * This function computes the mass, centroid and the inertia tensor relative to the calculated centroid.
+ * Source: https://www.cs.upc.edu/~virtual/SGI/docs/3.%20Further%20Reading/Polyhedral%20Mass%20Properties%20Revisited.pdf
+ *
+ * \tparam MeshType
+ * \param mesh      Triangular input mesh.
+ * \param density   Density of the mesh.
+ * \param centroid  Output centroid point.
+ * \param inertiaTensor Output inertia matrix.
+ * \param mass      Output mass.
+ * \attention The inertia tensor is computed relative to the centroid.
+ */
+template< typename MeshType >
+void computeMassProperties(const MeshType & mesh, typename MeshType::Scalar density,
+      Vector3<typename MeshType::Scalar>& centroid, Matrix3<typename MeshType::Scalar>& inertiaTensor,
+      typename MeshType::Scalar& mass) {
+   static_assert( MeshType::IsTriMesh == 1, "computeMassProperties only works with triangular meshes!" );
+
+   typedef typename MeshType::Point Point;
+   typedef typename MeshType::Scalar Scalar;
+
+   const Scalar mult[10] = {Scalar(1)/Scalar(6),
+                            Scalar(1)/Scalar(24), Scalar(1)/Scalar(24), Scalar(1)/Scalar(24),
+                            Scalar(1)/Scalar(60), Scalar(1)/Scalar(60), Scalar(1)/Scalar(60),
+                            Scalar(1)/Scalar(120), Scalar(1)/Scalar(120), Scalar(1)/Scalar(120)};
+
+   Scalar intg[10] = {0,0,0,0,0,0,0,0,0,0};
+
+   auto subExpr = [](Scalar& w0, Scalar& w1, Scalar& w2,
+         Scalar& f1, Scalar& f2, Scalar& f3,
+         Scalar& g0, Scalar& g1, Scalar& g2){
+      Scalar temp0 = w0+w1;
+      f1 = temp0 + w2;
+      Scalar temp1 = w0*w0;
+      Scalar temp2 = temp1 + w1*temp0;
+      f2 = temp2 + w2*f1;
+      f3 = w0*temp1 + w1*temp2 + w2*f2;
+      g0 = f2 + w0*(f1+w0);
+      g1 = f2 + w1*(f1+w1);
+      g2 = f2 + w2*(f1+w2);
+   };
+
+   Point v0, v1, v2;
+   for (auto fIt = mesh.faces_begin(); fIt != mesh.faces_end(); ++fIt) {
+      getVertexPositions(mesh, *fIt, v0, v1, v2);
+
+      Scalar a1 = v1[0]-v0[0];
+      Scalar b1 = v1[1]-v0[1];
+      Scalar c1 = v1[2]-v0[2];
+      Scalar a2 = v2[0]-v0[0];
+      Scalar b2 = v2[1]-v0[1];
+      Scalar c2 = v2[2]-v0[2];
+
+      Scalar d0 = b1*c2 - b2*c1;
+      Scalar d1 = a2*c1 - a1*c2;
+      Scalar d2 = a1*b2 - a2*b1;
+
+      Scalar f1x, f2x, f3x, g0x, g1x, g2x;
+      subExpr(v0[0], v1[0], v2[0], f1x, f2x, f3x, g0x, g1x, g2x);
+      Scalar f1y, f2y, f3y, g0y, g1y, g2y;
+      subExpr(v0[1], v1[1], v2[1], f1y, f2y, f3y, g0y, g1y, g2y);
+      Scalar f1z, f2z, f3z, g0z, g1z, g2z;
+      subExpr(v0[2], v1[2], v2[2], f1z, f2z, f3z, g0z, g1z, g2z);
+
+      intg[0] += d0*f1x;
+      intg[1] += d0*f2x; intg[2] += d1*f2y; intg[3] += d2*f2z;
+      intg[4] += d0*f3x; intg[5] += d1*f3y; intg[6] += d2*f3z;
+      intg[7] += d0*(v0[1]*g0x + v1[1]*g1x + v2[1]*g2x);
+      intg[8] += d1*(v0[2]*g0y + v1[2]*g1y + v2[2]*g2y);
+      intg[9] += d2*(v0[0]*g0z + v1[0]*g1z + v2[0]*g2z);
+   }
+
+   for (uint_t i = 0; i < 10; ++i) {
+      intg[i] *= mult[i];
+   }
+
+   mass = intg[0];
+
+   centroid[0] = intg[1] / mass;
+   centroid[1] = intg[2] / mass;
+   centroid[2] = intg[3] / mass;
+
+   inertiaTensor[0] = density * (intg[5] + intg[6] - mass*(centroid[1]*centroid[1] + centroid[2]*centroid[2])); //xx
+   inertiaTensor[4] = density * (intg[4] + intg[6] - mass*(centroid[2]*centroid[2] + centroid[0]*centroid[0])); // yy
+   inertiaTensor[8] = density * (intg[4] + intg[5] - mass*(centroid[0]*centroid[0] + centroid[1]*centroid[1])); // zz
+   inertiaTensor[1] = density * (-(intg[7] - mass * centroid[0]*centroid[1]));
+   inertiaTensor[5] = density * (-(intg[8] - mass * centroid[1]*centroid[2]));
+   inertiaTensor[2] = density * (-(intg[9] - mass * centroid[2]*centroid[0]));
+
+   mass *= density;
+}
+
 
 template< typename MeshType, typename InputIterator >
 std::vector< typename MeshType::VertexHandle > findConnectedVertices( const MeshType & mesh, const InputIterator facesBegin, const InputIterator facesEnd )
diff --git a/src/waLBerlaDefinitions.in.h b/src/waLBerlaDefinitions.in.h
index 654eb55ab235a5a583099cdbcfedb46e1844297a..360eea932f93c3c07f0dce034617142f6bf82cac 100644
--- a/src/waLBerlaDefinitions.in.h
+++ b/src/waLBerlaDefinitions.in.h
@@ -30,6 +30,7 @@
 #cmakedefine WALBERLA_BUILD_WITH_FFT
 
 #cmakedefine WALBERLA_BUILD_WITH_OPENMESH
+#cmakedefine WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE
 
 #cmakedefine WALBERLA_BUILD_WITH_CUDA
 #cmakedefine WALBERLA_CUDA_NVTX_AVAILABLE
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index 6397eb6f318b2ed08c38b85a486fa4e309809aaf..3742a2ef2dcabdd68d50fe0c436fd0b21d360435 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -188,3 +188,13 @@ waLBerla_execute_test( NAME   MESA_PD_Sorting )
 waLBerla_compile_test( NAME   MESA_PD_VTK_Outputs FILES vtk/VTKOutputs.cpp DEPENDS blockforest core vtk )
 waLBerla_execute_test( NAME   MESA_PD_VTK_Outputs PROCESSES 8 )
 
+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_Data_ConvexPolyhedron FILES data/ConvexPolyhedron.cpp DEPENDS core mesa_pd mesh_common )
+    waLBerla_execute_test( NAME   MESA_PD_Data_ConvexPolyhedron )
+
+    waLBerla_compile_test( NAME   MESA_PD_COLLISIONDETECTION_ConvexPolyhedron_GJK_EPA FILES collision_detection/ConvexPolyhedron_GJK_EPA.cpp DEPENDS core mesa_pd mesh_common )
+    waLBerla_execute_test( NAME   MESA_PD_COLLISIONDETECTION_ConvexPolyhedron_GJK_EPA )
+endif()
\ No newline at end of file
diff --git a/tests/mesa_pd/collision_detection/ConvexPolyhedron_GJK_EPA.cpp b/tests/mesa_pd/collision_detection/ConvexPolyhedron_GJK_EPA.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..037073c9cdbba45429d240388e8838d3ec7775bf
--- /dev/null
+++ b/tests/mesa_pd/collision_detection/ConvexPolyhedron_GJK_EPA.cpp
@@ -0,0 +1,181 @@
+//======================================================================================================================
+//
+//  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 CollisionTest.cpp
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#include "core/debug/TestSubsystem.h"
+
+#include <core/Environment.h>
+#include <core/all.h>
+#include <core/logging/Logging.h>
+#include <mesa_pd/collision_detection/GeneralContactDetection.h>
+#include <mesa_pd/data/Flags.h>
+#include <mesa_pd/data/shape/Ellipsoid.h>
+#include <mesa_pd/data/shape/Sphere.h>
+
+#include "mesa_pd/data/shape/ConvexPolyhedron.h"
+#include "mesh_common/MeshIO.h"
+#include "mesh_common/TriangleMeshes.h"
+
+namespace walberla {
+namespace mesa_pd {
+
+using namespace walberla::mesa_pd::collision_detection;
+using namespace walberla::mesa_pd::collision_detection::analytic;
+using namespace walberla::mesa_pd::data;
+
+bool gjkEPAcollideHybrid(Support &geom1, Support &geom2, Vec3& normal, Vec3& contactPoint, real_t& penetrationDepth) {
+   // This function is stolen from the general mesapd GJK-EPA test.
+
+   // For more information on hybrid GJK/EPA see page 166 in "Collision Detecton in Interactive 3D
+   // Environments" by Gino van den Bergen.
+
+   //1. Run GJK with considerably enlarged objects.
+   real_t margin = real_t(1e-6);
+   GJK gjk;
+   if(gjk.doGJKmargin(geom1, geom2, margin)){
+      //2. If collision is possible perform EPA.
+      //std::cerr << "Peforming EPA.";
+      EPA epa;
+      epa.useSphereOptimization( true );
+      return epa.doEPAmargin(geom1, geom2, gjk, normal, contactPoint, penetrationDepth, margin);
+   } else {
+      return false;
+   }
+}
+
+/**
+ * \brief Compare the collision results of a normal box with the results of a mesh resembling this box.
+ * \param testSupport The particle with which the box and the mesh are colliding.
+ * \param convPoly Convex polyhedron, contains a mesh resembling the box.
+ * \param box Box, reference for convex polyhedron.
+ */
+void runComparisonTest(Support& testSupport, Support& convPoly, Support& box) {
+   Vec3 polyNormal;
+   Vec3 polyContactPoint;
+   real_t polyPenetrationDepth;
+   bool polyCollided = gjkEPAcollideHybrid(convPoly, testSupport, polyNormal, polyContactPoint, polyPenetrationDepth);
+   //WALBERLA_LOG_INFO("conv poly: " << polyCollided << ", " << polyNormal << ", " << polyContactPoint << ", " << polyPenetrationDepth);
+
+   Vec3 boxNormal;
+   Vec3 boxContactPoint;
+   real_t boxPenetrationDepth;
+   bool boxCollided = gjkEPAcollideHybrid(box, testSupport, boxNormal, boxContactPoint, boxPenetrationDepth);
+   //WALBERLA_LOG_INFO("box:       " << boxCollided << ", " << boxNormal << ", " << boxContactPoint << ", " << boxPenetrationDepth);
+
+   WALBERLA_CHECK(polyCollided == boxCollided);
+   if (polyCollided) {
+      auto eps = real_t(1e-3);
+      WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(polyNormal, boxNormal, real_t(1e-2));
+      WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(polyContactPoint, boxContactPoint, eps);
+      WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(polyPenetrationDepth, boxPenetrationDepth, eps);
+   }
+}
+
+/** Test the GJK-EPA implementation for collision detection with convex polyhedrons vs. multiple other shapes.
+ */
+int main(int argc, char** argv){
+   walberla::debug::enterTestMode();
+   walberla::mpi::Environment env(argc, argv);
+   walberla::mpi::MPIManager::instance()->useWorldComm();
+
+   if (std::is_same<walberla::real_t, float>::value) {
+      WALBERLA_LOG_WARNING("waLBerla build in sp mode: skipping test due to low precision");
+      return EXIT_SUCCESS;
+   }
+
+   WALBERLA_LOG_INFO("ConvexPolyhedron <-> Sphere vs. Box <-> Sphere");
+   mesh::TriangleMesh cubeMesh;
+   std::string cubeInputFile = "../mesh/cube.obj";
+   mesh::readFromFile< mesh::TriangleMesh >(cubeInputFile, cubeMesh);
+   mesh::translate(cubeMesh, -mesh::toWalberla(mesh::computeCentroid(cubeMesh)));
+
+   ConvexPolyhedron convPoly_(cubeMesh);
+   Support convPoly(Vec3(0), Rot3(), convPoly_);
+
+   Box box_(Vec3(1));
+   Support box(Vec3(0), Rot3(), box_);
+
+   Sphere sphere_(real_t(0.5));
+   Support sphere(Vec3(real_t(0.5),real_t(0.5),real_t(0.9)), Rot3(), sphere_);
+   runComparisonTest(sphere, convPoly, box);
+
+   /// Fuzzy Testing
+
+   AABB collisionTestAABB(mesa_pd::Vec3(real_t(-1)), mesa_pd::Vec3(real_t(1)));
+   std::mt19937 rng;
+
+   WALBERLA_LOG_INFO("Fuzzy ConvexPolyhedron <-> Sphere vs. Box <-> Sphere");
+   for(uint_t i = 0; i < 500; ++i) {
+      auto rndPoint = collisionTestAABB.randomPoint(rng);
+
+      Sphere fuzzSphere_(real_t(0.5));
+      Support fuzzSphere(rndPoint, Rot3(), sphere_);
+
+      runComparisonTest(fuzzSphere, convPoly, box);
+   }
+
+   WALBERLA_LOG_INFO("Fuzzy ConvexPolyhedron <-> Ellipsoid vs. Box <-> Ellipsoid");
+   for(uint_t i = 0; i < 500; ++i) {
+      auto rndPoint = collisionTestAABB.randomPoint(rng);
+
+      Ellipsoid el_(rndPoint);
+      Support el(Vec3(real_t(1)), Rot3(), el_);
+
+      el.rot_.rotate(Vec3(math::realRandom(), math::realRandom(), math::realRandom()));
+
+      runComparisonTest(el, convPoly, box);
+   }
+
+   WALBERLA_LOG_INFO("Fuzzy ConvexPolyhedron <-> Turned Box vs. Box <-> Turned Box");
+   for(uint_t i = 0; i < 500; ++i) {
+      auto rndPoint = collisionTestAABB.randomPoint(rng);
+
+      Box fuzzBox_(Vec3(real_t(1),real_t(0.5),real_t(0.5)));
+      Support fuzzBox(rndPoint, Rot3(), fuzzBox_);
+
+      fuzzBox.rot_.rotate(Vec3(math::realRandom(), math::realRandom(), math::realRandom()));
+
+      runComparisonTest(fuzzBox, convPoly, box);
+   }
+
+   WALBERLA_LOG_INFO("Fuzzy Turned ConvexPolyhedron <-> Box vs. Turned Box <-> Box");
+   for(uint_t i = 0; i < 500; ++i) {
+      auto rndPoint = collisionTestAABB.randomPoint(rng);
+
+      Box fuzzBox_(Vec3(real_t(1),real_t(0.5),real_t(0.5)));
+      Support fuzzBox(rndPoint, Rot3(), fuzzBox_);
+
+      // rotate all the things
+      fuzzBox.rot_.rotate(Vec3(math::realRandom(), math::realRandom(), math::realRandom()));
+      Vec3 rotPhi(math::realRandom(), math::realRandom(), math::realRandom());
+      convPoly.rot_.rotate(rotPhi);
+      box.rot_.rotate(rotPhi);
+
+      runComparisonTest(fuzzBox, convPoly, box);
+   }
+
+   return EXIT_SUCCESS;
+}
+
+}
+}
+
+int main( int argc, char* argv[] ) {
+   return walberla::mesa_pd::main( argc, argv );
+}
\ No newline at end of file
diff --git a/tests/mesa_pd/data/ConvexPolyhedron.cpp b/tests/mesa_pd/data/ConvexPolyhedron.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ff9f2ab7690a45d5a8ceadc9ab853925cd3c6ea9
--- /dev/null
+++ b/tests/mesa_pd/data/ConvexPolyhedron.cpp
@@ -0,0 +1,182 @@
+//======================================================================================================================
+//
+//  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 MeshMesapdConvexPolyhedronTest.cpp
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#include "mesa_pd/data/shape/ConvexPolyhedron.h"
+
+#include "core/debug/TestSubsystem.h"
+#include "core/mpi/all.h"
+
+#include "mesa_pd/data/DataTypes.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/mpi/ShapePackUnpack.h"
+#include "mesh_common/MeshIO.h"
+#include "mesh_common/MeshOperations.h"
+#include "mesh_common/PolyMeshes.h"
+
+namespace walberla {
+
+int main( int argc, char** argv )
+{
+   walberla::debug::enterTestMode();
+   mpi::Environment env(argc, argv);
+   mpi::MPIManager::instance()->useWorldComm();
+
+   /// Load Mesh File
+   WALBERLA_LOG_INFO("LOAD MESH");
+   // Ideally this should be done on only one process and
+   // communicated to other processes, if or when required.
+   mesh::TriangleMesh bunnyMesh;
+   std::string meshInputFile = "../mesh/bunny.obj";
+   mesh::readFromFile<mesh::TriangleMesh>(meshInputFile, bunnyMesh);
+   WALBERLA_LOG_INFO("Read mesh file: " << meshInputFile << " (" << bunnyMesh.n_vertices() << " vertices, "
+      << bunnyMesh.n_faces() << " faces, volume = " << mesh::computeVolume(bunnyMesh) << ")" );
+   mesh::translate(bunnyMesh, -mesh::toWalberla(mesh::computeCentroid(bunnyMesh)));
+
+   /// MESAPD Data
+   WALBERLA_LOG_INFO("CREATE CONVEX POLYHEDRON");
+   auto ps = std::make_shared<mesa_pd::data::ParticleStorage>(3);
+   auto ss = std::make_shared<mesa_pd::data::ShapeStorage>();
+   auto ac = walberla::make_shared<mesa_pd::data::ParticleAccessorWithShape>(ps, ss);
+   auto bunnyShapeID = ss->create<mesa_pd::data::ConvexPolyhedron>(bunnyMesh);
+
+   auto bunnyParticle = ps->create();
+   bunnyParticle->setShapeID(bunnyShapeID);
+   auto mesh1ParticleUid = bunnyParticle->getUid();
+   auto mesh1ParticleIdx = ac->uidToIdx(mesh1ParticleUid);
+
+   auto mesh1ParticleShape = dynamic_cast<mesa_pd::data::ConvexPolyhedron*>(ac->getShape(mesh1ParticleIdx));
+
+   WALBERLA_CHECK_FLOAT_EQUAL(ac->getShape(mesh1ParticleIdx)->getVolume(), mesh::computeVolume(bunnyMesh));
+   WALBERLA_CHECK_EQUAL(mesh1ParticleShape->getMesh().n_vertices(), bunnyMesh.n_vertices());
+   WALBERLA_CHECK_EQUAL(mesh1ParticleShape->getMesh().n_faces(), bunnyMesh.n_faces());
+
+   /// Check: Pack - Unpack
+   WALBERLA_LOG_INFO("PACK - UNPACK");
+
+   shared_ptr<mesa_pd::data::BaseShape> bs0 = make_shared<mesa_pd::data::ConvexPolyhedron>(bunnyMesh);
+   mesh1ParticleShape->updateMassAndInertia(real_t(1));
+   std::shared_ptr<mesa_pd::data::BaseShape> bs1 = nullptr;
+
+   WALBERLA_LOG_INFO("Packing mesh shape");
+   mpi::SendBuffer sb;
+   sb << bs0;
+   WALBERLA_LOG_INFO("Unpacking mesh shape");
+   mpi::RecvBuffer rb(sb);
+   rb >> bs1;
+
+   //WALBERLA_CHECK_EQUAL(bs0->getShapeType(), mesa_pd::data::ConvexPolyhedron::SHAPE_TYPE);
+   //WALBERLA_CHECK_EQUAL(bs1->getShapeType(), mesa_pd::data::ConvexPolyhedron::SHAPE_TYPE);
+   WALBERLA_CHECK_IDENTICAL(bs0->getMass(), bs1->getMass());
+   WALBERLA_CHECK_IDENTICAL(bs0->getInvMass(), bs1->getInvMass());
+   WALBERLA_CHECK_IDENTICAL(bs0->getInertiaBF(), bs1->getInertiaBF());
+   WALBERLA_CHECK_IDENTICAL(bs0->getInvInertiaBF(), bs1->getInvInertiaBF());
+
+   auto bm0 = dynamic_cast<mesa_pd::data::ConvexPolyhedron*>(bs0.get());
+   auto bm1 = dynamic_cast<mesa_pd::data::ConvexPolyhedron*>(bs1.get());
+   WALBERLA_CHECK_EQUAL(bm0->getMesh().n_vertices(), bm1->getMesh().n_vertices());
+   WALBERLA_CHECK_EQUAL(bm0->getMesh().n_faces(), bm1->getMesh().n_faces());
+
+   /// Volume and Inertia
+   WALBERLA_LOG_INFO("VOLUME AND INERTIA");
+
+   mesh::TriangleMesh cubeMesh;
+   std::string cubeMeshInputFile = "../mesh/cube.obj";
+   mesh::readFromFile<mesh::TriangleMesh>(cubeMeshInputFile, cubeMesh);
+   mesh::translate(cubeMesh, -mesh::toWalberla(mesh::computeCentroid(cubeMesh)));
+
+   real_t cubeVolume = mesh::computeVolume(cubeMesh);
+   real_t cubeSideLen = std::cbrt(cubeVolume);
+   auto cubeDensity = real_t(123);
+   mesa_pd::data::Box boxCubeShape((Vector3<real_t>(cubeSideLen)));
+   boxCubeShape.updateMassAndInertia(cubeDensity);
+
+   mesa_pd::data::ConvexPolyhedron meshCubeShape(cubeMesh);
+   meshCubeShape.updateMassAndInertia(cubeDensity);
+
+   // Test mass properties of MESAPD convexpolyhedron body
+
+   WALBERLA_CHECK_FLOAT_EQUAL(boxCubeShape.getInertiaBF(), meshCubeShape.getInertiaBF());
+   WALBERLA_CHECK_FLOAT_EQUAL(boxCubeShape.getVolume(), meshCubeShape.getVolume());
+   WALBERLA_CHECK_FLOAT_EQUAL(boxCubeShape.getMass(), meshCubeShape.getMass());
+   WALBERLA_CHECK_FLOAT_EQUAL(boxCubeShape.getInvMass(), meshCubeShape.getInvMass());
+
+   // Test new mesh ops function to calculate all mass properties at once
+
+   Matrix3<real_t> massPropInertia;
+   Vector3<real_t> massPropCentroid;
+   real_t massPropMass;
+
+   mesh::computeMassProperties(cubeMesh, cubeDensity, massPropCentroid, massPropInertia, massPropMass);
+
+   WALBERLA_CHECK_FLOAT_EQUAL(boxCubeShape.getInertiaBF(), massPropInertia);
+   WALBERLA_CHECK_FLOAT_EQUAL(boxCubeShape.getMass(), massPropMass);
+   WALBERLA_CHECK_FLOAT_EQUAL(massPropCentroid, Vector3<real_t>(real_t(0)))
+
+   // Test new mesh ops inertia calculation against old one
+
+   auto oldInertia = mesh::computeInertiaTensor(cubeMesh)*cubeDensity;
+   WALBERLA_CHECK_FLOAT_EQUAL(massPropInertia, oldInertia);
+
+   /// Interaction radius
+   // for a cube: space diagonal length = side lengths * sqrt(3) -> bounding sphere radius = space diagonal / 2
+   WALBERLA_CHECK_FLOAT_EQUAL(meshCubeShape.getBoundingSphereRadius(), real_t(1) * std::sqrt(real_t(3)) * real_t(0.5));
+
+   auto bunnyShape = dynamic_cast<mesa_pd::data::ConvexPolyhedron*>(ss->shapes[bunnyShapeID].get());
+   real_t bunnyRadius = bunnyShape->getBoundingSphereRadius();
+   real_t maxSqRadius(0);
+   for(auto vh : bunnyMesh.vertices()) {
+      auto v = mesh::toWalberla(bunnyMesh.point(vh));
+      auto centroidToVSqr = v.sqrLength();
+
+      if (centroidToVSqr > maxSqRadius) {
+         maxSqRadius = centroidToVSqr;
+      }
+   }
+   WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(bunnyRadius, std::sqrt(maxSqRadius), real_t(1e-4));
+
+   /// Support
+   WALBERLA_LOG_INFO("SUPPORT");
+
+   for (int x = -1; x <= 1; x+=2) {
+      for (int y = -1; y <= 1; y+=2) {
+         for (int z = -1; z <= 1; z+=2) {
+            Vector3<real_t> d((real_t(x)), real_t(y), real_t(z));
+            WALBERLA_CHECK_FLOAT_EQUAL(boxCubeShape.support(d), meshCubeShape.support(d));
+         }
+      }
+   }
+
+   AABB supportTestAABB(mesa_pd::Vec3(real_t(-1)), mesa_pd::Vec3(real_t(1)));
+   std::mt19937 rng;
+   for(uint_t i = 0; i < 500; ++i) {
+      auto rndPoint = supportTestAABB.randomPoint(rng);
+      WALBERLA_CHECK_FLOAT_EQUAL(boxCubeShape.support(rndPoint), meshCubeShape.support(rndPoint));
+   }
+
+   return EXIT_SUCCESS;
+}
+}
+
+int main( int argc, char* argv[] )
+{
+   return walberla::main( argc, argv );
+}
\ No newline at end of file
diff --git a/tests/mesa_pd/vtk/ConvexPolyhedronVTKOutput.cpp b/tests/mesa_pd/vtk/ConvexPolyhedronVTKOutput.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..124493ca93043b9aa6955b5f3444fb4690845e99
--- /dev/null
+++ b/tests/mesa_pd/vtk/ConvexPolyhedronVTKOutput.cpp
@@ -0,0 +1,131 @@
+//======================================================================================================================
+//
+//  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 VTKOutput.cpp
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#include "blockforest/all.h"
+
+#include "core/debug/TestSubsystem.h"
+#include "core/math/Rot3.h"
+#include "core/mpi/all.h"
+
+#include "vtk/all.h"
+
+#include <mesa_pd/domain/BlockForestDomain.h>
+#include <mesa_pd/mpi/ShapePackUnpack.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/MeshParticleVTKOutput.h>
+#include <mesa_pd/vtk/ConvexPolyhedron/data_sources/SurfaceVelocityVertexDataSource.h>
+#include <mesh_common/PolyMeshes.h>
+
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+#include "mesh_common/MeshIO.h"
+#include "mesh_common/MeshOperations.h"
+
+namespace walberla {
+
+using namespace walberla::mesa_pd;
+
+int main(int argc, char** argv) {
+   walberla::debug::enterTestMode();
+   mpi::Environment env(argc, argv);
+   mpi::MPIManager::instance()->useWorldComm();
+
+   std::string vtk_out = "vtk_out";
+
+   /// BlockForest
+   Vector3< real_t > minCorner(real_t(0), real_t(0), real_t(0));
+   Vector3< real_t > maxCorner(real_t(4), real_t(4), real_t(4));
+   math::AABB domainSize(minCorner, maxCorner);
+   Vector3< bool > periodicity(true, true, true);
+   Vector3< uint_t > numBlocks(2, 2, 1);
+   shared_ptr< BlockForest > forest = blockforest::createBlockForest(domainSize, numBlocks, periodicity);
+   auto domain = std::make_shared< mesa_pd::domain::BlockForestDomain >(forest);
+
+   WALBERLA_LOG_INFO(domain->getNumLocalAABBs() << ": " << domain->getUnionOfLocalAABBs());
+
+   /// MESAPD Data
+   WALBERLA_LOG_INFO_ON_ROOT("CREATE MESAPD DATA STRUCTURES");
+   auto ps = std::make_shared< data::ParticleStorage >(4);
+   auto ss = std::make_shared< data::ShapeStorage >();
+   auto ac = walberla::make_shared< data::ParticleAccessorWithShape >(ps, ss);
+
+   /// Load Mesh File
+   WALBERLA_LOG_INFO_ON_ROOT("LOAD MESH");
+   // Ideally this should be done on only one process and
+   // communicated to other processes, if or when required.
+   mesh::TriangleMesh cubeMesh;
+   std::string cubeMeshInputFile = "../mesh/cube.obj";
+   mesh::readAndBroadcast< mesh::TriangleMesh >(cubeMeshInputFile, cubeMesh);
+   WALBERLA_LOG_INFO_ON_ROOT("Read mesh file: " << cubeMeshInputFile << " (" << cubeMesh.n_vertices() << " vertices, "
+                                                << cubeMesh.n_faces()
+                                                << " faces, volume = " << mesh::computeVolume(cubeMesh) << ")");
+   mesh::translate(cubeMesh, mesh::toWalberla(-mesh::computeCentroid(cubeMesh)));
+
+   /// Mesh Shape
+   auto cubeShapeID = ss->create< data::ConvexPolyhedron >(cubeMesh);
+
+   /// Mesh Particles
+   math::Rot3< real_t > cubeRotation(Vector3< real_t >(real_t(M_PI) / real_t(4), real_t(0), real_t(0)));
+   for (uint_t x = 0; x <= 1; ++x) {
+      for (uint_t y = 0; y <= 1; ++y) {
+         for (uint_t z = 0; z <= 1; ++z) {
+            // position of new mesh cube slightly shifted inwards
+            Vector3< real_t > pos((real_t(0.05) + real_t(x) * maxCorner[0] - real_t(x) * real_t(0.1)),
+                                  real_t(0.05) + real_t(y) * maxCorner[1] - real_t(y) * real_t(0.1),
+                                  real_t(0.05) + real_t(z) * maxCorner[2] - real_t(z) * real_t(0.1));
+
+            if (domain->isContainedInProcessSubdomain(uint_c(mpi::MPIManager::instance()->rank()), pos)) {
+               auto rotatedCubeParticle = ps->create();
+               rotatedCubeParticle->setShapeID(cubeShapeID);
+               rotatedCubeParticle->setPosition(pos);
+               rotatedCubeParticle->setRotation(cubeRotation);
+               rotatedCubeParticle->setOwner(mpi::MPIManager::instance()->rank());
+               rotatedCubeParticle->setLinearVelocity(pos);
+               WALBERLA_LOG_INFO("Created cube particle");
+            }
+         }
+      }
+   }
+
+   /// VTK Output
+
+   auto vtkDomainOutput = walberla::vtk::createVTKOutput_DomainDecomposition(forest, "domain_decomposition", 1, vtk_out, "simulation_step");
+   vtkDomainOutput->write();
+
+   mesa_pd::MeshParticleVTKOutput< mesh::PolyMesh > meshParticleVTK(ps, ss, "mesh", uint_t(1));
+   meshParticleVTK.addFaceOutput< data::SelectParticleUid >("UID");
+   meshParticleVTK.addVertexOutput< data::SelectParticleInteractionRadius >("InteractionRadius");
+   meshParticleVTK.addFaceOutput< data::SelectParticleLinearVelocity >("LinearVelocity");
+   meshParticleVTK.addVertexOutput< data::SelectParticlePosition >("Position");
+   auto surfaceVelDataSource = make_shared<
+      mesa_pd::SurfaceVelocityVertexDataSource< mesh::PolyMesh, mesa_pd::data::ParticleAccessorWithShape > >(
+      "SurfaceVelocity", *ac);
+   meshParticleVTK.addVertexDataSource(surfaceVelDataSource);
+   meshParticleVTK();
+
+   return EXIT_SUCCESS;
+}
+}
+
+int main( int argc, char* argv[] )
+{
+   return walberla::main( argc, argv );
+}
\ No newline at end of file
diff --git a/tests/mesh/CMakeLists.txt b/tests/mesh/CMakeLists.txt
index a8ba46bf91600553536f74698437ee66c1869eb0..11fee5c5468060d6b09744de5a5acc219b652b91 100644
--- a/tests/mesh/CMakeLists.txt
+++ b/tests/mesh/CMakeLists.txt
@@ -96,5 +96,4 @@ if ( WALBERLA_BUILD_WITH_OPENMESH )
 
    waLBerla_compile_test( FILES MeshMarshalling.cpp DEPENDS core mesh )
    waLBerla_execute_test( NAME  MeshMarshalling )
-
 endif()