From b67091edea71eb1dfcbcce486ea1af44dc2d8a6c Mon Sep 17 00:00:00 2001
From: Christoph Rettinger <christoph.rettinger@fau.de>
Date: Fri, 11 Mar 2022 13:58:14 +0100
Subject: [PATCH] Fix of integration of angular velocity of non-spherical
particles
---
.../GranularGas/MESA_PD_GranularGas.cpp | 7 +-
.../GranularGas/MESA_PD_KernelBenchmark.cpp | 5 +-
.../MESA_PD_KernelLoadBalancing.cpp | 5 +-
.../GranularGas/MESA_PD_LoadBalancing.cpp | 5 +-
.../IntegratorAccuracy/IntegratorAccuracy.cpp | 10 +-
apps/showcases/ParticlePacking/Evaluation.h | 1 +
.../ParticlePacking/ParticlePacking.cpp | 19 ++--
apps/showcases/ParticlePacking/Utility.h | 20 +---
python/mesa_pd/kernel/ExplicitEuler.py | 1 +
python/mesa_pd/kernel/SemiImplicitEuler.py | 1 +
python/mesa_pd/kernel/VelocityVerlet.py | 1 +
.../common/ParticleFunctions.templ.h | 15 +++
.../templates/kernel/ExplicitEuler.templ.h | 11 ++-
.../kernel/InitParticlesForHCSITS.templ.h | 5 +-
.../kernel/SemiImplicitEuler.templ.h | 11 ++-
.../templates/kernel/VelocityVerlet.templ.h | 22 +++--
src/mesa_pd/common/ParticleFunctions.h | 15 +++
.../data/ParticleAccessorWithBaseShape.h | 31 +++---
src/mesa_pd/data/ParticleAccessorWithShape.h | 2 +-
src/mesa_pd/kernel/ExplicitEuler.h | 11 ++-
src/mesa_pd/kernel/ExplicitEulerWithShape.h | 10 +-
src/mesa_pd/kernel/InitParticlesForHCSITS.h | 5 +-
src/mesa_pd/kernel/SemiImplicitEuler.h | 11 ++-
src/mesa_pd/kernel/VelocityVerlet.h | 22 +++--
tests/mesa_pd/ContactDetection.cpp | 21 +---
tests/mesa_pd/DropTestAnalytic.cpp | 21 +---
tests/mesa_pd/DropTestGeneral.cpp | 21 +---
.../AnalyticContactDetection.cpp | 23 +----
.../GeneralContactDetection.cpp | 19 +---
.../kernel/CoefficientOfRestitutionLSD.cpp | 22 +----
.../kernel/CoefficientOfRestitutionNLSD.cpp | 23 +----
.../kernel/CoefficientOfRestitutionSD.cpp | 23 +----
.../mesa_pd/kernel/DetectAndStoreContacts.cpp | 22 +----
tests/mesa_pd/kernel/ExplicitEuler.cpp | 2 +
.../kernel/GenerateAnalyticContacts.cpp | 23 +----
tests/mesa_pd/kernel/HCSITSKernels.cpp | 36 ++-----
tests/mesa_pd/kernel/IntegratorAccuracy.cpp | 9 +-
tests/mesa_pd/kernel/LinearSpringDashpot.cpp | 22 +----
.../kernel/LinkedCellsVsBruteForce.cpp | 21 +---
tests/mesa_pd/kernel/SemiImplicitEuler.cpp | 1 +
tests/mesa_pd/kernel/SpringDashpot.cpp | 21 +---
tests/mesa_pd/kernel/VelocityVerlet.cpp | 3 +-
.../kernel/cnt/SphericalSegmentAccessor.h | 3 +
.../ExplicitEulerInterfaceCheck.cpp | 3 +
.../ExplicitEulerWithShapeInterfaceCheck.cpp | 3 +
.../SemiImplicitEulerInterfaceCheck.cpp | 3 +
.../VelocityVerletInterfaceCheck.cpp | 3 +
.../VelocityVerletWithShapeInterfaceCheck.cpp | 99 -------------------
tests/mesa_pd/mpi/ReduceContactHistory.cpp | 2 -
49 files changed, 223 insertions(+), 472 deletions(-)
rename apps/benchmarks/GranularGas/Accessor.h => src/mesa_pd/data/ParticleAccessorWithBaseShape.h (58%)
delete mode 100644 tests/mesa_pd/kernel/interfaces/VelocityVerletWithShapeInterfaceCheck.cpp
diff --git a/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp b/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp
index 4ddbe5d8b..fe7b9ff8f 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 d15e0dfd5..6a7ce88b4 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 e3e11a43d..176388493 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 c327dbd26..b8360a14f 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 af2dc8a3a..7e8709ea8 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 930fec6ac..c1a531a6d 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 41f97f8b9..2c55d064f 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 5e0f0b4c6..86b6a24c7 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 1b51adb85..6d7b39fe1 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 f5b6632a1..6732a78cf 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 1a31554cb..1bb4ff445 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 599c8d58b..61629f059 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 d6893ffbe..4cb293517 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 9769d3999..48c167315 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 edc393536..528139501 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 08797d155..9d37ed22e 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 52b043d12..155ca3af1 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 1d34ce30e..e9de554b8 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 8c884d370..9f956839a 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 bf1bfa660..44afb7b9c 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 06ed0c232..982511b8c 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 2022a0b80..6f6f7b6e3 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 13963d4f5..bf8cec461 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 67a5f1dc2..93d64239d 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 651968772..e01eeb38d 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 827f0951f..01bea0bea 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 52495a885..08f2bb82f 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 911e7344d..3860c5801 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 1dcb49882..e56fdbb16 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 b425ece1b..2ce861a46 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 25e47423b..9535269c5 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 7bed1e985..c1adc3102 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 5e4c2db77..1bd7a14ee 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 b4673d189..e73592b6a 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 fcd2b82c8..81f245513 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 297f6c25b..422be960e 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 433ecb546..0641e3d71 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 bf14293e0..cc7f2abb8 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 8743d6ac3..40158448f 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 551d016f1..c12d8406c 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 1e90c9dd6..8d174af52 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 d82f349c0..21cd0a918 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 22c09b564..5b2090c59 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 741d3ab69..0731347af 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 db83474b4..76749788a 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 5b287377b..4bb7530c0 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 869b4bb89..bf171d854 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 869b4bb89..000000000
--- 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 0414ed75d..7d9ab3103 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>
--
GitLab