diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index 80b5bff0e6b6355b795ed1bb024f55cb75c80d08..599c0e075f0469d5783ed9d499e16cd026477152 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -98,6 +98,7 @@ if __name__ == '__main__':
mpd.add(kernel.NonLinearSpringDashpot())
mpd.add(kernel.SingleCast(ps))
mpd.add(kernel.SpringDashpot())
+ mpd.add(kernel.SpringDashpotSpring())
mpd.add(kernel.TemperatureIntegration())
mpd.add(kernel.VelocityVerlet())
mpd.add(kernel.VelocityVerletWithShape())
diff --git a/python/mesa_pd/kernel/SpringDashpotSpring.py b/python/mesa_pd/kernel/SpringDashpotSpring.py
new file mode 100644
index 0000000000000000000000000000000000000000..9399e156332b7541dfc6c445a860920c3ba42238
--- /dev/null
+++ b/python/mesa_pd/kernel/SpringDashpotSpring.py
@@ -0,0 +1,11 @@
+# -*- coding: utf-8 -*-
+
+from mesa_pd.utility import generate_file
+
+
+class SpringDashpotSpring:
+ def generate(self, module):
+ ctx = {'module': module}
+ ctx["parameters"] = ["stiffnessN", "dampingN", "stiffnessT", "coefficientOfFriction"]
+
+ generate_file(module['module_path'], 'kernel/SpringDashpotSpring.templ.h', ctx)
diff --git a/python/mesa_pd/kernel/__init__.py b/python/mesa_pd/kernel/__init__.py
index eef9907a547ca08e93d35a26eed6d89ef005580b..8e79c416daafd913a50262539f7f0fbe36a335c7 100644
--- a/python/mesa_pd/kernel/__init__.py
+++ b/python/mesa_pd/kernel/__init__.py
@@ -15,6 +15,7 @@ from .LinearSpringDashpot import LinearSpringDashpot
from .NonLinearSpringDashpot import NonLinearSpringDashpot
from .SingleCast import SingleCast
from .SpringDashpot import SpringDashpot
+from .SpringDashpotSpring import SpringDashpotSpring
from .TemperatureIntegration import TemperatureIntegration
from .VelocityVerlet import VelocityVerlet
from .VelocityVerletWithShape import VelocityVerletWithShape
@@ -35,6 +36,7 @@ __all__ = ['DoubleCast',
'NonLinearSpringDashpot',
'SingleCast',
'SpringDashpot',
+ 'SpringDashpotSpring',
'TemperatureIntegration',
'VelocityVerlet',
'VelocityVerletWithShape']
diff --git a/python/mesa_pd/templates/kernel/SpringDashpotSpring.templ.h b/python/mesa_pd/templates/kernel/SpringDashpotSpring.templ.h
new file mode 100644
index 0000000000000000000000000000000000000000..78d3d4863175141fc76c657828c9c448804671bd
--- /dev/null
+++ b/python/mesa_pd/templates/kernel/SpringDashpotSpring.templ.h
@@ -0,0 +1,204 @@
+//======================================================================================================================
+//
+// 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 SpringDashpotSpring.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/common/ParticleFunctions.h>
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+
+#include <core/math/Constants.h>
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+
+/**
+ * Basic DEM kernel
+ *
+ * DEM kernel with spring-dashpot interaction in normal direction and spring in tangential direction.
+ *
+ * \ingroup mesa_pd_kernel
+ */
+class SpringDashpotSpring
+{
+public:
+ SpringDashpotSpring(const uint_t numParticleTypes);
+ SpringDashpotSpring(const SpringDashpotSpring& other) = default;
+ SpringDashpotSpring(SpringDashpotSpring&& other) = default;
+ SpringDashpotSpring& operator=(const SpringDashpotSpring& other) = default;
+ SpringDashpotSpring& operator=(SpringDashpotSpring&& other) = default;
+
+ template <typename Accessor>
+ void operator()(const size_t p_idx1,
+ const size_t p_idx2,
+ Accessor& ac,
+ const Vec3& contactPoint,
+ const Vec3& contactNormal,
+ const real_t& penetrationDepth,
+ const real_t dt) const;
+
+ {% for param in parameters %}
+ /// assumes this parameter is symmetric
+ void set{{param | capFirst}}(const size_t type1, const size_t type2, const real_t& val);
+ {%- endfor %}
+
+ {% for param in parameters %}
+ real_t get{{param | capFirst}}(const size_t type1, const size_t type2) const;
+ {%- endfor %}
+
+ inline
+ real_t calcCoefficientOfRestitution(const size_t type1,
+ const size_t type2,
+ const real_t meff)
+ {
+ auto a = real_t(0.5) * getDampingN(type1, type2) / meff;
+ return std::exp(-a * math::pi / std::sqrt(getStiffnessN(type1, type2) / meff - a*a));
+ }
+
+ inline
+ real_t calcCollisionTime(const size_t type1,
+ const size_t type2,
+ const real_t meff)
+ {
+ auto a = real_t(0.5) * getDampingN(type1, type2) / meff;
+ return math::pi / std::sqrt( getStiffnessN(type1, type2)/meff - a*a);
+ }
+
+ inline
+ void setParametersFromCOR(const size_t type1,
+ const size_t type2,
+ const real_t cor,
+ const real_t collisionTime,
+ const real_t meff)
+ {
+ const real_t lnDryResCoeff = std::log(cor);
+ setStiffnessN(type1, type2, math::pi * math::pi * meff / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::pi * math::pi + lnDryResCoeff* lnDryResCoeff )) ));
+ setDampingN( type1, type2, - real_t(2) * std::sqrt( meff * getStiffnessN(type1, type2) ) * ( lnDryResCoeff / std::sqrt( math::pi * math::pi + ( lnDryResCoeff * lnDryResCoeff ) ) ));
+ }
+private:
+ uint_t numParticleTypes_;
+ {% for param in parameters %}
+ std::vector<real_t> {{param}}_ {};
+ {%- endfor %}
+};
+
+SpringDashpotSpring::SpringDashpotSpring(const uint_t numParticleTypes)
+{
+ numParticleTypes_ = numParticleTypes;
+ {% for param in parameters %}
+ {{param}}_.resize(numParticleTypes * numParticleTypes, real_t(0));
+ {%- endfor %}
+}
+
+{% for param in parameters %}
+inline void SpringDashpotSpring::set{{param | capFirst}}(const size_t type1, const size_t type2, const real_t& val)
+{
+ WALBERLA_ASSERT_LESS( type1, numParticleTypes_ );
+ WALBERLA_ASSERT_LESS( type2, numParticleTypes_ );
+ {{param}}_[numParticleTypes_*type1 + type2] = val;
+ {{param}}_[numParticleTypes_*type2 + type1] = val;
+}
+{%- endfor %}
+
+{% for param in parameters %}
+inline real_t SpringDashpotSpring::get{{param | capFirst}}(const size_t type1, const size_t type2) const
+{
+ WALBERLA_ASSERT_LESS( type1, numParticleTypes_ );
+ WALBERLA_ASSERT_LESS( type2, numParticleTypes_ );
+ WALBERLA_ASSERT_FLOAT_EQUAL( {{param}}_[numParticleTypes_*type1 + type2],
+ {{param}}_[numParticleTypes_*type2 + type1],
+ "parameter matrix for {{param}} not symmetric!");
+ return {{param}}_[numParticleTypes_*type1 + type2];
+}
+{%- endfor %}
+
+template <typename Accessor>
+inline void SpringDashpotSpring::operator()(const size_t p_idx1,
+ const size_t p_idx2,
+ Accessor& ac,
+ const Vec3& contactPoint,
+ const Vec3& contactNormal,
+ const real_t& penetrationDepth,
+ const real_t dt) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+ if (p_idx1 != p_idx2)
+ {
+ // skip if no penetration is present
+ real_t delta = -penetrationDepth;
+ if (delta < real_t(0)) return;
+
+ // calculate relative velocities
+ const Vec3 relVel ( -(getVelocityAtWFPoint(p_idx1, ac, contactPoint) - getVelocityAtWFPoint(p_idx2, ac, contactPoint)) );
+ const real_t relVelN( math::dot(relVel, contactNormal) );
+ const Vec3 relVelT( relVel - ( relVelN * contactNormal ) );
+ const Vec3 contactTangent = relVelT.getNormalizedOrZero();
+
+ // Calculating the normal force based on a linear spring-dashpot force model
+ real_t fNabs = getStiffnessN(ac.getType(p_idx1), ac.getType(p_idx2)) * delta + getDampingN(ac.getType(p_idx1), ac.getType(p_idx2)) * relVelN;
+ const Vec3 fN = fNabs * contactNormal;
+
+ // get tangential displacement from contact history
+ auto tangentialDisplacement = Vec3(real_t(0));
+ auto contactHistory = ac.getOldContactHistoryRef(p_idx1).find(ac.getUid(p_idx2));
+ if(contactHistory != ac.getOldContactHistoryRef(p_idx1).end())
+ {
+ // get infos from the contact history
+ tangentialDisplacement = dot(contactHistory->second.getTangentialSpringDisplacement(), contactTangent) * contactTangent;
+ }
+
+ // accumulate tangential displacement
+ tangentialDisplacement += relVelT * dt;
+
+ // Calculating the tangential force
+ const auto maxTangentialForce = fNabs * getCoefficientOfFriction(ac.getType(p_idx1), ac.getType(p_idx2));
+ Vec3 fT = getStiffnessT(ac.getType(p_idx1), ac.getType(p_idx2)) * tangentialDisplacement;
+ if (length(fT) > maxTangentialForce)
+ fT = maxTangentialForce * fT.getNormalizedOrZero();
+
+ // store new tangential displacements
+ auto& ch1 = ac.getNewContactHistoryRef(p_idx1)[ac.getUid(p_idx2)];
+ ch1.setTangentialSpringDisplacement(tangentialDisplacement);
+
+ auto& ch2 = ac.getNewContactHistoryRef(p_idx2)[ac.getUid(p_idx1)];
+ ch2.setTangentialSpringDisplacement(tangentialDisplacement);
+
+ // Add normal force at contact point
+ addForceAtWFPosAtomic( p_idx1, ac, fN, contactPoint );
+ addForceAtWFPosAtomic( p_idx2, ac, -fN, contactPoint );
+
+ // Add tangential force at contact point
+ addForceAtWFPosAtomic( p_idx1, ac, fT, contactPoint );
+ addForceAtWFPosAtomic( p_idx2, ac, -fT, contactPoint );
+ }
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file
diff --git a/src/mesa_pd/kernel/SpringDashpotSpring.h b/src/mesa_pd/kernel/SpringDashpotSpring.h
new file mode 100644
index 0000000000000000000000000000000000000000..f2f87671a4e9bf08c87f17100f3b47a482f385e3
--- /dev/null
+++ b/src/mesa_pd/kernel/SpringDashpotSpring.h
@@ -0,0 +1,261 @@
+//======================================================================================================================
+//
+// 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 SpringDashpotSpring.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/common/ParticleFunctions.h>
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+
+#include <core/math/Constants.h>
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+
+/**
+ * Basic DEM kernel
+ *
+ * DEM kernel with spring-dashpot interaction in normal direction and spring in tangential direction.
+ *
+ * \ingroup mesa_pd_kernel
+ */
+class SpringDashpotSpring
+{
+public:
+ SpringDashpotSpring(const uint_t numParticleTypes);
+ SpringDashpotSpring(const SpringDashpotSpring& other) = default;
+ SpringDashpotSpring(SpringDashpotSpring&& other) = default;
+ SpringDashpotSpring& operator=(const SpringDashpotSpring& other) = default;
+ SpringDashpotSpring& operator=(SpringDashpotSpring&& other) = default;
+
+ template <typename Accessor>
+ void operator()(const size_t p_idx1,
+ const size_t p_idx2,
+ Accessor& ac,
+ const Vec3& contactPoint,
+ const Vec3& contactNormal,
+ const real_t& penetrationDepth,
+ const real_t dt) const;
+
+
+ /// assumes this parameter is symmetric
+ void setStiffnessN(const size_t type1, const size_t type2, const real_t& val);
+ /// assumes this parameter is symmetric
+ void setDampingN(const size_t type1, const size_t type2, const real_t& val);
+ /// assumes this parameter is symmetric
+ void setStiffnessT(const size_t type1, const size_t type2, const real_t& val);
+ /// assumes this parameter is symmetric
+ void setCoefficientOfFriction(const size_t type1, const size_t type2, const real_t& val);
+
+
+ real_t getStiffnessN(const size_t type1, const size_t type2) const;
+ real_t getDampingN(const size_t type1, const size_t type2) const;
+ real_t getStiffnessT(const size_t type1, const size_t type2) const;
+ real_t getCoefficientOfFriction(const size_t type1, const size_t type2) const;
+
+ inline
+ real_t calcCoefficientOfRestitution(const size_t type1,
+ const size_t type2,
+ const real_t meff)
+ {
+ auto a = real_t(0.5) * getDampingN(type1, type2) / meff;
+ return std::exp(-a * math::pi / std::sqrt(getStiffnessN(type1, type2) / meff - a*a));
+ }
+
+ inline
+ real_t calcCollisionTime(const size_t type1,
+ const size_t type2,
+ const real_t meff)
+ {
+ auto a = real_t(0.5) * getDampingN(type1, type2) / meff;
+ return math::pi / std::sqrt( getStiffnessN(type1, type2)/meff - a*a);
+ }
+
+ inline
+ void setParametersFromCOR(const size_t type1,
+ const size_t type2,
+ const real_t cor,
+ const real_t collisionTime,
+ const real_t meff)
+ {
+ const real_t lnDryResCoeff = std::log(cor);
+ setStiffnessN(type1, type2, math::pi * math::pi * meff / ( collisionTime * collisionTime * ( real_t(1) - lnDryResCoeff * lnDryResCoeff / ( math::pi * math::pi + lnDryResCoeff* lnDryResCoeff )) ));
+ setDampingN( type1, type2, - real_t(2) * std::sqrt( meff * getStiffnessN(type1, type2) ) * ( lnDryResCoeff / std::sqrt( math::pi * math::pi + ( lnDryResCoeff * lnDryResCoeff ) ) ));
+ }
+private:
+ uint_t numParticleTypes_;
+
+ std::vector<real_t> stiffnessN_ {};
+ std::vector<real_t> dampingN_ {};
+ std::vector<real_t> stiffnessT_ {};
+ std::vector<real_t> coefficientOfFriction_ {};
+};
+
+SpringDashpotSpring::SpringDashpotSpring(const uint_t numParticleTypes)
+{
+ numParticleTypes_ = numParticleTypes;
+
+ stiffnessN_.resize(numParticleTypes * numParticleTypes, real_t(0));
+ dampingN_.resize(numParticleTypes * numParticleTypes, real_t(0));
+ stiffnessT_.resize(numParticleTypes * numParticleTypes, real_t(0));
+ coefficientOfFriction_.resize(numParticleTypes * numParticleTypes, real_t(0));
+}
+
+
+inline void SpringDashpotSpring::setStiffnessN(const size_t type1, const size_t type2, const real_t& val)
+{
+ WALBERLA_ASSERT_LESS( type1, numParticleTypes_ );
+ WALBERLA_ASSERT_LESS( type2, numParticleTypes_ );
+ stiffnessN_[numParticleTypes_*type1 + type2] = val;
+ stiffnessN_[numParticleTypes_*type2 + type1] = val;
+}
+inline void SpringDashpotSpring::setDampingN(const size_t type1, const size_t type2, const real_t& val)
+{
+ WALBERLA_ASSERT_LESS( type1, numParticleTypes_ );
+ WALBERLA_ASSERT_LESS( type2, numParticleTypes_ );
+ dampingN_[numParticleTypes_*type1 + type2] = val;
+ dampingN_[numParticleTypes_*type2 + type1] = val;
+}
+inline void SpringDashpotSpring::setStiffnessT(const size_t type1, const size_t type2, const real_t& val)
+{
+ WALBERLA_ASSERT_LESS( type1, numParticleTypes_ );
+ WALBERLA_ASSERT_LESS( type2, numParticleTypes_ );
+ stiffnessT_[numParticleTypes_*type1 + type2] = val;
+ stiffnessT_[numParticleTypes_*type2 + type1] = val;
+}
+inline void SpringDashpotSpring::setCoefficientOfFriction(const size_t type1, const size_t type2, const real_t& val)
+{
+ WALBERLA_ASSERT_LESS( type1, numParticleTypes_ );
+ WALBERLA_ASSERT_LESS( type2, numParticleTypes_ );
+ coefficientOfFriction_[numParticleTypes_*type1 + type2] = val;
+ coefficientOfFriction_[numParticleTypes_*type2 + type1] = val;
+}
+
+
+inline real_t SpringDashpotSpring::getStiffnessN(const size_t type1, const size_t type2) const
+{
+ WALBERLA_ASSERT_LESS( type1, numParticleTypes_ );
+ WALBERLA_ASSERT_LESS( type2, numParticleTypes_ );
+ WALBERLA_ASSERT_FLOAT_EQUAL( stiffnessN_[numParticleTypes_*type1 + type2],
+ stiffnessN_[numParticleTypes_*type2 + type1],
+ "parameter matrix for stiffnessN not symmetric!");
+ return stiffnessN_[numParticleTypes_*type1 + type2];
+}
+inline real_t SpringDashpotSpring::getDampingN(const size_t type1, const size_t type2) const
+{
+ WALBERLA_ASSERT_LESS( type1, numParticleTypes_ );
+ WALBERLA_ASSERT_LESS( type2, numParticleTypes_ );
+ WALBERLA_ASSERT_FLOAT_EQUAL( dampingN_[numParticleTypes_*type1 + type2],
+ dampingN_[numParticleTypes_*type2 + type1],
+ "parameter matrix for dampingN not symmetric!");
+ return dampingN_[numParticleTypes_*type1 + type2];
+}
+inline real_t SpringDashpotSpring::getStiffnessT(const size_t type1, const size_t type2) const
+{
+ WALBERLA_ASSERT_LESS( type1, numParticleTypes_ );
+ WALBERLA_ASSERT_LESS( type2, numParticleTypes_ );
+ WALBERLA_ASSERT_FLOAT_EQUAL( stiffnessT_[numParticleTypes_*type1 + type2],
+ stiffnessT_[numParticleTypes_*type2 + type1],
+ "parameter matrix for stiffnessT not symmetric!");
+ return stiffnessT_[numParticleTypes_*type1 + type2];
+}
+inline real_t SpringDashpotSpring::getCoefficientOfFriction(const size_t type1, const size_t type2) const
+{
+ WALBERLA_ASSERT_LESS( type1, numParticleTypes_ );
+ WALBERLA_ASSERT_LESS( type2, numParticleTypes_ );
+ WALBERLA_ASSERT_FLOAT_EQUAL( coefficientOfFriction_[numParticleTypes_*type1 + type2],
+ coefficientOfFriction_[numParticleTypes_*type2 + type1],
+ "parameter matrix for coefficientOfFriction not symmetric!");
+ return coefficientOfFriction_[numParticleTypes_*type1 + type2];
+}
+
+template <typename Accessor>
+inline void SpringDashpotSpring::operator()(const size_t p_idx1,
+ const size_t p_idx2,
+ Accessor& ac,
+ const Vec3& contactPoint,
+ const Vec3& contactNormal,
+ const real_t& penetrationDepth,
+ const real_t dt) const
+{
+ static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+ if (p_idx1 != p_idx2)
+ {
+ // skip if no penetration is present
+ real_t delta = -penetrationDepth;
+ if (delta < real_t(0)) return;
+
+ // calculate relative velocities
+ const Vec3 relVel ( -(getVelocityAtWFPoint(p_idx1, ac, contactPoint) - getVelocityAtWFPoint(p_idx2, ac, contactPoint)) );
+ const real_t relVelN( math::dot(relVel, contactNormal) );
+ const Vec3 relVelT( relVel - ( relVelN * contactNormal ) );
+ const Vec3 contactTangent = relVelT.getNormalizedOrZero();
+
+ // Calculating the normal force based on a linear spring-dashpot force model
+ real_t fNabs = getStiffnessN(ac.getType(p_idx1), ac.getType(p_idx2)) * delta + getDampingN(ac.getType(p_idx1), ac.getType(p_idx2)) * relVelN;
+ const Vec3 fN = fNabs * contactNormal;
+
+ // get tangential displacement from contact history
+ auto tangentialDisplacement = Vec3(real_t(0));
+ auto contactHistory = ac.getOldContactHistoryRef(p_idx1).find(ac.getUid(p_idx2));
+ if(contactHistory != ac.getOldContactHistoryRef(p_idx1).end())
+ {
+ // get infos from the contact history
+ tangentialDisplacement = dot(contactHistory->second.getTangentialSpringDisplacement(), contactTangent) * contactTangent;
+ }
+
+ // accumulate tangential displacement
+ tangentialDisplacement += relVelT * dt;
+
+ // Calculating the tangential force
+ const auto maxTangentialForce = fNabs * getCoefficientOfFriction(ac.getType(p_idx1), ac.getType(p_idx2));
+ Vec3 fT = getStiffnessT(ac.getType(p_idx1), ac.getType(p_idx2)) * tangentialDisplacement;
+ if (length(fT) > maxTangentialForce)
+ fT = maxTangentialForce * fT.getNormalizedOrZero();
+
+ // store new tangential displacements
+ auto& ch1 = ac.getNewContactHistoryRef(p_idx1)[ac.getUid(p_idx2)];
+ ch1.setTangentialSpringDisplacement(tangentialDisplacement);
+
+ auto& ch2 = ac.getNewContactHistoryRef(p_idx2)[ac.getUid(p_idx1)];
+ ch2.setTangentialSpringDisplacement(tangentialDisplacement);
+
+ // Add normal force at contact point
+ addForceAtWFPosAtomic( p_idx1, ac, fN, contactPoint );
+ addForceAtWFPosAtomic( p_idx2, ac, -fN, contactPoint );
+
+ // Add tangential force at contact point
+ addForceAtWFPosAtomic( p_idx1, ac, fT, contactPoint );
+ addForceAtWFPosAtomic( p_idx2, ac, -fT, contactPoint );
+ }
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index 81a8a74cad16373f4b5f39856e7c25a4d13a0f2f..ffc72839000545a515160c4a50591265f49aa23b 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -134,6 +134,9 @@ waLBerla_execute_test( NAME MESA_PD_Kernel_LinkedCellsVsBruteForce PROCESSES 2
waLBerla_compile_test( NAME MESA_PD_Kernel_SingleCast FILES kernel/SingleCast.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_SingleCast )
+waLBerla_compile_test( NAME MESA_PD_Kernel_SpherePile FILES kernel/SpherePile.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_Kernel_SpherePile )
+
waLBerla_compile_test( NAME MESA_PD_Kernel_SpringDashpot FILES kernel/SpringDashpot.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_SpringDashpot )
diff --git a/tests/mesa_pd/kernel/SpherePile.cpp b/tests/mesa_pd/kernel/SpherePile.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e326a1110ad8f0ef4d0e32d2a275313f11c02879
--- /dev/null
+++ b/tests/mesa_pd/kernel/SpherePile.cpp
@@ -0,0 +1,211 @@
+//======================================================================================================================
+//
+// 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 SpherePile.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
+#include "mesa_pd/common/ParticleFunctions.h"
+
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+
+#include "mesa_pd/kernel/DoubleCast.h"
+#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/SpringDashpot.h"
+#include "mesa_pd/kernel/SpringDashpotSpring.h"
+#include "mesa_pd/mpi/ReduceContactHistory.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+#include "core/Environment.h"
+#include "core/logging/Logging.h"
+#include "vtk/VTKOutput.h"
+
+#include <iostream>
+
+namespace walberla {
+namespace mesa_pd {
+
+auto createSphere(data::ParticleStorage &ps, const Vec3 &pos)
+{
+ auto p = ps.create();
+ p->setPosition(pos);
+ p->setType(0);
+ return p;
+}
+
+/*
+ * Simulates oblique sphere-wall collision and checks rebound angle, i.e. the tangential part of the collision model.
+ *
+ */
+int main(int argc, char **argv)
+{
+ walberla::mpi::Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+ walberla::mpi::MPIManager::instance()->useWorldComm();
+
+ auto dt = real_t(0.001);
+ auto cor = real_t(0.2);
+ auto ct = dt * real_t(20);
+ auto radius = real_t(1);
+ auto density = real_t(2700);
+ auto simSteps = 100000;
+
+ //init data structures
+ auto ps = walberla::make_shared<data::ParticleStorage>(2);
+ auto ss = walberla::make_shared<data::ShapeStorage>();
+ auto ac = walberla::make_shared<data::ParticleAccessorWithShape>(ps, ss);
+
+ auto sphereShape = ss->create<data::Sphere>(radius);
+ ss->shapes[sphereShape]->updateMassAndInertia(density);
+
+ // create sphere
+ createSphere(*ps, Vec3(0, 0, 0));
+ createSphere(*ps, Vec3(2, 0, 0));
+ auto sp = createSphere(*ps, Vec3(1, 0, std::sqrt(real_t(4) - real_t(1))));
+
+ // create plane
+ data::Particle &&p0 = *ps->create(true);
+ p0.setPosition(Vec3(0, 0, -1));
+ p0.setShapeID(ss->create<data::HalfSpace>(Vector3<real_t>(0, 0, 1)));
+ p0.setType(0);
+ data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::INFINITE);
+ data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::FIXED);
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** VTK ***");
+ auto vtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps) ;
+ auto vtkWriter = walberla::vtk::createVTKOutput_PointData(vtkOutput,
+ "Bodies",
+ 1,
+ "vtk",
+ "simulation_step",
+ false,
+ false);
+
+ // explicit euler
+ kernel::ExplicitEulerWithShape explEuler(dt);
+ collision_detection::AnalyticContactDetection acd;
+ kernel::DoubleCast double_cast;
+ kernel::SpringDashpot sd(1);
+ kernel::SpringDashpotSpring sds(1);
+ mpi::ReduceContactHistory rch;
+
+ WALBERLA_LOG_DEVEL_VAR(sp->getPosition());
+
+ for (auto i = 0; i < simSteps; ++i)
+ {
+ ps->forEachParticle(false,
+ kernel::SelectLocal(),
+ *ac,
+ [&](const size_t idx1)
+ {
+ if (ac->getShape(idx1)->getShapeType() == data::Sphere::SHAPE_TYPE)
+ {
+ ac->setForce(idx1, Vec3(0,0,-1) * real_t(9.81) / ac->getInvMass(idx1));
+ }
+ });
+ ps->forEachParticlePairHalf(false,
+ kernel::SelectAll(),
+ *ac,
+ [&](const size_t idx1, const size_t idx2)
+ {
+ if ((ac->getShape(idx1)->getShapeType() != data::Sphere::SHAPE_TYPE) &&
+ (ac->getShape(idx2)->getShapeType() != data::Sphere::SHAPE_TYPE))
+ {
+ //skip plane - plane collision
+ } else
+ {
+ if (double_cast(idx1, idx2, *ac, acd, *ac))
+ {
+ //WALBERLA_LOG_DEVEL_VAR(acd.getPenetrationDepth());
+ auto meff = real_t(1) / (ac->getInvMass(idx1) + ac->getInvMass(idx2));
+ sds.setParametersFromCOR(0, 0, cor, ct, meff);
+ sds.setCoefficientOfFriction(0,0, real_t(0.4));
+ sds.setStiffnessT(0,0, sds.getStiffnessN(0,0));
+ sds(acd.getIdx1(), acd.getIdx2(), *ac, acd.getContactPoint(),
+ acd.getContactNormal(), acd.getPenetrationDepth(), dt);
+ }
+ }
+ });
+
+ rch(*ps);
+ ps->forEachParticle(false,
+ kernel::SelectLocal(),
+ *ac,
+ explEuler,
+ *ac);
+ }
+ WALBERLA_LOG_DEVEL_VAR(sp->getPosition());
+ WALBERLA_CHECK_GREATER(sp->getPosition()[2], real_t(1));
+
+ for (auto i = 0; i < simSteps; ++i)
+ {
+ ps->forEachParticle(false,
+ kernel::SelectLocal(),
+ *ac,
+ [&](const size_t idx1)
+ {
+ if (ac->getShape(idx1)->getShapeType() == data::Sphere::SHAPE_TYPE)
+ {
+ ac->setForce(idx1, Vec3(0,0,-1) * real_t(9.81) / ac->getInvMass(idx1));
+ }
+ });
+ ps->forEachParticlePairHalf(false,
+ kernel::SelectAll(),
+ *ac,
+ [&](const size_t idx1, const size_t idx2)
+ {
+ if ((ac->getShape(idx1)->getShapeType() != data::Sphere::SHAPE_TYPE) &&
+ (ac->getShape(idx2)->getShapeType() != data::Sphere::SHAPE_TYPE))
+ {
+ //skip plane - plane collision
+ } else
+ {
+ if (double_cast(idx1, idx2, *ac, acd, *ac))
+ {
+ //WALBERLA_LOG_DEVEL_VAR(acd.getPenetrationDepth());
+ auto meff = real_t(1) / (ac->getInvMass(idx1) + ac->getInvMass(idx2));
+ sd.setParametersFromCOR(0, 0, cor, ct, meff);
+ sd.setFriction(0,0, real_t(0.4));
+ sd.setDampingT(0,0, sd.getDampingN(0,0));
+ sd(acd.getIdx1(), acd.getIdx2(), *ac, acd.getContactPoint(),
+ acd.getContactNormal(), acd.getPenetrationDepth());
+ }
+ }
+ });
+
+ rch(*ps);
+ ps->forEachParticle(false,
+ kernel::SelectLocal(),
+ *ac,
+ explEuler,
+ *ac);
+ }
+ WALBERLA_LOG_DEVEL_VAR(sp->getPosition());
+ WALBERLA_CHECK_LESS(sp->getPosition()[2], real_t(1));
+
+ return EXIT_SUCCESS;
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
+int main(int argc, char **argv)
+{
+ return walberla::mesa_pd::main(argc, argv);
+}