introduced new spring-dashpot spring kernel

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

python/mesa_pd/kernel/SpringDashpotSpring.py

+# -*- 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)

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']

python/mesa_pd/templates/kernel/SpringDashpotSpring.templ.h

+//======================================================================================================================
+//
+//  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

src/mesa_pd/kernel/SpringDashpotSpring.h

+//======================================================================================================================
+//
+//  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

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 )
 

tests/mesa_pd/kernel/SpherePile.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file   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);
+}