Merge branch 'mr_mesapd_shape_integration_fix' into 'master'

Fix of integration of angular velocity of non-spherical particles

See merge request walberla/walberla!527

apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp

@@ -18,7 +18,6 @@
 //
 //======================================================================================================================
 
-#include "Accessor.h"
 #include "check.h"
 #include "Contact.h"
 #include "CreateParticles.h"
@@ -32,7 +31,7 @@
 
 #include <mesa_pd/collision_detection/AnalyticContactDetection.h>
 #include <mesa_pd/data/LinkedCells.h>
-#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleAccessorWithShape.h>
 #include <mesa_pd/data/ParticleStorage.h>
 #include <mesa_pd/data/ShapeStorage.h>
 #include <mesa_pd/domain/BlockForestDomain.h>
@@ -90,7 +89,7 @@ public:
    }
 
    inline
-   void operator()(const size_t idx1, const size_t idx2, ParticleAccessorWithShape& ac)
+   void operator()(const size_t idx1, const size_t idx2, data::ParticleAccessorWithShape& ac)
    {
 
       ++contactsChecked_;
@@ -194,7 +193,7 @@ int main( int argc, char ** argv )
    //init data structures
    auto ps = std::make_shared<data::ParticleStorage>(100);
    auto ss = std::make_shared<data::ShapeStorage>();
-   ParticleAccessorWithShape accessor(ps, ss);
+   data::ParticleAccessorWithShape accessor(ps, ss);
    data::LinkedCells     lc(localDomain.getExtended(params.spacing), params.spacing );
 
    auto  smallSphere = ss->create<data::Sphere>( params.radius );

apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp

@@ -18,7 +18,6 @@
 //
 //======================================================================================================================
 
-#include "Accessor.h"
 #include "check.h"
 #include "Contact.h"
 #include "CreateParticles.h"
@@ -32,7 +31,7 @@
 
 #include <mesa_pd/collision_detection/AnalyticContactDetection.h>
 #include <mesa_pd/data/LinkedCells.h>
-#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleAccessorWithShape.h>
 #include <mesa_pd/data/ParticleStorage.h>
 #include <mesa_pd/data/ShapeStorage.h>
 #include <mesa_pd/domain/BlockForestDomain.h>
@@ -209,7 +208,7 @@ int main( int argc, char ** argv )
    //init data structures
    auto ss = std::make_shared<data::ShapeStorage>();
    auto ps = std::make_shared<data::ParticleStorage>(100);
-   ParticleAccessorWithShape accessor(ps, ss);
+   data::ParticleAccessorWithShape accessor(ps, ss);
    data::LinkedCells     lc(localDomain.getExtended(params.spacing), params.spacing );
 
    auto  smallSphere = ss->create<data::Sphere>( params.radius );

apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp

@@ -18,7 +18,6 @@
 //
 //======================================================================================================================
 
-#include "Accessor.h"
 #include "check.h"
 #include "Contact.h"
 #include "CreateParticles.h"
@@ -32,7 +31,7 @@
 
 #include <mesa_pd/collision_detection/AnalyticContactDetection.h>
 #include <mesa_pd/data/LinkedCells.h>
-#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleAccessorWithShape.h>
 #include <mesa_pd/data/ParticleStorage.h>
 #include <mesa_pd/data/ShapeStorage.h>
 #include <mesa_pd/data/SparseLinkedCells.h>
@@ -198,7 +197,7 @@ int main( int argc, char ** argv )
    //init data structures
    auto ps = std::make_shared<data::ParticleStorage>(100);
    auto ss = std::make_shared<data::ShapeStorage>();
-   ParticleAccessorWithShape accessor(ps, ss);
+   data::ParticleAccessorWithShape accessor(ps, ss);
    auto lc = std::make_shared<data::LinkedCells>(domain->getUnionOfLocalAABBs().getExtended(params.spacing), params.spacing );
    forest->addBlockData(domain::createBlockForestDataHandling(ps), "Storage");
 

apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp

@@ -18,7 +18,6 @@
 //
 //======================================================================================================================
 
-#include "Accessor.h"
 #include "check.h"
 #include "Contact.h"
 #include "CreateParticles.h"
@@ -32,7 +31,7 @@
 
 #include <mesa_pd/collision_detection/AnalyticContactDetection.h>
 #include <mesa_pd/data/SparseLinkedCells.h>
-#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleAccessorWithShape.h>
 #include <mesa_pd/data/ParticleStorage.h>
 #include <mesa_pd/data/ShapeStorage.h>
 #include <mesa_pd/data/STLOverloads.h>
@@ -236,7 +235,7 @@ int main( int argc, char ** argv )
    //init data structures
    auto ps = std::make_shared<data::ParticleStorage>(100);
    auto ss = std::make_shared<data::ShapeStorage>();
-   ParticleAccessorWithShape accessor(ps, ss);
+   data::ParticleAccessorWithShape accessor(ps, ss);
    auto lc = std::make_shared<data::SparseLinkedCells>(domain->getUnionOfLocalAABBs().getExtended(params.spacing), params.spacing );
    forest->addBlockData(domain::createBlockForestDataHandling(ps), "Storage");
 

apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.cpp

@@ -43,15 +43,15 @@ public:
    const auto &getInvMass(const size_t /*p_idx*/) const
    { return invMass_; }
 
-   void setInvInertiaBF(const size_t /*p_idx*/, const Mat3 &val)
-   { invInertiaBF_ = val; }
+   const auto &getInvInertiaBF(const size_t /*p_idx*/) const // dummy
+   { return dummyI_; }
 
-   const auto &getInvInertiaBF(const size_t /*p_idx*/) const
-   { return invInertiaBF_; }
+   const auto &getInertiaBF(const size_t /*p_idx*/) const // dummy
+   { return dummyI_; }
 
 private:
    real_t invMass_;
-   Mat3 invInertiaBF_;
+   Mat3 dummyI_{0_r};
 };
 
 struct Oscillator

apps/showcases/ParticlePacking/Evaluation.h

@@ -20,6 +20,7 @@
 
 #pragma once
 
+#include "mesa_pd/data/ContactAccessor.h"
 #include "mesa_pd/data/ParticleStorage.h"
 
 #include "Utility.h"

apps/showcases/ParticlePacking/ParticlePacking.cpp

@@ -35,6 +35,7 @@
 
 #include "mesa_pd/data/ContactStorage.h"
 #include "mesa_pd/data/ContactAccessor.h"
+#include "mesa_pd/data/ParticleAccessorWithBaseShape.h"
 #include "mesa_pd/data/ParticleStorage.h"
 #include "mesa_pd/data/HashGrids.h"
 
@@ -418,7 +419,7 @@ int main(int argc, char **argv) {
    /// MESAPD Data
    auto particleStorage = std::make_shared<data::ParticleStorage>(1);
    auto contactStorage = std::make_shared<data::ContactStorage>(1);
-   ParticleAccessorWithShape particleAccessor(particleStorage);
+   data::ParticleAccessorWithBaseShape particleAccessor(particleStorage);
    data::ContactAccessor contactAccessor(contactStorage);
 
    // configure shape creation
@@ -626,7 +627,7 @@ int main(int argc, char **argv) {
 
    math::DistributedSample diameterSample;
    particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
-                                    [&diameterSample](const size_t idx, ParticleAccessorWithShape& ac){diameterSample.insert(real_c(2)*ac.getInteractionRadius(idx));}, particleAccessor);
+                                    [&diameterSample](const size_t idx, data::ParticleAccessorWithBaseShape& ac){diameterSample.insert(real_c(2)*ac.getInteractionRadius(idx));}, particleAccessor);
 
    diameterSample.mpiAllGather();
    WALBERLA_LOG_INFO_ON_ROOT("Statistics of initially created particles' interaction diameters: " << diameterSample.format());
@@ -713,7 +714,7 @@ int main(int argc, char **argv) {
    meshParticleVTK.addFaceOutput< data::SelectParticleLinearVelocity >("LinearVelocity");
    meshParticleVTK.addVertexOutput< data::SelectParticlePosition >("Position");
    meshParticleVTK.addVertexOutput< data::SelectParticleNumContacts >("numContacts");
-   auto surfaceVelDataSource = make_shared<mesa_pd::SurfaceVelocityVertexDataSource< mesh::PolyMesh, ParticleAccessorWithShape > >("SurfaceVelocity", particleAccessor);
+   auto surfaceVelDataSource = make_shared<mesa_pd::SurfaceVelocityVertexDataSource< mesh::PolyMesh, data::ParticleAccessorWithBaseShape > >("SurfaceVelocity", particleAccessor);
    meshParticleVTK.setParticleSelector(vtkParticleSelector);
    meshParticleVTK.addVertexDataSource(surfaceVelDataSource);
 
@@ -850,7 +851,7 @@ int main(int argc, char **argv) {
 
          timing.start("Contact detection");
          hashGrids.forEachParticlePairHalf(useOpenMP, kernel::ExcludeInfiniteInfinite(), particleAccessor,
-                                           [domain, contactStorage](size_t idx1, size_t idx2, ParticleAccessorWithShape &ac){
+                                           [domain, contactStorage](size_t idx1, size_t idx2, data::ParticleAccessorWithBaseShape &ac){
                                               kernel::DoubleCast double_cast;
                                               mpi::ContactFilter contact_filter;
                                               collision_detection::GeneralContactDetection contactDetection;
@@ -890,7 +891,7 @@ int main(int argc, char **argv) {
          } else
          {
             linkedCells.forEachParticlePairHalf(useOpenMP, kernel::ExcludeInfiniteInfinite(), particleAccessor,
-                                                [domain, contactStorage](size_t idx1, size_t idx2, ParticleAccessorWithShape &ac){
+                                                [domain, contactStorage](size_t idx1, size_t idx2, data::ParticleAccessorWithBaseShape &ac){
 
                                                    collision_detection::GeneralContactDetection contactDetection;
                                                    //Attention: does not use contact threshold in general case (GJK)
@@ -926,9 +927,9 @@ int main(int argc, char **argv) {
 
       timing.start("Contact eval");
       particleStorage->forEachParticle(useOpenMP, kernel::SelectAll(), particleAccessor,
-                                       [](size_t p_idx, ParticleAccessorWithShape& ac){ac.setNumContacts(p_idx,0);}, particleAccessor);
+                                       [](size_t p_idx, data::ParticleAccessorWithBaseShape& ac){ac.setNumContacts(p_idx,0);}, particleAccessor);
       contactStorage->forEachContact(useOpenMP, kernel::SelectAll(), contactAccessor,
-                                     [](size_t c, data::ContactAccessor &ca, ParticleAccessorWithShape &pa) {
+                                     [](size_t c, data::ContactAccessor &ca, data::ParticleAccessorWithBaseShape &pa) {
                                         auto idx1 = ca.getId1(c);
                                         auto idx2 = ca.getId2(c);
                                         pa.getNumContactsRef(idx1)++;
@@ -996,7 +997,7 @@ int main(int argc, char **argv) {
          timing.start("DEM");
          timing.start("Collision");
          contactStorage->forEachContact(useOpenMP, kernel::SelectAll(), contactAccessor,
-                                        [&dem_collision, coefficientOfRestitution, dem_collisionTime, dem_kappa, dt](size_t c, data::ContactAccessor &ca, ParticleAccessorWithShape &pa){
+                                        [&dem_collision, coefficientOfRestitution, dem_collisionTime, dem_kappa, dt](size_t c, data::ContactAccessor &ca, data::ParticleAccessorWithBaseShape &pa){
                                            auto idx1 = ca.getId1(c);
                                            auto idx2 = ca.getId2(c);
                                            auto meff = real_t(1) / (pa.getInvMass(idx1) + pa.getInvMass(idx2));
@@ -1125,7 +1126,7 @@ int main(int argc, char **argv) {
 
          // apply damping
          particleStorage->forEachParticle(useOpenMP, kernel::SelectAll(), particleAccessor,
-                                          [velocityDampingFactor](size_t idx, ParticleAccessorWithShape &ac){
+                                          [velocityDampingFactor](size_t idx, data::ParticleAccessorWithBaseShape &ac){
                                              ac.getLinearVelocityRef(idx) *= velocityDampingFactor;
                                              ac.getAngularVelocityRef(idx) *= velocityDampingFactor;},
                                            particleAccessor);

apps/showcases/ParticlePacking/Utility.h

@@ -20,7 +20,11 @@
 
 #pragma once
 
+#include "core/mpi/MPITextFile.h"
+
 #include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/shape/HalfSpace.h"
+#include "mesa_pd/data/shape/CylindricalBoundary.h"
 
 #include <iterator>
 #include <algorithm>
@@ -29,22 +33,6 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& particleStorage)
-         : ParticleAccessor(particleStorage)
-   {}
-
-   walberla::real_t const & getInvMass(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInvMass();}
-
-   Mat3 const & getInvInertiaBF(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInvInertiaBF();}
-   Mat3 const & getInertiaBF(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInertiaBF();}
-   data::BaseShape* getShape(const size_t p_idx) const {return ps_->getBaseShape(p_idx).get();}
-};
-
-
-
 template< typename T>
 std::vector<T> parseStringToVector(std::string inputStr)
 {

python/mesa_pd/kernel/ExplicitEuler.py

@@ -20,6 +20,7 @@ class ExplicitEuler:
             self.context['interface'].append(create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"))
             self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"))
             self.context['interface'].append(create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"))
+            self.context['interface'].append(create_access("inertiaBF", "walberla::mesa_pd::Mat3", access="g"))
             self.context['interface'].append(create_access("torque", "walberla::mesa_pd::Vec3", access="gs"))
 
     def generate(self, module):

python/mesa_pd/kernel/SemiImplicitEuler.py

@@ -20,6 +20,7 @@ class SemiImplicitEuler:
             self.context['interface'].append(create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"))
             self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"))
             self.context['interface'].append(create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"))
+            self.context['interface'].append(create_access("inertiaBF", "walberla::mesa_pd::Mat3", access="g"))
             self.context['interface'].append(create_access("torque", "walberla::mesa_pd::Vec3", access="gs"))
 
     def generate(self, module):

python/mesa_pd/kernel/VelocityVerlet.py

@@ -17,6 +17,7 @@ class VelocityVerlet:
             self.context['interface'].append(create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"))
             self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"))
             self.context['interface'].append(create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"))
+            self.context['interface'].append(create_access("inertiaBF", "walberla::mesa_pd::Mat3", access="g"))
             self.context['interface'].append(create_access("torque", "walberla::mesa_pd::Vec3", access="gs"))
             self.context['interface'].append(create_access("oldTorque", "walberla::mesa_pd::Vec3", access="gs"))
 

python/mesa_pd/templates/common/ParticleFunctions.templ.h

@@ -67,6 +67,21 @@ inline Vec3 transformVectorFromBFtoWF(const size_t p_idx, Accessor& ac, const Ve
    return ac.getRotation(p_idx).getMatrix() * vectorBF;
 }
 
+/**
+ * Transform (inverse) particle's moment of inertia from body frame coordinates (as stored by shape) to world frame.
+ */
+template <typename Accessor>
+inline Mat3 getInvInertia(const size_t p_idx, Accessor& ac)
+{
+   return math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(), ac.getInvInertiaBF(p_idx));
+}
+
+template <typename Accessor>
+inline Mat3 getInertia(const size_t p_idx, Accessor& ac)
+{
+   return math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(), ac.getInertiaBF(p_idx));
+}
+
 /**
  * Force is applied at the center of mass.
  */

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

@@ -28,6 +28,9 @@
 
 #include <mesa_pd/data/DataTypes.h>
 #include <mesa_pd/data/IAccessor.h>
+{%- if bIntegrateRotation %}
+#include <mesa_pd/common/ParticleFunctions.h>
+{%- endif %}
 
 namespace walberla {
 namespace mesa_pd {
@@ -87,8 +90,12 @@ inline void ExplicitEuler::operator()(const size_t idx,
                                 ac.getLinearVelocity(idx));
 
       {%- if bIntegrateRotation %}
-      const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
-                                                  ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
+      // computation done in body frame: d(omega)/ dt = J^-1 ((J*omega) x omega + T), update in world frame
+      // see Wachs, 2019, doi:10.1007/s00707-019-02389-9, Eq. 27
+      const auto omegaBF = transformVectorFromWFtoBF(idx, ac, ac.getAngularVelocity(idx));
+      const auto torqueBF = transformVectorFromWFtoBF(idx, ac, ac.getTorque(idx));
+      const Vec3 wdotBF = ac.getInvInertiaBF(idx) * ( ( ac.getInertiaBF(idx) * omegaBF ) % omegaBF + torqueBF );
+      const Vec3 wdot = transformVectorFromBFtoWF(idx, ac, wdotBF);
 
       // Calculating the rotation angle
       const Vec3 phi( 0.5_r * wdot * dt_ * dt_ +

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

@@ -87,7 +87,8 @@ inline void InitParticlesForHCSITS::operator()(size_t j, Accessor& ac, real_t dt
       initializeVelocityCorrections( ac, j, ac.getDvRef(j), ac.getDwRef(j), dt ); // use applied external forces to calculate starting velocity
       if(!isSet(particle_flags, FIXED)){ // Update velocities with global acceleration and angular velocity with euler eqn
          ac.getLinearVelocityRef(j) = ac.getLinearVelocity(j) + getGlobalAcceleration() * dt;
-         ac.getAngularVelocityRef(j) = ac.getAngularVelocity(j) + dt * ( ac.getInvInertiaBF(j) * ( ( ac.getInertiaBF(j) * ac.getAngularVelocity(j) ) % ac.getAngularVelocity(j) ) );
+         const auto omegaBF = transformVectorFromWFtoBF(j, ac, ac.getAngularVelocity(j));
+         ac.getAngularVelocityRef(j) = ac.getAngularVelocity(j) + dt * transformVectorFromBFtoWF(j, ac, ( ac.getInvInertiaBF(j) * ( ( ac.getInertiaBF(j) * omegaBF ) % omegaBF ) ) );
       }
    }
 }
@@ -110,7 +111,7 @@ template <typename Accessor>
 inline void InitParticlesForHCSITS::initializeVelocityCorrections(Accessor& ac, size_t body, Vec3& dv, Vec3& dw, real_t dt ) const
 {
    dv = ( ac.getInvMass(body) * dt ) * ac.getForce(body);
-   dw = dt * ( ac.getInvInertiaBF(body) * ac.getTorque(body) );
+   dw = dt * ( getInvInertia(body, ac) * ac.getTorque(body) );
 
    ac.getForceRef(body) = Vec3();
    ac.getTorqueRef(body) = Vec3();

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

@@ -28,6 +28,9 @@
 
 #include <mesa_pd/data/DataTypes.h>
 #include <mesa_pd/data/IAccessor.h>
+{%- if bIntegrateRotation %}
+#include <mesa_pd/common/ParticleFunctions.h>
+{%- endif %}
 
 namespace walberla {
 namespace mesa_pd {
@@ -81,8 +84,12 @@ inline void SemiImplicitEuler::operator()(const size_t idx,
                                  ac.getPosition(idx));
 
       {%- if bIntegrateRotation %}
-      const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
-                                                  ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
+      // computation done in body frame: d(omega)/ dt = J^-1 ((J*omega) x omega + T), update in world frame
+      // see Wachs, 2019, doi:10.1007/s00707-019-02389-9, Eq. 27
+      const auto omegaBF = transformVectorFromWFtoBF(idx, ac, ac.getAngularVelocity(idx));
+      const auto torqueBF = transformVectorFromWFtoBF(idx, ac, ac.getTorque(idx));
+      const Vec3 wdotBF = ac.getInvInertiaBF(idx) * ( ( ac.getInertiaBF(idx) * omegaBF ) % omegaBF + torqueBF );
+      const Vec3 wdot = transformVectorFromBFtoWF(idx, ac, wdotBF);
 
       ac.setAngularVelocity(idx, wdot * dt_ +
                                  ac.getAngularVelocity(idx));

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

@@ -28,6 +28,9 @@
 
 #include <mesa_pd/data/DataTypes.h>
 #include <mesa_pd/data/IAccessor.h>
+{%- if bIntegrateRotation %}
+#include <mesa_pd/common/ParticleFunctions.h>
+{%- endif %}
 
 namespace walberla {
 namespace mesa_pd {
@@ -95,8 +98,14 @@ inline void VelocityVerletPreForceUpdate::operator()(const size_t p_idx, Accesso
                             real_t(0.5) * ac.getInvMass(p_idx) * ac.getOldForce(p_idx) * dt_ * dt_);
 
       {%- if bIntegrateRotation %}
-      const Vec3 wdot = math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(),
-                                                  ac.getInvInertiaBF(p_idx)) * ac.getOldTorque(p_idx);
+      // computation done in body frame: d(omega)/ dt = J^-1 ((J*omega) x omega + T), update in world frame
+      // see Wachs, 2019, doi:10.1007/s00707-019-02389-9, Eq. 27
+      // note: this implementation (pre and post) is experimental as it is in principle unclear in which order
+      //       angular velocities and rotations (affect again the transformations WF - BF) have to be carried out
+      const auto omegaBF = transformVectorFromWFtoBF(p_idx, ac, ac.getAngularVelocity(p_idx));
+      const auto torqueBF = transformVectorFromWFtoBF(p_idx, ac, ac.getOldTorque(p_idx));
+      const Vec3 wdotBF = ac.getInvInertiaBF(p_idx) * ( ( ac.getInertiaBF(p_idx) * omegaBF ) % omegaBF + torqueBF );
+      const Vec3 wdot = transformVectorFromBFtoWF(p_idx, ac, wdotBF);
 
       // Calculating the rotation angle
       const Vec3 phi( ac.getAngularVelocity(p_idx) * dt_ + real_t(0.5) * wdot * dt_ * dt_);
@@ -120,12 +129,13 @@ inline void VelocityVerletPostForceUpdate::operator()(const size_t p_idx, Access
                                   real_t(0.5) * ac.getInvMass(p_idx) * (ac.getOldForce(p_idx) + ac.getForce(p_idx)) * dt_);
 
       {%- if bIntegrateRotation %}
-      const auto torque = ac.getOldTorque(p_idx) + ac.getTorque(p_idx);
-      const Vec3 wdot = math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(),
-                                                  ac.getInvInertiaBF(p_idx)) * torque;
+      const auto omegaBF = transformVectorFromWFtoBF(p_idx, ac, ac.getAngularVelocity(p_idx));
+      const auto torqueBF = transformVectorFromWFtoBF(p_idx, ac, 0.5_r * (ac.getOldTorque(p_idx) + ac.getTorque(p_idx)));
+      const Vec3 wdotBF = ac.getInvInertiaBF(p_idx) * ( ( ac.getInertiaBF(p_idx) * omegaBF ) % omegaBF + torqueBF );
+      const Vec3 wdot = transformVectorFromBFtoWF(p_idx, ac, wdotBF);
 
       ac.setAngularVelocity(p_idx, ac.getAngularVelocity(p_idx) +
-                                   real_t(0.5) * wdot * dt_ );
+                                   wdot * dt_ );
       {%- endif %}
    }
 

src/mesa_pd/common/ParticleFunctions.h

@@ -67,6 +67,21 @@ inline Vec3 transformVectorFromBFtoWF(const size_t p_idx, Accessor& ac, const Ve
    return ac.getRotation(p_idx).getMatrix() * vectorBF;
 }
 
+/**
+ * Transform (inverse) particle's moment of inertia from body frame coordinates (as stored by shape) to world frame.
+ */
+template <typename Accessor>
+inline Mat3 getInvInertia(const size_t p_idx, Accessor& ac)
+{
+   return math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(), ac.getInvInertiaBF(p_idx));
+}
+
+template <typename Accessor>
+inline Mat3 getInertia(const size_t p_idx, Accessor& ac)
+{
+   return math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(), ac.getInertiaBF(p_idx));
+}
+
 /**
  * Force is applied at the center of mass.
  */

apps/benchmarks/GranularGas/Accessor.h → src/mesa_pd/data/ParticleAccessorWithBaseShape.h

@@ -13,33 +13,38 @@
 //  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   Accessor.h
-//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//! \file
+//! \author Christoph Rettinger <sebastian.eibl@fau.de>
 //
 //======================================================================================================================
 
+#pragma once
+
 #include <mesa_pd/data/ParticleAccessor.h>
 #include <mesa_pd/data/ParticleStorage.h>
-#include <mesa_pd/data/ShapeStorage.h>
 
 namespace walberla {
 namespace mesa_pd {
+namespace data {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
+class ParticleAccessorWithBaseShape : public data::ParticleAccessor
 {
 public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
+   ParticleAccessorWithBaseShape(std::shared_ptr<data::ParticleStorage>& ps)
       : ParticleAccessor(ps)
-      , ss_(ss)
    {}
 
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
+   const auto& getInvMass(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInvMass();}
+   const auto& getMass(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getMass();}
+
+   const auto& getInvInertiaBF(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInvInertiaBF();}
+   const auto& getInertiaBF(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInertiaBF();}
+
+   auto getVolume(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getVolume();}
 
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
+   data::BaseShape* getShape(const size_t p_idx) const {return ps_->getBaseShape(p_idx).get();}
 };
 
-} // namespace mesa_pd
-} // namespace walberla
+} //namespace data
+} //namespace mesa_pd
+} //namespace walberla

src/mesa_pd/data/ParticleAccessorWithShape.h

@@ -40,7 +40,7 @@ public:
    const auto& getMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getMass();}
    const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
    const auto& getInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInertiaBF();}
-   const auto getVolume(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getVolume();}
+   auto getVolume(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getVolume();}
 
    data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)].get();}
 private:

src/mesa_pd/kernel/ExplicitEuler.h

@@ -28,6 +28,7 @@
 
 #include <mesa_pd/data/DataTypes.h>
 #include <mesa_pd/data/IAccessor.h>
+#include <mesa_pd/common/ParticleFunctions.h>
 
 namespace walberla {
 namespace mesa_pd {
@@ -64,6 +65,8 @@ namespace kernel {
  *
  * const walberla::mesa_pd::Mat3& getInvInertiaBF(const size_t p_idx) const;
  *
+ * const walberla::mesa_pd::Mat3& getInertiaBF(const size_t p_idx) const;
+ *
  * const walberla::mesa_pd::Vec3& getTorque(const size_t p_idx) const;
  * void setTorque(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
  *
@@ -97,8 +100,12 @@ inline void ExplicitEuler::operator()(const size_t idx,
                                 ac.getPosition(idx));
       ac.setLinearVelocity(idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ +
                                 ac.getLinearVelocity(idx));
-      const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
-                                                  ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
+      // computation done in body frame: d(omega)/ dt = J^-1 ((J*omega) x omega + T), update in world frame
+      // see Wachs, 2019, doi:10.1007/s00707-019-02389-9, Eq. 27
+      const auto omegaBF = transformVectorFromWFtoBF(idx, ac, ac.getAngularVelocity(idx));
+      const auto torqueBF = transformVectorFromWFtoBF(idx, ac, ac.getTorque(idx));
+      const Vec3 wdotBF = ac.getInvInertiaBF(idx) * ( ( ac.getInertiaBF(idx) * omegaBF ) % omegaBF + torqueBF );
+      const Vec3 wdot = transformVectorFromBFtoWF(idx, ac, wdotBF);
 
       // Calculating the rotation angle
       const Vec3 phi( 0.5_r * wdot * dt_ * dt_ +