diff --git a/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp b/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp
index 4ddbe5d8b7bb80288207a95ebbb9a6b0c1dbd00e..fe7b9ff8f3e8556b67000cdc3181b7ec76341d7c 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp
+++ b/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 );
diff --git a/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp b/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp
index d15e0dfd58670b59fcf091da026f600dec977864..6a7ce88b44efad1cd7b88223230adb52f80e1779 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp
+++ b/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 );
diff --git a/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp b/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp
index e3e11a43dcf4f224d44def73093e5d115110a307..1763884931f5784d8a472997af75e5c8a9a0262b 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp
+++ b/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");
 
diff --git a/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp b/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp
index c327dbd26541f73dfe3b3ddbb8b2cf08e59cd713..b8360a14f655be9258dcbe47fb67d22a069d2875 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp
+++ b/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");
 
diff --git a/apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.cpp b/apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.cpp
index af2dc8a3aaa54f805822eabcbd4139fb49d3cdf7..7e8709ea80d4336bb2407e8549e6aa69fcdcb9f7 100644
--- a/apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.cpp
+++ b/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
diff --git a/apps/showcases/ParticlePacking/Evaluation.h b/apps/showcases/ParticlePacking/Evaluation.h
index 930fec6ac51ccbbc510237e58fc58b5b86e56d46..c1a531a6d1e9e0603158b3e91019e582b9d10bb0 100644
--- a/apps/showcases/ParticlePacking/Evaluation.h
+++ b/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"
diff --git a/apps/showcases/ParticlePacking/ParticlePacking.cpp b/apps/showcases/ParticlePacking/ParticlePacking.cpp
index 41f97f8b9822f9f6fa18499c64f6a0b2fd5d565e..2c55d064f867d6541fbb48bb93941a088f6a46fd 100644
--- a/apps/showcases/ParticlePacking/ParticlePacking.cpp
+++ b/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);
diff --git a/apps/showcases/ParticlePacking/Utility.h b/apps/showcases/ParticlePacking/Utility.h
index 5e0f0b4c60dac22fe89db64973b3caba4c8b072d..86b6a24c77491f32aef87d3dc88760ba0ca78cac 100644
--- a/apps/showcases/ParticlePacking/Utility.h
+++ b/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)
 {
diff --git a/python/mesa_pd/kernel/ExplicitEuler.py b/python/mesa_pd/kernel/ExplicitEuler.py
index 1b51adb853ddbcf7174c1b2c3a23c48c641ff619..6d7b39fe131db0aa7f96262d29e7aa04b77d6c12 100644
--- a/python/mesa_pd/kernel/ExplicitEuler.py
+++ b/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):
diff --git a/python/mesa_pd/kernel/SemiImplicitEuler.py b/python/mesa_pd/kernel/SemiImplicitEuler.py
index f5b6632a1f57d7dba5631896a3bd26b5a5cce177..6732a78cf3240d247150d444a332623ae6d209b7 100644
--- a/python/mesa_pd/kernel/SemiImplicitEuler.py
+++ b/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):
diff --git a/python/mesa_pd/kernel/VelocityVerlet.py b/python/mesa_pd/kernel/VelocityVerlet.py
index 1a31554cba85349fe7fba8fd755fbff62f540915..1bb4ff445a361df45af7d07c284f653f384b17c4 100644
--- a/python/mesa_pd/kernel/VelocityVerlet.py
+++ b/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"))
 
diff --git a/python/mesa_pd/templates/common/ParticleFunctions.templ.h b/python/mesa_pd/templates/common/ParticleFunctions.templ.h
index 599c8d58b45939d9ced44fee0326531f55191800..61629f05909a60069ccb2b91d69bf79c8a034f2f 100644
--- a/python/mesa_pd/templates/common/ParticleFunctions.templ.h
+++ b/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.
  */
diff --git a/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h b/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h
index d6893ffbe2ae8f6845963812fbdfda1af694a6bf..4cb2935173b8db0274381dcb8ce2e0251fa27741 100644
--- a/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h
+++ b/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_ +
diff --git a/python/mesa_pd/templates/kernel/InitParticlesForHCSITS.templ.h b/python/mesa_pd/templates/kernel/InitParticlesForHCSITS.templ.h
index 9769d399943bc8ff4141af9aa89449d735822458..48c167315501684f5e380839e544083344a0859d 100644
--- a/python/mesa_pd/templates/kernel/InitParticlesForHCSITS.templ.h
+++ b/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();
diff --git a/python/mesa_pd/templates/kernel/SemiImplicitEuler.templ.h b/python/mesa_pd/templates/kernel/SemiImplicitEuler.templ.h
index edc393536385bbd471d153215976eae2dc07affb..528139501e26b5a300b8f9dfe191333aed87c9b5 100644
--- a/python/mesa_pd/templates/kernel/SemiImplicitEuler.templ.h
+++ b/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));
diff --git a/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h b/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h
index 08797d155282f6ec6b0b7b0ded6818a79a46d7f3..9d37ed22e94ccedcf8a7763950be78507d4914ba 100644
--- a/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h
+++ b/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 %}
    }
 
diff --git a/src/mesa_pd/common/ParticleFunctions.h b/src/mesa_pd/common/ParticleFunctions.h
index 52b043d1289ae7fdd3778512f5b93f4e697a827c..155ca3af1a8b2fd9627f5005f6def49c8b98092d 100644
--- a/src/mesa_pd/common/ParticleFunctions.h
+++ b/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.
  */
diff --git a/apps/benchmarks/GranularGas/Accessor.h b/src/mesa_pd/data/ParticleAccessorWithBaseShape.h
similarity index 58%
rename from apps/benchmarks/GranularGas/Accessor.h
rename to src/mesa_pd/data/ParticleAccessorWithBaseShape.h
index 1d34ce30eb3530257c916014fb46edc915737a1f..e9de554b8d23db088b18c608fb7a76da5ff53243 100644
--- a/apps/benchmarks/GranularGas/Accessor.h
+++ b/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
diff --git a/src/mesa_pd/data/ParticleAccessorWithShape.h b/src/mesa_pd/data/ParticleAccessorWithShape.h
index 8c884d3704251d54382207d1aae96c28272ca802..9f956839ac700575b9b5cbcffa78f68d7e3b2b17 100644
--- a/src/mesa_pd/data/ParticleAccessorWithShape.h
+++ b/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:
diff --git a/src/mesa_pd/kernel/ExplicitEuler.h b/src/mesa_pd/kernel/ExplicitEuler.h
index bf1bfa6609731a66d48a05cd5ab9995c30b4994e..44afb7b9c898ec133d754ba6d60bb2820a9b9885 100644
--- a/src/mesa_pd/kernel/ExplicitEuler.h
+++ b/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_ +
diff --git a/src/mesa_pd/kernel/ExplicitEulerWithShape.h b/src/mesa_pd/kernel/ExplicitEulerWithShape.h
index 06ed0c23211f051dbc81fd481c749f454bc7f509..982511b8cdcaa8fd6cdff7262dd577da6341be8f 100644
--- a/src/mesa_pd/kernel/ExplicitEulerWithShape.h
+++ b/src/mesa_pd/kernel/ExplicitEulerWithShape.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 {
@@ -89,8 +90,13 @@ inline void ExplicitEuler::operator()(const size_t idx,
    {
       ac.setPosition      (idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ * dt_ + ac.getLinearVelocity(idx) * dt_ + 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( ac.getAngularVelocity(idx) * dt_ + wdot * dt_ * dt_);
diff --git a/src/mesa_pd/kernel/InitParticlesForHCSITS.h b/src/mesa_pd/kernel/InitParticlesForHCSITS.h
index 2022a0b804b4911442e3d8a7ba7f54966e8dad48..6f6f7b6e34b5af38740696bf5254f349f914fe27 100644
--- a/src/mesa_pd/kernel/InitParticlesForHCSITS.h
+++ b/src/mesa_pd/kernel/InitParticlesForHCSITS.h
@@ -83,7 +83,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 ) ) );
       }
    }
 }
@@ -106,7 +107,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();
diff --git a/src/mesa_pd/kernel/SemiImplicitEuler.h b/src/mesa_pd/kernel/SemiImplicitEuler.h
index 13963d4f59c22f0a03d865eceac18e8b5183f3e2..bf8cec461b1e2bea8fd722f0ac4b36c7f218b640 100644
--- a/src/mesa_pd/kernel/SemiImplicitEuler.h
+++ b/src/mesa_pd/kernel/SemiImplicitEuler.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 {
@@ -59,6 +60,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);
  *
@@ -91,8 +94,12 @@ inline void SemiImplicitEuler::operator()(const size_t idx,
                                  ac.getLinearVelocity(idx));
       ac.setPosition      ( idx, ac.getLinearVelocity(idx) * dt_ +
                                  ac.getPosition(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);
 
       ac.setAngularVelocity(idx, wdot * dt_ +
                                  ac.getAngularVelocity(idx));
diff --git a/src/mesa_pd/kernel/VelocityVerlet.h b/src/mesa_pd/kernel/VelocityVerlet.h
index 67a5f1dc2c58bcbc2a79b432768304b060fbab09..93d64239d337264f479ba0272c2e6d9b4835c0f5 100644
--- a/src/mesa_pd/kernel/VelocityVerlet.h
+++ b/src/mesa_pd/kernel/VelocityVerlet.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 {
@@ -67,6 +68,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);
  *
@@ -111,8 +114,14 @@ inline void VelocityVerletPreForceUpdate::operator()(const size_t p_idx, Accesso
       ac.setPosition(p_idx, ac.getPosition(p_idx) +
                             ac.getLinearVelocity(p_idx) * dt_ +
                             real_t(0.5) * ac.getInvMass(p_idx) * ac.getOldForce(p_idx) * dt_ * dt_);
-      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_);
@@ -133,12 +142,13 @@ inline void VelocityVerletPostForceUpdate::operator()(const size_t p_idx, Access
    {
       ac.setLinearVelocity(p_idx, ac.getLinearVelocity(p_idx) +
                                   real_t(0.5) * ac.getInvMass(p_idx) * (ac.getOldForce(p_idx) + ac.getForce(p_idx)) * dt_);
-      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_ );
    }
 
    ac.setOldForce(p_idx,       ac.getForce(p_idx));
diff --git a/tests/mesa_pd/ContactDetection.cpp b/tests/mesa_pd/ContactDetection.cpp
index 65196877218ce561dcd8fe17bd83dac8a27d9fd3..e01eeb38d293309782344f0d5cad3ed0164b6a6f 100644
--- a/tests/mesa_pd/ContactDetection.cpp
+++ b/tests/mesa_pd/ContactDetection.cpp
@@ -21,7 +21,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>
@@ -45,23 +45,6 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
 /*
  * Tests linked cells
  * Fully periodic, simple cubic sphere array is created
@@ -119,7 +102,7 @@ int main( const int particlesPerAxisPerProcess = 2 )
    //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);
    const real_t linkedCellSize = real_c(linkedCellMultipleOfGenerationSpacing) * generationSpacing;
    data::LinkedCells lc(localDomain.getExtended(linkedCellSize), linkedCellSize );
 
diff --git a/tests/mesa_pd/DropTestAnalytic.cpp b/tests/mesa_pd/DropTestAnalytic.cpp
index 827f0951f2b51294b3c3f486fd54a777164f885e..01bea0bea488874f35fe58734a708d172439e79e 100644
--- a/tests/mesa_pd/DropTestAnalytic.cpp
+++ b/tests/mesa_pd/DropTestAnalytic.cpp
@@ -20,7 +20,7 @@
 
 #include <mesa_pd/collision_detection/AnalyticContactDetection.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>
 
@@ -42,23 +42,6 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
 int main( int argc, char ** argv )
 {
    using namespace walberla::timing;
@@ -75,7 +58,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);
+   mesa_pd::data::ParticleAccessorWithShape accessor(ps, ss);
 
    auto p                       = ps->create();
    p->getPositionRef()          = Vec3(real_t(0), real_t(0), real_t(0));
diff --git a/tests/mesa_pd/DropTestGeneral.cpp b/tests/mesa_pd/DropTestGeneral.cpp
index 52495a885a843f77480fee07e918a746603fdb17..08f2bb82f64bde82048468d852e31ad9d6a220f1 100644
--- a/tests/mesa_pd/DropTestGeneral.cpp
+++ b/tests/mesa_pd/DropTestGeneral.cpp
@@ -20,7 +20,7 @@
 
 #include <mesa_pd/collision_detection/GeneralContactDetection.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>
 
@@ -40,23 +40,6 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
 int main( int argc, char ** argv )
 {
    using namespace walberla::timing;
@@ -73,7 +56,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);
+   mesa_pd::data::ParticleAccessorWithShape accessor(ps, ss);
 
    auto p                       = ps->create();
    p->getPositionRef()          = Vec3(real_t(0), real_t(0), real_t(0));
diff --git a/tests/mesa_pd/collision_detection/AnalyticContactDetection.cpp b/tests/mesa_pd/collision_detection/AnalyticContactDetection.cpp
index 911e7344df4608c7de743e27f05bc0427b7dbf35..3860c58017289379536dc0c3a08cebb56df874b6 100644
--- a/tests/mesa_pd/collision_detection/AnalyticContactDetection.cpp
+++ b/tests/mesa_pd/collision_detection/AnalyticContactDetection.cpp
@@ -19,7 +19,7 @@
 //======================================================================================================================
 
 #include <mesa_pd/collision_detection/AnalyticContactDetection.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/kernel/DoubleCast.h>
@@ -34,23 +34,6 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
 void checkRotation( const Vec3& from, const Vec3& to )
 {
    WALBERLA_LOG_DEVEL_VAR(from);
@@ -75,7 +58,7 @@ void checkSphereSphereCollision( )
    //init data structures
    auto ps = std::make_shared<data::ParticleStorage>(100);
    auto ss = std::make_shared<data::ShapeStorage>();
-   ParticleAccessorWithShape ac(ps, ss);
+   data::ParticleAccessorWithShape ac(ps, ss);
 
    //initialize particles
    const real_t radius  = real_t(0.5);
@@ -125,7 +108,7 @@ void checkSphereHalfSpaceCollision( )
    //init data structures
    auto ps = std::make_shared<data::ParticleStorage>(100);
    auto ss = std::make_shared<data::ShapeStorage>();
-   ParticleAccessorWithShape ac(ps, ss);
+   data::ParticleAccessorWithShape ac(ps, ss);
 
    //initialize particles
    const real_t radius  = real_t(1.0);
diff --git a/tests/mesa_pd/collision_detection/GeneralContactDetection.cpp b/tests/mesa_pd/collision_detection/GeneralContactDetection.cpp
index 1dcb49882b02fd7cae711158f4bbc1d4817b0e60..e56fdbb16c3bb5ae854e94c282814e7c979f775a 100644
--- a/tests/mesa_pd/collision_detection/GeneralContactDetection.cpp
+++ b/tests/mesa_pd/collision_detection/GeneralContactDetection.cpp
@@ -19,40 +19,25 @@
 //======================================================================================================================
 
 #include <mesa_pd/collision_detection/GeneralContactDetection.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/kernel/DoubleCast.h>
 
 #include <core/Abort.h>
 #include <core/Environment.h>
-#include <core/logging/Logging.h>
-#include <core/waLBerlaBuildInfo.h>
 
 #include <memory>
 
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-      : ParticleAccessor(ps)
-      , ss_(ss)
-   {}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
 void generalContactDetection()
 {
    //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 e0 = ps->create();
    e0->setPosition(Vec3(real_t(0),real_t(0),real_t(0)));
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
index b425ece1be55911c6733ec3f09442459e8cf6029..2ce861a46d29d01dfe02c406861696039726ee5e 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
@@ -20,7 +20,7 @@
 
 #include "mesa_pd/collision_detection/AnalyticContactDetection.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"
 
@@ -40,24 +40,6 @@ using namespace walberla;
 using namespace walberla::mesa_pd;
 
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
-
 /*
  * Tests the integrator accuracy for a DEM simulation by comparing the given coefficient of restitution to the simulated one.
  * For that, the velocity after a single sphere-wall collision is divided by the initial velocity before the simulation.
@@ -99,7 +81,7 @@ int main( int argc, char** argv )
    //init data structures
    auto ps = walberla::make_shared<data::ParticleStorage>(2);
    auto ss = walberla::make_shared<data::ShapeStorage>();
-   using ParticleAccessor_T = ParticleAccessorWithShape;
+   using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
    auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
 
    auto sphereShape = ss->create<data::Sphere>( radius );
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
index 25e47423b4c9f8a4d39c39913c6a7a3796f78047..9535269c5b505495e97655e9f9119a1fa986c8c0 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
@@ -20,7 +20,7 @@
 
 #include "mesa_pd/collision_detection/AnalyticContactDetection.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"
 
@@ -40,25 +40,6 @@ namespace dem_integrator_accuracy {
 using namespace walberla;
 using namespace walberla::mesa_pd;
 
-
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
-
 /*
  * Tests the integrator accuracy for a DEM simulation by comparing the given coefficient of restitution to the simulated one.
  * For that, the velocity after a single sphere-wall collision is divided by the initial velocity before the simulation.
@@ -100,7 +81,7 @@ int main( int argc, char** argv )
    //init data structures
    auto ps = walberla::make_shared<data::ParticleStorage>(2);
    auto ss = walberla::make_shared<data::ShapeStorage>();
-   using ParticleAccessor_T = ParticleAccessorWithShape;
+   using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
    auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
 
    auto sphereShape = ss->create<data::Sphere>( radius );
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
index 7bed1e98502c000fc8def474b688b5144dc77808..c1adc3102ef9d775e8b1bf952c87a9acdbcd73af 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
@@ -20,7 +20,7 @@
 
 #include "mesa_pd/collision_detection/AnalyticContactDetection.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"
 
@@ -39,25 +39,6 @@ namespace dem_integrator_accuracy {
 using namespace walberla;
 using namespace walberla::mesa_pd;
 
-
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
-
 /*
  * Tests the integrator accuracy for a DEM simulation by comparing the given coefficient of restitution to the simulated one.
  * For that, the velocity after a single sphere-wall collision is divided by the initial velocity before the simulation.
@@ -99,7 +80,7 @@ int main( int argc, char** argv )
    //init data structures
    auto ps = walberla::make_shared<data::ParticleStorage>(2);
    auto ss = walberla::make_shared<data::ShapeStorage>();
-   using ParticleAccessor_T = ParticleAccessorWithShape;
+   using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
    auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
 
    auto sphereShape = ss->create<data::Sphere>( radius );
diff --git a/tests/mesa_pd/kernel/DetectAndStoreContacts.cpp b/tests/mesa_pd/kernel/DetectAndStoreContacts.cpp
index 5e4c2db774a8f64fc54a455583defae9759f922d..1bd7a14ee63b84395a21683abc4064eda5dd743f 100644
--- a/tests/mesa_pd/kernel/DetectAndStoreContacts.cpp
+++ b/tests/mesa_pd/kernel/DetectAndStoreContacts.cpp
@@ -29,7 +29,7 @@
 #include <mesa_pd/domain/InfiniteDomain.h>
 #include <mesa_pd/kernel/ParticleSelector.h>
 
-#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleAccessorWithShape.h>
 #include <core/Environment.h>
 #include <core/logging/Logging.h>
 
@@ -38,24 +38,6 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
-
 int main( int argc, char ** argv )
 {
    Environment env(argc, argv);
@@ -67,7 +49,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  smallSphere = ss->create<data::Sphere>( real_t(1.2) );
 
diff --git a/tests/mesa_pd/kernel/ExplicitEuler.cpp b/tests/mesa_pd/kernel/ExplicitEuler.cpp
index b4673d1895067c3e6aa46b0d50a6de53cfa952e3..e73592b6ad7f5696c274f3ce2acd3ec1fef1e3b4 100644
--- a/tests/mesa_pd/kernel/ExplicitEuler.cpp
+++ b/tests/mesa_pd/kernel/ExplicitEuler.cpp
@@ -39,7 +39,9 @@ public:
    void setInvMass(const size_t /*p_idx*/, const walberla::real_t& v) { invMass_ = v;}
    const walberla::mesa_pd::Mat3& getInvInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_;}
    void setInvInertiaBF(const size_t /*p_idx*/, const walberla::mesa_pd::Mat3& v) { invInertiaBF_ = v;}
+   walberla::mesa_pd::Mat3 getInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_.getInverse();}
 
+private:
    walberla::real_t        invMass_;
    walberla::mesa_pd::Mat3 invInertiaBF_;
 };
diff --git a/tests/mesa_pd/kernel/GenerateAnalyticContacts.cpp b/tests/mesa_pd/kernel/GenerateAnalyticContacts.cpp
index fcd2b82c8ebaef80843686cada57d5d96a764a5a..81f245513e12f87cef763d066ac75e6f31726c85 100644
--- a/tests/mesa_pd/kernel/GenerateAnalyticContacts.cpp
+++ b/tests/mesa_pd/kernel/GenerateAnalyticContacts.cpp
@@ -20,7 +20,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>
@@ -36,8 +36,6 @@
 #include <blockforest/Initialization.h>
 #include <core/Environment.h>
 #include <core/grid_generator/SCIterator.h>
-#include <core/logging/Logging.h>
-#include <core/mpi/Reduce.h>
 
 #include <iostream>
 #include <memory>
@@ -45,23 +43,6 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
 int main( int argc, char ** argv )
 {
    Environment env(argc, argv);
@@ -78,7 +59,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(math::AABB(-1,-1,-1,4,4,4), real_t(1.3));
 
    //initialize particles
diff --git a/tests/mesa_pd/kernel/HCSITSKernels.cpp b/tests/mesa_pd/kernel/HCSITSKernels.cpp
index 297f6c25b46cddef4d6c3fbe15ef8c5023a63052..422be960e65afde6934c6cb5ebeab04e88a23fca 100644
--- a/tests/mesa_pd/kernel/HCSITSKernels.cpp
+++ b/tests/mesa_pd/kernel/HCSITSKernels.cpp
@@ -43,7 +43,7 @@
 #include "mesa_pd/mpi/notifications/VelocityUpdateNotification.h"
 #include "mesa_pd/mpi/notifications/VelocityCorrectionNotification.h"
 
-#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleAccessorWithShape.h>
 #include <core/Environment.h>
 #include <core/logging/Logging.h>
 
@@ -52,24 +52,6 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-      : ParticleAccessor(ps)
-      , ss_(ss)
-   {}
-
-   const real_t& getMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getMass();}
-   const real_t& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
-
-   const Mat3& getInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInertiaBF();}
-   const Mat3& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
 
 template<typename PStorage, typename CStorage, typename PAccessor, typename CAccessor>
 class TestHCSITSKernel {
@@ -141,7 +123,7 @@ void normalReactionTest(kernel::HCSITSRelaxationStep::RelaxationModel model)
    auto ps = std::make_shared<data::ParticleStorage>(100);
    auto cs = std::make_shared<data::ContactStorage>(100);
    auto ss = std::make_shared<data::ShapeStorage>();
-   ParticleAccessorWithShape paccessor(ps, ss);
+   data::ParticleAccessorWithShape paccessor(ps, ss);
    data::ContactAccessor caccessor(cs);
    auto density = real_t(7.874);
    auto radius = real_t(1.1);
@@ -171,7 +153,7 @@ void normalReactionTest(kernel::HCSITSRelaxationStep::RelaxationModel model)
    data::particle_flags::set(p2->getFlagsRef(), data::particle_flags::FIXED);
    data::particle_flags::set(p2->getFlagsRef(), data::particle_flags::GLOBAL);
 
-   TestHCSITSKernel<data::ParticleStorage, data::ContactStorage, ParticleAccessorWithShape, data::ContactAccessor> testHCSITS(*ps, *cs, paccessor, caccessor);
+   TestHCSITSKernel<data::ParticleStorage, data::ContactStorage, data::ParticleAccessorWithShape, data::ContactAccessor> testHCSITS(*ps, *cs, paccessor, caccessor);
    testHCSITS.model = model;
 
    WALBERLA_LOG_INFO(paccessor.getInvMass(0))
@@ -271,7 +253,7 @@ void SphereSphereTest(kernel::HCSITSRelaxationStep::RelaxationModel model){
    auto ps = std::make_shared<data::ParticleStorage>(100);
    auto cs = std::make_shared<data::ContactStorage>(100);
    auto ss = std::make_shared<data::ShapeStorage>();
-   ParticleAccessorWithShape paccessor(ps, ss);
+   data::ParticleAccessorWithShape paccessor(ps, ss);
    data::ContactAccessor caccessor(cs);
    auto density = real_t(7.874);
    auto radius = real_t(1.1);
@@ -294,7 +276,7 @@ void SphereSphereTest(kernel::HCSITSRelaxationStep::RelaxationModel model){
    p2->getOwnerRef()             = walberla::mpi::MPIManager::instance()->rank();
    p2->getLinearVelocityRef() = Vec3(real_t(-1), real_t(0), real_t(0));
    p2->getTypeRef()              = 0;
-   TestHCSITSKernel<data::ParticleStorage, data::ContactStorage, ParticleAccessorWithShape, data::ContactAccessor> testHCSITS(*ps, *cs, paccessor, caccessor);
+   TestHCSITSKernel<data::ParticleStorage, data::ContactStorage, data::ParticleAccessorWithShape, data::ContactAccessor> testHCSITS(*ps, *cs, paccessor, caccessor);
    testHCSITS.model = model;
    testHCSITS(dt);
 
@@ -322,7 +304,7 @@ void SphereSeperationTest(kernel::HCSITSRelaxationStep::RelaxationModel model){
    auto ps = std::make_shared<data::ParticleStorage>(100);
    auto cs = std::make_shared<data::ContactStorage>(100);
    auto ss = std::make_shared<data::ShapeStorage>();
-   ParticleAccessorWithShape paccessor(ps, ss);
+   data::ParticleAccessorWithShape paccessor(ps, ss);
    data::ContactAccessor caccessor(cs);
    auto density = real_t(7.874);
    auto radius = real_t(1.1);
@@ -345,7 +327,7 @@ void SphereSeperationTest(kernel::HCSITSRelaxationStep::RelaxationModel model){
    p2->getOwnerRef()             = walberla::mpi::MPIManager::instance()->rank();
    p2->getLinearVelocityRef()    = Vec3(real_t(-1), real_t(0), real_t(0));
    p2->getTypeRef()              = 0;
-   TestHCSITSKernel<data::ParticleStorage, data::ContactStorage, ParticleAccessorWithShape, data::ContactAccessor> testHCSITS(*ps, *cs, paccessor, caccessor);
+   TestHCSITSKernel<data::ParticleStorage, data::ContactStorage, data::ParticleAccessorWithShape, data::ContactAccessor> testHCSITS(*ps, *cs, paccessor, caccessor);
 
    int solveCount = 0;
    testHCSITS.model = model;
@@ -379,7 +361,7 @@ void SlidingSphereFrictionalReactionTest(kernel::HCSITSRelaxationStep::Relaxatio
    auto ps = std::make_shared<data::ParticleStorage>(100);
    auto cs = std::make_shared<data::ContactStorage>(100);
    auto ss = std::make_shared<data::ShapeStorage>();
-   ParticleAccessorWithShape paccessor(ps, ss);
+   data::ParticleAccessorWithShape paccessor(ps, ss);
    data::ContactAccessor caccessor(cs);
    auto density = real_t(1);
    auto radius = real_t(1);
@@ -406,7 +388,7 @@ void SlidingSphereFrictionalReactionTest(kernel::HCSITSRelaxationStep::Relaxatio
    data::particle_flags::set(p2->getFlagsRef(), data::particle_flags::FIXED);
    data::particle_flags::set(p2->getFlagsRef(), data::particle_flags::GLOBAL);
 
-   TestHCSITSKernel<data::ParticleStorage, data::ContactStorage, ParticleAccessorWithShape, data::ContactAccessor> testHCSITS(*ps, *cs, paccessor, caccessor);
+   TestHCSITSKernel<data::ParticleStorage, data::ContactStorage, data::ParticleAccessorWithShape, data::ContactAccessor> testHCSITS(*ps, *cs, paccessor, caccessor);
    testHCSITS.model = model;
    testHCSITS.globalAcc = Vec3(0,0,-10);
    int solveCount = 0;
diff --git a/tests/mesa_pd/kernel/IntegratorAccuracy.cpp b/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
index 433ecb5467e60032b94ed0b3ebf81586c98b44ec..0641e3d710265a3ec21176c3e61d31f92e65073d 100644
--- a/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
+++ b/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
@@ -43,15 +43,14 @@ 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 Mat3(real_t(0)); }
 
-   const auto &getInvInertiaBF(const size_t /*p_idx*/) const
-   { return invInertiaBF_; }
+   const auto getInertiaBF(const size_t /*p_idx*/) const // dummy
+   { return Mat3(real_t(0)); }
 
 private:
    real_t invMass_;
-   Mat3 invInertiaBF_;
 };
 
 struct Oscillator
diff --git a/tests/mesa_pd/kernel/LinearSpringDashpot.cpp b/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
index bf14293e0e1a27fcb6caf981ab80d50d6293e98d..cc7f2abb8411e47bda40aff8007b8f0df83a474a 100644
--- a/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
+++ b/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
@@ -21,7 +21,7 @@
 #include "mesa_pd/collision_detection/AnalyticContactDetection.h"
 #include "mesa_pd/common/ParticleFunctions.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"
 
@@ -39,24 +39,6 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
-
 /*
  * Simulates oblique sphere-wall collision and checks rebound angle, i.e. the tangential part of the collision model.
  *
@@ -106,7 +88,7 @@ int main( int argc, char ** argv )
    //init data structures
    auto ps = walberla::make_shared<data::ParticleStorage>(2);
    auto ss = walberla::make_shared<data::ShapeStorage>();
-   using ParticleAccessor_T = ParticleAccessorWithShape;
+   using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
    auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
 
    auto sphereShape = ss->create<data::Sphere>( radius );
diff --git a/tests/mesa_pd/kernel/LinkedCellsVsBruteForce.cpp b/tests/mesa_pd/kernel/LinkedCellsVsBruteForce.cpp
index 8743d6ac3bcf21b8a2d9cf82c6d3defed3405c01..40158448f305522a4372f3141412d9807ead8970 100644
--- a/tests/mesa_pd/kernel/LinkedCellsVsBruteForce.cpp
+++ b/tests/mesa_pd/kernel/LinkedCellsVsBruteForce.cpp
@@ -20,7 +20,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>
@@ -45,23 +45,6 @@
 namespace walberla {
 namespace mesa_pd {
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
 class comp
 {
 public:
@@ -104,7 +87,7 @@ int main( int argc, char ** argv )
    std::vector<collision_detection::AnalyticContactDetection> cs1(100);
    std::vector<collision_detection::AnalyticContactDetection> cs2(100);
 
-   ParticleAccessorWithShape accessor(ps, ss);
+   data::ParticleAccessorWithShape accessor(ps, ss);
 
    //initialize particles
    const real_t radius  = real_t(0.5);
diff --git a/tests/mesa_pd/kernel/SemiImplicitEuler.cpp b/tests/mesa_pd/kernel/SemiImplicitEuler.cpp
index 551d016f11c958d04f92db784f8f4daa9c94a35a..c12d8406c4c748cd2e0edfc2c0cafb56ad1eaafe 100644
--- a/tests/mesa_pd/kernel/SemiImplicitEuler.cpp
+++ b/tests/mesa_pd/kernel/SemiImplicitEuler.cpp
@@ -39,6 +39,7 @@ public:
    void setInvMass(const size_t /*p_idx*/, const walberla::real_t& v) { invMass_ = v;}
    const walberla::mesa_pd::Mat3& getInvInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_;}
    void setInvInertiaBF(const size_t /*p_idx*/, const walberla::mesa_pd::Mat3& v) { invInertiaBF_ = v;}
+   const walberla::mesa_pd::Mat3 getInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_.getInverse();}
 
    walberla::real_t        invMass_;
    walberla::mesa_pd::Mat3 invInertiaBF_;
diff --git a/tests/mesa_pd/kernel/SpringDashpot.cpp b/tests/mesa_pd/kernel/SpringDashpot.cpp
index 1e90c9dd6d6083ebb9e2f39c82c2e4e2e4ec5672..8d174af5209f048d1cd92646c146fe6f1f4ba7bb 100644
--- a/tests/mesa_pd/kernel/SpringDashpot.cpp
+++ b/tests/mesa_pd/kernel/SpringDashpot.cpp
@@ -20,7 +20,7 @@
 
 #include <mesa_pd/collision_detection/AnalyticContactDetection.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>
 
@@ -36,23 +36,6 @@ namespace walberla {
 
 using namespace walberla::mesa_pd;
 
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
-   ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
-         : ParticleAccessor(ps)
-         , ss_(ss)
-   {}
-
-   const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
-
-   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
-
-   data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
-private:
-   std::shared_ptr<data::ShapeStorage> ss_;
-};
-
 int main( int argc, char ** argv )
 {
    Environment env(argc, argv);
@@ -72,7 +55,7 @@ int main( int argc, char ** argv )
    auto smallSphere = ss->create<data::Sphere>( real_t(2) );
    ss->shapes[smallSphere]->updateMassAndInertia(real_t(2500));
 
-   ParticleAccessorWithShape ac(ps, ss);
+   mesa_pd::data::ParticleAccessorWithShape ac(ps, ss);
 
    data::Particle&& p1 = *ps->create();
    p1.getPositionRef() = Vec3(0,0,0);
diff --git a/tests/mesa_pd/kernel/VelocityVerlet.cpp b/tests/mesa_pd/kernel/VelocityVerlet.cpp
index d82f349c0c72a593685cf5ed226d04d59e87a984..21cd0a9184384590b391ec49272d7531830a57a3 100644
--- a/tests/mesa_pd/kernel/VelocityVerlet.cpp
+++ b/tests/mesa_pd/kernel/VelocityVerlet.cpp
@@ -18,7 +18,7 @@
 //
 //======================================================================================================================
 
-#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleAccessorWithShape.h>
 #include <mesa_pd/data/ParticleStorage.h>
 #include <mesa_pd/data/shape/Sphere.h>
 
@@ -48,6 +48,7 @@ public:
    const auto& getInvMass(const size_t /*p_idx*/) const {return sp.getInvMass();}
 
    const auto& getInvInertiaBF(const size_t /*p_idx*/) const {return sp.getInvInertiaBF();}
+   const auto& getInertiaBF(const size_t /*p_idx*/) const {return sp.getInertiaBF();}
 
    data::BaseShape* getShape(const size_t /*p_idx*/) {return &sp;}
 private:
diff --git a/tests/mesa_pd/kernel/cnt/SphericalSegmentAccessor.h b/tests/mesa_pd/kernel/cnt/SphericalSegmentAccessor.h
index 22c09b56486a3d4c6b9332b0b56c0a63bb66fdea..5b2090c599098955c233edc0f6531c511fc440aa 100644
--- a/tests/mesa_pd/kernel/cnt/SphericalSegmentAccessor.h
+++ b/tests/mesa_pd/kernel/cnt/SphericalSegmentAccessor.h
@@ -38,10 +38,13 @@ public:
    constexpr auto getInvMass(const size_t /*p_idx*/) const {return 1_r / kernel::cnt::mass_T;}
 
    constexpr auto& getInvInertiaBF(const size_t /*p_idx*/) const {return invI;}
+
+   constexpr auto& getInertiaBF(const size_t /*p_idx*/) const {return I;}
 private:
    //  - sphere   :  I = (2/5)*mass*radius^2
    static constexpr auto Ia = 0.4_r * kernel::cnt::mass_T * kernel::cnt::inner_radius * kernel::cnt::inner_radius;
    static constexpr auto invI = Mat3(1_r/Ia, 0_r, 0_r, 0_r, 1_r/Ia, 0_r, 0_r, 0_r, 1_r/Ia);
+   static constexpr auto I = Mat3(Ia, 0_r, 0_r, 0_r, Ia, 0_r, 0_r, 0_r, Ia);
 };
 
 } //namespace mesa_pd
diff --git a/tests/mesa_pd/kernel/interfaces/ExplicitEulerInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/ExplicitEulerInterfaceCheck.cpp
index 741d3ab6976744b4f2edd0ad083171d7d9b2472d..0731347afcf6c7124e883b7db5105bc69664f960 100644
--- a/tests/mesa_pd/kernel/interfaces/ExplicitEulerInterfaceCheck.cpp
+++ b/tests/mesa_pd/kernel/interfaces/ExplicitEulerInterfaceCheck.cpp
@@ -59,6 +59,8 @@ public:
    
    const walberla::mesa_pd::Mat3& getInvInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_;}
    
+   const walberla::mesa_pd::Mat3& getInertiaBF(const size_t /*p_idx*/) const {return inertiaBF_;}
+   
    const walberla::mesa_pd::Vec3& getTorque(const size_t /*p_idx*/) const {return torque_;}
    void setTorque(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { torque_ = v;}
    
@@ -81,6 +83,7 @@ private:
    walberla::mesa_pd::Rot3 rotation_;
    walberla::mesa_pd::Vec3 angularVelocity_;
    walberla::mesa_pd::Mat3 invInertiaBF_;
+   walberla::mesa_pd::Mat3 inertiaBF_;
    walberla::mesa_pd::Vec3 torque_;
 };
 
diff --git a/tests/mesa_pd/kernel/interfaces/ExplicitEulerWithShapeInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/ExplicitEulerWithShapeInterfaceCheck.cpp
index db83474b4ff801832004045124b2fa7759a5245f..76749788ae6c9d9ebc15f816ad9b58fb61767278 100644
--- a/tests/mesa_pd/kernel/interfaces/ExplicitEulerWithShapeInterfaceCheck.cpp
+++ b/tests/mesa_pd/kernel/interfaces/ExplicitEulerWithShapeInterfaceCheck.cpp
@@ -56,6 +56,8 @@ public:
    void setAngularVelocity(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { angularVelocity_ = v;}
    
    const walberla::mesa_pd::Mat3& getInvInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_;}
+
+   const walberla::mesa_pd::Mat3& getInertiaBF(const size_t /*p_idx*/) const {return inertiaBF_;}
    
    const walberla::mesa_pd::Vec3& getTorque(const size_t /*p_idx*/) const {return torque_;}
    void setTorque(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { torque_ = v;}
@@ -80,6 +82,7 @@ private:
    walberla::mesa_pd::Rot3 rotation_;
    walberla::mesa_pd::Vec3 angularVelocity_;
    walberla::mesa_pd::Mat3 invInertiaBF_;
+   walberla::mesa_pd::Mat3 inertiaBF_;
    walberla::mesa_pd::Vec3 torque_;
    walberla::mesa_pd::data::particle_flags::FlagT flags_;
 };
diff --git a/tests/mesa_pd/kernel/interfaces/SemiImplicitEulerInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/SemiImplicitEulerInterfaceCheck.cpp
index 5b287377b68bdea2ec67e2a1127f09d528ba3d69..4bb7530c0292a9af710817590fe9fffd5f247428 100644
--- a/tests/mesa_pd/kernel/interfaces/SemiImplicitEulerInterfaceCheck.cpp
+++ b/tests/mesa_pd/kernel/interfaces/SemiImplicitEulerInterfaceCheck.cpp
@@ -59,6 +59,8 @@ public:
    
    const walberla::mesa_pd::Mat3& getInvInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_;}
    
+   const walberla::mesa_pd::Mat3& getInertiaBF(const size_t /*p_idx*/) const {return inertiaBF_;}
+   
    const walberla::mesa_pd::Vec3& getTorque(const size_t /*p_idx*/) const {return torque_;}
    void setTorque(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { torque_ = v;}
    
@@ -81,6 +83,7 @@ private:
    walberla::mesa_pd::Rot3 rotation_;
    walberla::mesa_pd::Vec3 angularVelocity_;
    walberla::mesa_pd::Mat3 invInertiaBF_;
+   walberla::mesa_pd::Mat3 inertiaBF_;
    walberla::mesa_pd::Vec3 torque_;
 };
 
diff --git a/tests/mesa_pd/kernel/interfaces/VelocityVerletInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/VelocityVerletInterfaceCheck.cpp
index 869b4bb897cec19d051e2d46727c2f7cb1a9408d..bf171d8546e77e907bfae4680e9a4c50f76954f9 100644
--- a/tests/mesa_pd/kernel/interfaces/VelocityVerletInterfaceCheck.cpp
+++ b/tests/mesa_pd/kernel/interfaces/VelocityVerletInterfaceCheck.cpp
@@ -60,6 +60,8 @@ public:
    
    const walberla::mesa_pd::Mat3& getInvInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_;}
    
+   const walberla::mesa_pd::Mat3& getInertiaBF(const size_t /*p_idx*/) const {return inertiaBF_;}
+   
    const walberla::mesa_pd::Vec3& getTorque(const size_t /*p_idx*/) const {return torque_;}
    void setTorque(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { torque_ = v;}
    
@@ -87,6 +89,7 @@ private:
    walberla::mesa_pd::Rot3 rotation_;
    walberla::mesa_pd::Vec3 angularVelocity_;
    walberla::mesa_pd::Mat3 invInertiaBF_;
+   walberla::mesa_pd::Mat3 inertiaBF_;
    walberla::mesa_pd::Vec3 torque_;
    walberla::mesa_pd::Vec3 oldTorque_;
    walberla::mesa_pd::data::particle_flags::FlagT flags_;
diff --git a/tests/mesa_pd/kernel/interfaces/VelocityVerletWithShapeInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/VelocityVerletWithShapeInterfaceCheck.cpp
deleted file mode 100644
index 869b4bb897cec19d051e2d46727c2f7cb1a9408d..0000000000000000000000000000000000000000
--- a/tests/mesa_pd/kernel/interfaces/VelocityVerletWithShapeInterfaceCheck.cpp
+++ /dev/null
@@ -1,99 +0,0 @@
-//======================================================================================================================
-//
-//  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
-//! \author Sebastian Eibl <sebastian.eibl@fau.de>
-//
-//======================================================================================================================
-
-//======================================================================================================================
-//
-//  THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
-//
-//======================================================================================================================
-
-#include <mesa_pd/data/IAccessor.h>
-#include <mesa_pd/kernel/VelocityVerlet.h>
-
-#include <core/UniqueID.h>
-
-#include <map>
-
-namespace walberla {
-namespace mesa_pd {
-
-class Accessor : public data::IAccessor
-{
-public:
-   ~Accessor() override = default;
-   const walberla::mesa_pd::Vec3& getPosition(const size_t /*p_idx*/) const {return position_;}
-   void setPosition(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { position_ = v;}
-   
-   const walberla::mesa_pd::Vec3& getLinearVelocity(const size_t /*p_idx*/) const {return linearVelocity_;}
-   void setLinearVelocity(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { linearVelocity_ = v;}
-   
-   const walberla::real_t& getInvMass(const size_t /*p_idx*/) const {return invMass_;}
-   
-   const walberla::mesa_pd::Vec3& getForce(const size_t /*p_idx*/) const {return force_;}
-   void setForce(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { force_ = v;}
-   
-   const walberla::mesa_pd::Vec3& getOldForce(const size_t /*p_idx*/) const {return oldForce_;}
-   void setOldForce(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { oldForce_ = v;}
-   
-   const walberla::mesa_pd::Rot3& getRotation(const size_t /*p_idx*/) const {return rotation_;}
-   void setRotation(const size_t /*p_idx*/, const walberla::mesa_pd::Rot3& v) { rotation_ = v;}
-   
-   const walberla::mesa_pd::Vec3& getAngularVelocity(const size_t /*p_idx*/) const {return angularVelocity_;}
-   void setAngularVelocity(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { angularVelocity_ = v;}
-   
-   const walberla::mesa_pd::Mat3& getInvInertiaBF(const size_t /*p_idx*/) const {return invInertiaBF_;}
-   
-   const walberla::mesa_pd::Vec3& getTorque(const size_t /*p_idx*/) const {return torque_;}
-   void setTorque(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { torque_ = v;}
-   
-   const walberla::mesa_pd::Vec3& getOldTorque(const size_t /*p_idx*/) const {return oldTorque_;}
-   void setOldTorque(const size_t /*p_idx*/, const walberla::mesa_pd::Vec3& v) { oldTorque_ = v;}
-   
-   const walberla::mesa_pd::data::particle_flags::FlagT& getFlags(const size_t /*p_idx*/) const {return flags_;}
-   
-
-   id_t getInvalidUid() const {return UniqueID<int>::invalidID();}
-   size_t getInvalidIdx() const {return std::numeric_limits<size_t>::max();}
-   /**
-   * @brief Returns the index of particle specified by uid.
-   * @param uid unique id of the particle to be looked up
-   * @return the index of the particle or std::numeric_limits<size_t>::max() if the particle is not found
-   */
-   size_t uidToIdx(const id_t& /*uid*/) const {return 0;}
-   size_t size() const { return 1; }
-private:
-   walberla::mesa_pd::Vec3 position_;
-   walberla::mesa_pd::Vec3 linearVelocity_;
-   walberla::real_t invMass_;
-   walberla::mesa_pd::Vec3 force_;
-   walberla::mesa_pd::Vec3 oldForce_;
-   walberla::mesa_pd::Rot3 rotation_;
-   walberla::mesa_pd::Vec3 angularVelocity_;
-   walberla::mesa_pd::Mat3 invInertiaBF_;
-   walberla::mesa_pd::Vec3 torque_;
-   walberla::mesa_pd::Vec3 oldTorque_;
-   walberla::mesa_pd::data::particle_flags::FlagT flags_;
-};
-
-template void kernel::VelocityVerletPreForceUpdate::operator()(const size_t p_idx1, Accessor& ac) const;
-template void kernel::VelocityVerletPostForceUpdate::operator()(const size_t p_idx1, Accessor& ac) const;
-
-} //namespace mesa_pd
-} //namespace walberla
\ No newline at end of file
diff --git a/tests/mesa_pd/mpi/ReduceContactHistory.cpp b/tests/mesa_pd/mpi/ReduceContactHistory.cpp
index 0414ed75ddaf1ae0b77df03d39b682494cdae1e1..7d9ab310329b00b1e903f7170b7bee0b33524c95 100644
--- a/tests/mesa_pd/mpi/ReduceContactHistory.cpp
+++ b/tests/mesa_pd/mpi/ReduceContactHistory.cpp
@@ -21,10 +21,8 @@
 #include <mesa_pd/collision_detection/AnalyticContactDetection.h>
 #include <mesa_pd/data/ParticleAccessor.h>
 #include <mesa_pd/data/ParticleStorage.h>
-#include <mesa_pd/data/ShapeStorage.h>
 #include <mesa_pd/domain/BlockForestDomain.h>
 #include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
 #include <mesa_pd/kernel/ParticleSelector.h>
 #include <mesa_pd/mpi/ContactFilter.h>
 #include <mesa_pd/mpi/ReduceContactHistory.h>