diff --git a/src/mesh/pe/raytracing/Intersects.h b/src/mesh/pe/raytracing/Intersects.h
new file mode 100644
index 0000000000000000000000000000000000000000..52761d65ef6be955d2d505ef6e4e114e2b9fec54
--- /dev/null
+++ b/src/mesh/pe/raytracing/Intersects.h
@@ -0,0 +1,83 @@
+//======================================================================================================================
+//
+// 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 Intersects.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesh/MatrixVectorOperations.h>
+#include <mesh/pe/rigid_body/ConvexPolyhedron.h>
+#include <mesh/TriangleMeshes.h>
+#include <pe/raytracing/Intersects.h>
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+// Implemented following the description in
+// "FAST RAY - CONVEX POLYHEDRON INTERSECTION" by Eric Haines
+inline bool intersects(const mesh::pe::ConvexPolyhedronID poly, const Ray& ray, real_t& t_near, Vec3& n)
+{
+ Ray transformedRay = ray.transformedToBF(poly);
+
+ t_near = std::numeric_limits<real_t>::min();
+ real_t t_far = std::numeric_limits<real_t>::max();
+ const mesh::TriangleMesh& mesh = poly->getMesh();
+
+ for(auto fh : mesh.faces())
+ {
+ //plane equation x*Pn + d = 0
+ const Vector3<real_t> Pn = mesh::toWalberla(mesh.normal(fh)); // Plane normal
+ const real_t d = - Pn * mesh::toWalberla(mesh.point(mesh.to_vertex_handle(mesh.halfedge_handle(fh))));
+ const real_t vn = Pn * transformedRay.getOrigin() + d;
+ const real_t vd = Pn * transformedRay.getDirection();
+
+ if ( floatIsEqual(vd, real_t(0)) )
+ {
+ if (vn > real_t(0)) return false;
+ continue;
+ }
+
+ const real_t t = -vn / vd;
+
+ if (vd > real_t(0))
+ {
+ // back-facing
+ if (t < real_t(0)) return false;
+ if (t < t_far) t_far=t;
+ } else
+ {
+ // front-facing
+ if (t > t_near)
+ {
+ t_near = t;
+ n = Pn;
+ }
+ if (t_near > t_far) return false;
+ }
+ }
+
+ if (t_near < 0) return false;
+
+ n = poly->vectorFromBFtoWF(n);
+ return true;
+}
+
+} //namespace raytracing
+} //namespace pe
+} //namespace walberla
diff --git a/tests/mesh/CMakeLists.txt b/tests/mesh/CMakeLists.txt
index 3eebf6ed6ba80f3caeed8d6597b0be18d417358a..932e96ec634405d4aea068b24bb2129bc2c1a06b 100644
--- a/tests/mesh/CMakeLists.txt
+++ b/tests/mesh/CMakeLists.txt
@@ -53,6 +53,9 @@ if ( WALBERLA_BUILD_WITH_OPENMESH )
waLBerla_compile_test( FILES PeVTKMeshWriterTest.cpp DEPENDS mesh )
waLBerla_execute_test( NAME PeVTKMeshWriterTest COMMAND $<TARGET_FILE:PeVTKMeshWriterTest> )
+ waLBerla_compile_test( FILES MeshPeRaytracing.cpp DEPENDS mesh pe )
+ waLBerla_execute_test( NAME MeshPeRaytracing )
+
waLBerla_compile_test( FILES MeshAABBIntersectionTest.cpp DEPENDS mesh )
waLBerla_execute_test( NAME MeshAABBIntersectionTestDbg COMMAND $<TARGET_FILE:MeshAABBIntersectionTest> 10 )
waLBerla_execute_test( NAME MeshAABBIntersectionTest COMMAND $<TARGET_FILE:MeshAABBIntersectionTest> 10000 CONFIGURATIONS Release RelWithDbgInfo )
diff --git a/tests/mesh/MeshPeRaytracing.cpp b/tests/mesh/MeshPeRaytracing.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0cde05aa066621ae749391fa372a9757bd89f7b6
--- /dev/null
+++ b/tests/mesh/MeshPeRaytracing.cpp
@@ -0,0 +1,129 @@
+//======================================================================================================================
+//
+// 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 MeshPeRaytracing.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <core/debug/TestSubsystem.h>
+#include <core/logging/Logging.h>
+#include <core/mpi/Environment.h>
+
+#include <core/DataTypes.h>
+#include <core/math/Vector3.h>
+#include <mesh/pe/raytracing/Intersects.h>
+#include <mesh/pe/rigid_body/ConvexPolyhedron.h>
+#include <mesh/QHull.h>
+#include <pe/Materials.h>
+#include <pe/raytracing/Intersects.h>
+#include <pe/rigidbody/Box.h>
+
+#include <vector>
+
+namespace walberla{
+namespace mesh {
+namespace pe {
+
+int CpRayIntersectionTest(const int resolution = 10)
+{
+ using namespace walberla::math;
+ using namespace walberla::pe::raytracing;
+
+ std::vector<Vector3<real_t>> points;
+ points.push_back( Vector3<real_t>( real_t(-1), real_t(-1), real_t(-1) ) );
+ points.push_back( Vector3<real_t>( real_t(-1), real_t(-1), real_t( 1) ) );
+ points.push_back( Vector3<real_t>( real_t(-1), real_t( 1), real_t(-1) ) );
+ points.push_back( Vector3<real_t>( real_t(-1), real_t( 1), real_t( 1) ) );
+ points.push_back( Vector3<real_t>( real_t( 1), real_t(-1), real_t(-1) ) );
+ points.push_back( Vector3<real_t>( real_t( 1), real_t(-1), real_t( 1) ) );
+ points.push_back( Vector3<real_t>( real_t( 1), real_t( 1), real_t(-1) ) );
+ points.push_back( Vector3<real_t>( real_t( 1), real_t( 1), real_t( 1) ) );
+
+ shared_ptr< TriangleMesh > mesh = make_shared<TriangleMesh>();
+ mesh::QHull<TriangleMesh> qhull( points, mesh );
+ qhull.run();
+
+ const Vec3 center(1,2,3);
+
+ ConvexPolyhedron cp(0, 0, center, Vec3(0,0,0), Quat(), *mesh, Material::find("iron"), false, true, true);
+ cp.rotate(real_t(1), real_t(2), real_t(3));
+ Box bx(0, 0, center, Vec3(0,0,0), Quat(), Vec3(2,2,2), Material::find("iron"), false, true, true);
+ bx.rotate(real_t(1), real_t(2), real_t(3));
+
+ real_t dx = real_t(1.0) / static_cast<real_t>(resolution);
+ //rays pointed at center of body
+ for (int x = 0; x < resolution; ++x)
+ {
+ WALBERLA_LOG_INFO("[" << x+1 << " / " << resolution << "]" );
+ const real_t rand1 = real_c(x) * dx;
+ for (int y = 0; y < resolution; ++y)
+ {
+ const real_t rand2 = real_c(y) * dx;
+ real_t theta = 2 * M_PI * rand1;
+ real_t phi = acos(1 - 2 * rand2);
+ Vec3 dir(sin(phi) * cos(theta), sin(phi) * sin(theta), cos(phi));
+
+ Ray ray( center + dir*real_t(5), -dir);
+ real_t bx_t, cp_t;
+ Vec3 bx_n, cp_n;
+ WALBERLA_CHECK( intersects(&bx, ray, bx_t, bx_n) );
+ WALBERLA_CHECK( intersects(&cp, ray, cp_t, cp_n) );
+ WALBERLA_CHECK_FLOAT_EQUAL(bx_t, cp_t);
+ WALBERLA_CHECK_FLOAT_EQUAL(bx_n, cp_n);
+ }
+ }
+
+ //rays emitted form a point outside of the body
+ for (int x = 0; x < resolution; ++x)
+ {
+ WALBERLA_LOG_INFO("[" << x+1 << " / " << resolution << "]" );
+ const real_t rand1 = real_c(x) * dx;
+ for (int y = 0; y < resolution; ++y)
+ {
+ const real_t rand2 = real_c(y) * dx;
+ real_t theta = 2 * M_PI * rand1;
+ real_t phi = acos(1 - 2 * rand2);
+ Vec3 dir(sin(phi) * cos(theta), sin(phi) * sin(theta), cos(phi));
+
+ Ray ray( Vec3(5,5,5), -dir);
+ real_t bx_t, cp_t;
+ Vec3 bx_n, cp_n;
+ const bool bx_intersects = intersects(&bx, ray, bx_t, bx_n);
+ const bool cp_intersects = intersects(&cp, ray, cp_t, cp_n);
+ WALBERLA_CHECK_EQUAL( bx_intersects, cp_intersects );
+ if (bx_intersects)
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL(bx_t, cp_t);
+ WALBERLA_CHECK_FLOAT_EQUAL(bx_n, cp_n);
+ }
+ }
+ }
+
+ return EXIT_SUCCESS;
+}
+
+} //namespace pe
+} //namespace mesh
+} //namespace walberla
+
+int main( int argc, char * argv[] )
+{
+ walberla::debug::enterTestMode();
+ walberla::mpi::Environment mpiEnv( argc, argv );
+ walberla::mpi::MPIManager::instance()->useWorldComm();
+
+ return walberla::mesh::pe::CpRayIntersectionTest(10);
+}