diff --git a/apps/benchmarks/CMakeLists.txt b/apps/benchmarks/CMakeLists.txt
index a81420b475e6378a3f453ebd4ae16c2c0c895e8d..31e6417ae85f60488e8938511b4db0a2f2533af2 100644
--- a/apps/benchmarks/CMakeLists.txt
+++ b/apps/benchmarks/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory( ComplexGeometry )
+add_subdirectory( DEM )
add_subdirectory( MeshDistance )
add_subdirectory( CouetteFlow )
add_subdirectory( ForcesOnSphereNearPlaneInShearFlow )
diff --git a/apps/benchmarks/DEM/CMakeLists.txt b/apps/benchmarks/DEM/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..138e92026d74c900b1f25cd172cef82de1b2bff2
--- /dev/null
+++ b/apps/benchmarks/DEM/CMakeLists.txt
@@ -0,0 +1 @@
+waLBerla_add_executable( NAME DEM FILES DEM.cpp DEPENDS blockforest core pe )
diff --git a/apps/benchmarks/DEM/DEM.cpp b/apps/benchmarks/DEM/DEM.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..52b6e856a4fb47403d3fb0e9db448c3b1756e460
--- /dev/null
+++ b/apps/benchmarks/DEM/DEM.cpp
@@ -0,0 +1,153 @@
+//======================================================================================================================
+//
+// 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 DEM.cpp
+//! \brief demonstration of basic functionality of DEM
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "pe/basic.h"
+
+#include <blockforest/Initialization.h>
+#include <core/DataTypes.h>
+
+#include <string>
+
+namespace walberla {
+
+namespace dem {
+real_t calcCoefficientOfRestitution(const real_t k, const real_t gamma, const real_t meff)
+{
+ auto a = real_t(0.5) * gamma / meff;
+ return std::exp(-a * math::PI / std::sqrt(k / meff - a*a));
+}
+
+real_t calcCollisionTime(const real_t k, const real_t gamma, const real_t meff)
+{
+ auto a = real_t(0.5) * gamma / meff;
+ return math::PI / std::sqrt( k/meff - a*a);
+}
+}
+
+int main( int argc, char** argv )
+{
+ using namespace walberla;
+ using namespace walberla::pe;
+
+ typedef boost::tuple<Sphere, Plane> BodyTuple ;
+
+ walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+
+ real_t dt = real_c(0.0001); //!< integration time
+ real_t radius = real_c(1); //!< particle radius
+ real_t density = real_c(2707); //!< particle density
+ real_t m = Sphere::calcMass( radius, density); //!< particle mass
+ std::string model = "kg"; //!< input model
+ real_t k = real_c(8.11e6); //!< linear spring stiffness
+ real_t gamma = real_c(6.86e1); //!< damper
+ real_t e = dem::calcCoefficientOfRestitution(k, gamma, m); //!< coefficient of restitution
+ real_t t = dem::calcCollisionTime(k, gamma, m); //!< collision time
+
+ for( int i = 1; i < argc; ++i )
+ {
+ if( std::strcmp( argv[i], "-dt" ) == 0 ) dt = real_c( std::stod( argv[++i] ) );
+ else if( std::strcmp( argv[i], "-radius" ) == 0 ) radius = real_c( std::stod( argv[++i] ) );
+ else if( std::strcmp( argv[i], "-density" ) == 0 ) density = real_c( std::stod( argv[++i] ) );
+ else if( std::strcmp( argv[i], "-kg" ) == 0 )
+ {
+ k = real_c( std::stod( argv[++i] ) );
+ gamma = real_c( std::stod( argv[++i] ) );
+ model = "kg";
+ }
+ else if( std::strcmp( argv[i], "-et" ) == 0 )
+ {
+ e = real_c( std::atof( argv[++i] ) );
+ t = real_c( std::atof( argv[++i] ) );
+ model = "et";
+ }
+ else WALBERLA_ABORT("Found invalid command line argument: \"" << argv[i] << "\" - aborting...");
+ }
+
+ m = Sphere::calcMass( radius, density);
+ if (model=="kg")
+ {
+ e = dem::calcCoefficientOfRestitution(k, gamma, m); //!< coefficient of restitution
+ t = dem::calcCollisionTime(k, gamma, m); //!< collision time
+ } else
+ WALBERLA_ABORT("unsupported model");
+
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+
+ // create blocks
+ shared_ptr< StructuredBlockForest > forest = blockforest::createUniformBlockGrid(
+ math::AABB(-5,-5,0,5,5,10),
+ uint_c( 1), uint_c( 1), uint_c( 1), // number of blocks in x,y,z direction
+ uint_c( 1), uint_c( 1), uint_c( 1), // how many cells per block (x,y,z)
+ true, // max blocks per process
+ false, false, false, // full periodicity
+ false);
+
+ SetBodyTypeIDs<BodyTuple>::execute();
+
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+ auto hccdID = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "HCCD");
+ auto fcdID = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
+ cr::DEM cr(globalBodyStorage, forest->getBlockStoragePointer(), storageID, hccdID, fcdID);
+
+ const real_t static_cof ( real_t(0.4) / real_t(2) ); // Coefficient of static friction. Roughly 0.85 with high variation depending on surface roughness for low stresses. Note: pe doubles the input coefficient of friction for material-material contacts.
+ const real_t dynamic_cof ( static_cof ); // Coefficient of dynamic friction. Similar to static friction for low speed friction.
+ MaterialID material = createMaterial( "granular", density, e, static_cof, dynamic_cof, real_t( 0.5 ), 1, k, gamma, 0 );
+
+ pe::createPlane( *globalBodyStorage, 0, Vec3(0, 0, 1), forest->getDomain().minCorner(), material );
+
+ SphereID sp = pe::createSphere(
+ *globalBodyStorage,
+ forest->getBlockStorage(),
+ storageID,
+ 999999999,
+ Vec3(0,0,1),
+ real_c(1.0),
+ material);
+ WALBERLA_CHECK_NOT_NULLPTR(sp);
+ sp->setLinearVel( Vec3(0,0,-1) );
+
+ uint_t steps = 0;
+ do
+ {
+ cr.timestep( dt );
+ ++steps;
+ } while (cr.getNumberOfContacts() != 0);
+
+ WALBERLA_LOG_RESULT(std::setw(30) << "steps: " << steps );
+ WALBERLA_LOG_RESULT(std::setw(30) << "final velocity: " << sp->getLinearVel()[2]);
+ WALBERLA_LOG_RESULT(std::setw(30) << "final position: " << sp->getPosition()[2]);
+ WALBERLA_LOG_RESULT(std::setw(30) << "integration time: " << dt);
+ WALBERLA_LOG_RESULT(std::setw(30) << "particle radius: " << radius);
+ WALBERLA_LOG_RESULT(std::setw(30) << "particle density: " << density);
+ WALBERLA_LOG_RESULT(std::setw(30) << "particle mass: " << m);
+ WALBERLA_LOG_RESULT(std::setw(30) << "linear spring stiffness: " << k);
+ WALBERLA_LOG_RESULT(std::setw(30) << "damper: " << gamma);
+ WALBERLA_LOG_RESULT(std::setw(30) << "coefficient of restitution: " << e);
+ WALBERLA_LOG_RESULT(std::setw(30) << "collision time: " << t);
+
+ return EXIT_SUCCESS;
+}
+}
+
+int main( int argc, char** argv )
+{
+ return walberla::main(argc, argv);
+}