From a283e38fa2842869f8f42c315bd6582db5b7d140 Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Mon, 30 Mar 2020 13:30:24 +0200
Subject: [PATCH] [API] merge integrators with/without shape
To remove code duplication the distinction is now moved
into the generator. Default kernels are "with shape".
---
.../ObliqueDryCollision.cpp | 10 +-
.../ObliqueWetCollision.cpp | 10 +-
.../SettlingSphereInBox.cpp | 10 +-
.../SphereMovingWithPrescribedVelocity.cpp | 4 +-
.../SphereWallCollision.cpp | 10 +-
apps/benchmarks/GranularGas/GenerateModule.py | 3 +-
.../GranularGas/MESA_PD_GranularGas.cpp | 6 +-
.../GranularGas/MESA_PD_KernelBenchmark.cpp | 6 +-
.../MESA_PD_KernelLoadBalancing.cpp | 8 +-
.../GranularGas/MESA_PD_LoadBalancing.cpp | 8 +-
apps/tutorials/mesa_pd/01_MESAPD.cpp | 6 +-
python/mesa_pd.py | 2 -
python/mesa_pd/kernel/ExplicitEuler.py | 20 +-
.../mesa_pd/kernel/ExplicitEulerWithShape.py | 28 ---
python/mesa_pd/kernel/VelocityVerlet.py | 12 +-
.../mesa_pd/kernel/VelocityVerletWithShape.py | 39 ----
python/mesa_pd/kernel/__init__.py | 6 +-
.../templates/kernel/ExplicitEuler.templ.h | 26 ++-
.../kernel/ExplicitEulerWithShape.templ.h | 102 ---------
.../templates/kernel/VelocityVerlet.templ.h | 29 ++-
.../kernel/VelocityVerletWithShape.templ.h | 137 ------------
src/mesa_pd/kernel/ExplicitEuler.h | 32 ++-
src/mesa_pd/kernel/ExplicitEulerWithShape.h | 14 +-
src/mesa_pd/kernel/VelocityVerlet.h | 34 ++-
src/mesa_pd/kernel/VelocityVerletWithShape.h | 155 --------------
.../SettlingSphere.cpp | 10 +-
tests/mesa_pd/CMakeLists.txt | 6 -
tests/mesa_pd/DropTestAnalytic.cpp | 6 +-
tests/mesa_pd/DropTestGeneral.cpp | 6 +-
.../IntegratorWithShapeEvaluation.ipynb | 2 +-
.../kernel/CoefficientOfRestitutionLSD.cpp | 4 +-
.../kernel/CoefficientOfRestitutionNLSD.cpp | 4 +-
.../kernel/CoefficientOfRestitutionSD.cpp | 4 +-
tests/mesa_pd/kernel/ExplicitEuler.cpp | 9 +
.../mesa_pd/kernel/ExplicitEulerWithShape.cpp | 125 -----------
tests/mesa_pd/kernel/GenerateLinkedCells.cpp | 9 +-
tests/mesa_pd/kernel/IntegratorAccuracy.cpp | 8 +-
tests/mesa_pd/kernel/LinearSpringDashpot.cpp | 10 +-
tests/mesa_pd/kernel/SpherePile.cpp | 4 +-
tests/mesa_pd/kernel/VelocityVerlet.cpp | 111 ++++++++--
.../kernel/VelocityVerletWithShape.cpp | 197 ------------------
.../ExplicitEulerInterfaceCheck.cpp | 15 ++
.../ExplicitEulerWithShapeInterfaceCheck.cpp | 6 +-
.../VelocityVerletInterfaceCheck.cpp | 19 ++
.../VelocityVerletWithShapeInterfaceCheck.cpp | 8 +-
45 files changed, 364 insertions(+), 916 deletions(-)
delete mode 100644 python/mesa_pd/kernel/ExplicitEulerWithShape.py
delete mode 100644 python/mesa_pd/kernel/VelocityVerletWithShape.py
delete mode 100644 python/mesa_pd/templates/kernel/ExplicitEulerWithShape.templ.h
delete mode 100644 python/mesa_pd/templates/kernel/VelocityVerletWithShape.templ.h
delete mode 100644 src/mesa_pd/kernel/VelocityVerletWithShape.h
delete mode 100644 tests/mesa_pd/kernel/ExplicitEulerWithShape.cpp
delete mode 100644 tests/mesa_pd/kernel/VelocityVerletWithShape.cpp
diff --git a/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp b/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp
index 18adc7713..9314699f4 100644
--- a/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp
@@ -26,8 +26,8 @@
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/kernel/DoubleCast.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
-#include "mesa_pd/kernel/VelocityVerletWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
@@ -130,11 +130,11 @@ int main( int argc, char ** argv )
data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::FIXED);
// velocity verlet
- kernel::VelocityVerletWithShapePreForceUpdate vvPreForce( dt );
- kernel::VelocityVerletWithShapePostForceUpdate vvPostForce( dt );
+ kernel::VelocityVerletPreForceUpdate vvPreForce( dt );
+ kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
// explicit euler
- kernel::ExplicitEulerWithShape explEuler( dt );
+ kernel::ExplicitEuler explEuler( dt );
// collision response
collision_detection::AnalyticContactDetection acd;
diff --git a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
index dceb28ec8..cac87d6c7 100644
--- a/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/ObliqueWetCollision.cpp
@@ -80,10 +80,10 @@
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/domain/BlockForestDataHandling.h"
#include "mesa_pd/kernel/DoubleCast.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h"
-#include "mesa_pd/kernel/VelocityVerletWithShape.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
@@ -853,9 +853,9 @@ int main( int argc, char **argv )
syncCall();
real_t timeStepSizeRPD = real_t(1)/real_t(numRPDSubCycles);
- mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(timeStepSizeRPD);
- mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
- mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
+ mesa_pd::kernel::ExplicitEuler explEulerIntegrator(timeStepSizeRPD);
+ mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
+ mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
// linear model
mesa_pd::kernel::LinearSpringDashpot linearCollisionResponse(1);
diff --git a/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp b/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
index a0423c979..f1d6ac21e 100644
--- a/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SettlingSphereInBox.cpp
@@ -72,9 +72,9 @@
#include "mesa_pd/data/shape/Sphere.h"
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/kernel/DoubleCast.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/ParticleSelector.h"
-#include "mesa_pd/kernel/VelocityVerletWithShape.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
#include "mesa_pd/mpi/ContactFilter.h"
@@ -639,9 +639,9 @@ int main( int argc, char **argv )
syncCall();
- mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
- mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles));
- mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles));
+ mesa_pd::kernel::ExplicitEuler explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
+ mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles));
+ mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::mpi::ReduceProperty reduceProperty;
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp b/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
index 86dad9611..1be5bf122 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereMovingWithPrescribedVelocity.cpp
@@ -76,7 +76,7 @@
#include "mesa_pd/data/shape/Sphere.h"
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/kernel/DoubleCast.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
@@ -597,7 +597,7 @@ int main( int argc, char **argv )
syncCall();
- mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
+ mesa_pd::kernel::ExplicitEuler explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::mpi::ReduceProperty reduceProperty;
diff --git a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
index 676390ebe..cd84f67e4 100644
--- a/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
+++ b/apps/benchmarks/FluidParticleCoupling/SphereWallCollision.cpp
@@ -80,11 +80,11 @@
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/domain/BlockForestDataHandling.h"
#include "mesa_pd/kernel/DoubleCast.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/kernel/NonLinearSpringDashpot.h"
#include "mesa_pd/kernel/ParticleSelector.h"
-#include "mesa_pd/kernel/VelocityVerletWithShape.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/ClearNextNeighborSync.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
@@ -962,9 +962,9 @@ int main( int argc, char **argv )
syncCall();
real_t timeStepSizeRPD = real_t(1)/real_t(numRPDSubCycles);
- mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(timeStepSizeRPD);
- mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
- mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
+ mesa_pd::kernel::ExplicitEuler explEulerIntegrator(timeStepSizeRPD);
+ mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(timeStepSizeRPD);
+ mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(timeStepSizeRPD);
// linear model
mesa_pd::kernel::LinearSpringDashpot linearCollisionResponse(1);
diff --git a/apps/benchmarks/GranularGas/GenerateModule.py b/apps/benchmarks/GranularGas/GenerateModule.py
index 8971a0ad5..2ac9ebfb4 100755
--- a/apps/benchmarks/GranularGas/GenerateModule.py
+++ b/apps/benchmarks/GranularGas/GenerateModule.py
@@ -50,7 +50,6 @@ if __name__ == '__main__':
mpd.add(kernel.DoubleCast(ps))
mpd.add(kernel.ExplicitEuler())
- mpd.add(kernel.ExplicitEulerWithShape())
mpd.add(kernel.ForceLJ())
mpd.add(kernel.HeatConduction())
mpd.add(kernel.InsertParticleIntoLinkedCells())
@@ -60,7 +59,7 @@ if __name__ == '__main__':
mpd.add(kernel.SpringDashpot())
mpd.add(kernel.TemperatureIntegration())
mpd.add(kernel.VelocityVerlet())
- mpd.add(kernel.VelocityVerletWithShape())
+ mpd.add(kernel.VelocityVerlet())
mpd.add(mpi.BroadcastProperty())
mpd.add(mpi.ClearGhostOwnerSync())
diff --git a/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp b/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp
index 6712c5584..4ddbe5d8b 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp
+++ b/apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp
@@ -38,7 +38,7 @@
#include <mesa_pd/domain/BlockForestDomain.h>
#include <mesa_pd/kernel/AssocToBlock.h>
#include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
+#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
@@ -246,7 +246,7 @@ int main( int argc, char ** argv )
WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
// Init kernels
- kernel::ExplicitEulerWithShape explicitEulerWithShape( params.dt );
+ kernel::ExplicitEuler explicitEuler( params.dt );
DEM dem(domain, params.dt, ss->shapes[smallSphere]->getMass());
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::AssocToBlock assoc(forest);
@@ -316,7 +316,7 @@ int main( int argc, char ** argv )
tp["Euler"].start();
//ps->forEachParticle(false, [&](const size_t idx){WALBERLA_CHECK_EQUAL(ps->getForce(idx), Vec3(0,0,0), *(*ps)[idx] << "\n" << idx);});
- ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
+ ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
if (params.bBarrier) WALBERLA_MPI_BARRIER();
tp["Euler"].end();
diff --git a/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp b/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp
index 4523db176..d15e0dfd5 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp
+++ b/apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp
@@ -38,7 +38,7 @@
#include <mesa_pd/domain/BlockForestDomain.h>
#include <mesa_pd/kernel/AssocToBlock.h>
#include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
+#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
@@ -263,7 +263,7 @@ int main( int argc, char ** argv )
WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
// Init kernels
- kernel::ExplicitEulerWithShape explicitEulerWithShape( params.dt );
+ kernel::ExplicitEuler explicitEuler( params.dt );
kernel::AssocToBlock assoc(forest);
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::SpringDashpot sd(1);
@@ -398,7 +398,7 @@ int main( int argc, char ** argv )
tp["Euler"].start();
for (int64_t i=0; i < params.simulationSteps; ++i)
{
- ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
+ ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
}
tp["Euler"].end();
LIKWID_MARKER_STOP("Euler");
diff --git a/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp b/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp
index dfff27199..e3e11a43d 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp
+++ b/apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp
@@ -41,7 +41,7 @@
#include <mesa_pd/domain/InfoCollection.h>
#include <mesa_pd/kernel/AssocToBlock.h>
#include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
+#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
@@ -236,7 +236,7 @@ int main( int argc, char ** argv )
WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
// Init kernels
- kernel::ExplicitEulerWithShape explicitEulerWithShape( params.dt );
+ kernel::ExplicitEuler explicitEuler( params.dt );
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::SpringDashpot dem(1);
dem.setStiffness(0, 0, real_t(0));
@@ -334,7 +334,7 @@ int main( int argc, char ** argv )
tpImbalanced["Euler"].start();
for (int64_t i=0; i < params.simulationSteps; ++i)
{
- ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
+ ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
}
tpImbalanced["Euler"].end();
@@ -479,7 +479,7 @@ int main( int argc, char ** argv )
tpBalanced["Euler"].start();
for (int64_t i=0; i < params.simulationSteps; ++i)
{
- ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
+ ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
}
tpBalanced["Euler"].end();
diff --git a/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp b/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp
index 5c658a301..87b8dbaea 100644
--- a/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp
+++ b/apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp
@@ -41,7 +41,7 @@
#include <mesa_pd/domain/InfoCollection.h>
#include <mesa_pd/kernel/AssocToBlock.h>
#include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
+#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
@@ -230,7 +230,7 @@ int main( int argc, char ** argv )
WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
// Init kernels
- kernel::ExplicitEulerWithShape explicitEulerWithShape( params.dt );
+ kernel::ExplicitEuler explicitEuler( params.dt );
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::SpringDashpot dem(1);
dem.setStiffness(0, 0, real_t(0));
@@ -320,7 +320,7 @@ int main( int argc, char ** argv )
tpImbalanced["Euler"].start();
//ps->forEachParticle(false, [&](const size_t idx){WALBERLA_CHECK_EQUAL(ps->getForce(idx), Vec3(0,0,0), *(*ps)[idx] << "\n" << idx);});
- ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
+ ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
if (params.bBarrier) WALBERLA_MPI_BARRIER();
tpImbalanced["Euler"].end();
@@ -418,7 +418,7 @@ int main( int argc, char ** argv )
tpBalanced["Euler"].start();
//ps->forEachParticle(false, [&](const size_t idx){WALBERLA_CHECK_EQUAL(ps->getForce(idx), Vec3(0,0,0), *(*ps)[idx] << "\n" << idx);});
- ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
+ ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
if (params.bBarrier) WALBERLA_MPI_BARRIER();
tpBalanced["Euler"].end();
diff --git a/apps/tutorials/mesa_pd/01_MESAPD.cpp b/apps/tutorials/mesa_pd/01_MESAPD.cpp
index 31d21bead..9e952cc8f 100644
--- a/apps/tutorials/mesa_pd/01_MESAPD.cpp
+++ b/apps/tutorials/mesa_pd/01_MESAPD.cpp
@@ -36,7 +36,7 @@
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/domain/BlockForestDomain.h>
#include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
+#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
@@ -173,7 +173,7 @@ int main( int argc, char ** argv )
WALBERLA_LOG_INFO_ON_ROOT("*** KERNELS ***");
//! [Kernels]
- kernel::ExplicitEulerWithShape explicitEulerWithShape( dt );
+ kernel::ExplicitEuler explicitEuler( dt );
kernel::InsertParticleIntoLinkedCells ipilc;
kernel::SpringDashpot dem( 2 );
auto meffSphereSphere = real_t(0.5) * (real_c(4.0)/real_c(3.0) * math::pi) * radius * radius * radius * density;
@@ -218,7 +218,7 @@ int main( int argc, char ** argv )
RP.operator()<ForceTorqueNotification>(*ps);
- ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
+ ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
SNN(*ps, domain);
}
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index 599c0e075..785d9a4ff 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -85,7 +85,6 @@ if __name__ == '__main__':
mpd.add(kernel.DetectAndStoreContacts())
mpd.add(kernel.DoubleCast(ps))
mpd.add(kernel.ExplicitEuler())
- mpd.add(kernel.ExplicitEulerWithShape())
mpd.add(kernel.ForceLJ())
mpd.add(kernel.HCSITSRelaxationStep())
mpd.add(kernel.HeatConduction())
@@ -101,7 +100,6 @@ if __name__ == '__main__':
mpd.add(kernel.SpringDashpotSpring())
mpd.add(kernel.TemperatureIntegration())
mpd.add(kernel.VelocityVerlet())
- mpd.add(kernel.VelocityVerletWithShape())
mpd.add(mpi.BroadcastProperty())
mpd.add(mpi.ClearGhostOwnerSync())
diff --git a/python/mesa_pd/kernel/ExplicitEuler.py b/python/mesa_pd/kernel/ExplicitEuler.py
index b8cbb0eb2..6a18b0e50 100644
--- a/python/mesa_pd/kernel/ExplicitEuler.py
+++ b/python/mesa_pd/kernel/ExplicitEuler.py
@@ -5,13 +5,19 @@ from mesa_pd.utility import generate_file
class ExplicitEuler:
- def __init__(self):
- self.context = dict()
- self.context['interface'] = [create_access("position", "walberla::mesa_pd::Vec3", access="gs"),
- create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"),
- create_access("invMass", "walberla::real_t", access="g"),
- create_access("force", "walberla::mesa_pd::Vec3", access="gs"),
- create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g")]
+ def __init__(self, integrate_rotation=True):
+ self.context = {'bIntegrateRotation': integrate_rotation, 'interface': []}
+ self.context['interface'].append(create_access("position", "walberla::mesa_pd::Vec3", access="gs"))
+ self.context['interface'].append(create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"))
+ self.context['interface'].append(create_access("invMass", "walberla::real_t", access="g"))
+ self.context['interface'].append(create_access("force", "walberla::mesa_pd::Vec3", access="gs"))
+ self.context['interface'].append(create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g"))
+
+ if integrate_rotation:
+ 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("torque", "walberla::mesa_pd::Vec3", access="gs"))
def generate(self, module):
ctx = {'module': module, **self.context}
diff --git a/python/mesa_pd/kernel/ExplicitEulerWithShape.py b/python/mesa_pd/kernel/ExplicitEulerWithShape.py
deleted file mode 100644
index 9179e7b70..000000000
--- a/python/mesa_pd/kernel/ExplicitEulerWithShape.py
+++ /dev/null
@@ -1,28 +0,0 @@
-# -*- coding: utf-8 -*-
-
-from mesa_pd.accessor import create_access
-from mesa_pd.utility import generate_file
-
-
-class ExplicitEulerWithShape:
- def __init__(self):
- self.context = dict()
- self.context['interface'] = [create_access("position", "walberla::mesa_pd::Vec3", access="gs"),
- create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"),
- create_access("invMass", "walberla::real_t", access="g"),
- create_access("force", "walberla::mesa_pd::Vec3", access="gs"),
- create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"),
- create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"),
- create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"),
- create_access("torque", "walberla::mesa_pd::Vec3", access="gs"),
- create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g")]
-
- def generate(self, module):
- ctx = {'module': module, **self.context}
- generate_file(module['module_path'], 'kernel/ExplicitEulerWithShape.templ.h', ctx)
-
- ctx["InterfaceTestName"] = "ExplicitEulerWithShapeInterfaceCheck"
- ctx["KernelInclude"] = "kernel/ExplicitEulerWithShape.h"
- ctx["ExplicitInstantiation"] = "template void kernel::ExplicitEulerWithShape::operator()(const size_t p_idx1, Accessor& ac) const;"
- generate_file(module['test_path'], 'tests/CheckInterface.templ.cpp', ctx,
- 'kernel/interfaces/ExplicitEulerWithShapeInterfaceCheck.cpp')
diff --git a/python/mesa_pd/kernel/VelocityVerlet.py b/python/mesa_pd/kernel/VelocityVerlet.py
index 8d2d42a35..74eca6c13 100644
--- a/python/mesa_pd/kernel/VelocityVerlet.py
+++ b/python/mesa_pd/kernel/VelocityVerlet.py
@@ -5,13 +5,21 @@ from mesa_pd.utility import generate_file
class VelocityVerlet:
- def __init__(self):
- self.context = {'interface': []}
+ def __init__(self, integrate_rotation = True):
+ self.context = {'bIntegrateRotation': integrate_rotation, 'interface': []}
self.context['interface'].append(create_access("position", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("invMass", "walberla::real_t", access="g"))
self.context['interface'].append(create_access("force", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("oldForce", "walberla::mesa_pd::Vec3", access="gs"))
+
+ if integrate_rotation:
+ 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("torque", "walberla::mesa_pd::Vec3", access="gs"))
+ self.context['interface'].append(create_access("oldTorque", "walberla::mesa_pd::Vec3", access="gs"))
+
self.context['interface'].append(
create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g"))
diff --git a/python/mesa_pd/kernel/VelocityVerletWithShape.py b/python/mesa_pd/kernel/VelocityVerletWithShape.py
deleted file mode 100644
index e3d21c595..000000000
--- a/python/mesa_pd/kernel/VelocityVerletWithShape.py
+++ /dev/null
@@ -1,39 +0,0 @@
-# -*- coding: utf-8 -*-
-
-from mesa_pd.accessor import create_access
-from mesa_pd.utility import generate_file
-
-
-class VelocityVerletWithShape:
- def __init__(self):
- self.context = {'interface': []}
- self.context['interface'].append(create_access("position", "walberla::mesa_pd::Vec3", access="gs"))
- self.context['interface'].append(create_access("linearVelocity", "walberla::mesa_pd::Vec3", access="gs"))
- self.context['interface'].append(create_access("invMass", "walberla::real_t", access="g"))
- self.context['interface'].append(create_access("force", "walberla::mesa_pd::Vec3", access="gs"))
- self.context['interface'].append(create_access("oldForce", "walberla::mesa_pd::Vec3", access="gs"))
-
- 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("torque", "walberla::mesa_pd::Vec3", access="gs"))
- self.context['interface'].append(create_access("oldTorque", "walberla::mesa_pd::Vec3", access="gs"))
-
- self.context['interface'].append(
- create_access("flags", "walberla::mesa_pd::data::particle_flags::FlagT", access="g"))
-
- def getRequirements(self):
- return self.accessor
-
- def generate(self, module):
- ctx = {'module': module, **self.context}
-
- generate_file(module['module_path'], 'kernel/VelocityVerletWithShape.templ.h', ctx)
-
- ctx["InterfaceTestName"] = "VelocityVerletWithShapeInterfaceCheck"
- ctx["KernelInclude"] = "kernel/VelocityVerletWithShape.h"
- ctx[
- "ExplicitInstantiation"] = "template void kernel::VelocityVerletWithShapePreForceUpdate::operator()(const size_t p_idx1, Accessor& ac) const;\n" + \
- "template void kernel::VelocityVerletWithShapePostForceUpdate::operator()(const size_t p_idx1, Accessor& ac) const;"
- generate_file(module['test_path'], 'tests/CheckInterface.templ.cpp', ctx,
- 'kernel/interfaces/VelocityVerletWithShapeInterfaceCheck.cpp')
diff --git a/python/mesa_pd/kernel/__init__.py b/python/mesa_pd/kernel/__init__.py
index 8e79c416d..317646d1e 100644
--- a/python/mesa_pd/kernel/__init__.py
+++ b/python/mesa_pd/kernel/__init__.py
@@ -2,7 +2,6 @@
from .DetectAndStoreContacts import DetectAndStoreContacts
from .DoubleCast import DoubleCast
from .ExplicitEuler import ExplicitEuler
-from .ExplicitEulerWithShape import ExplicitEulerWithShape
from .ForceLJ import ForceLJ
from .HCSITSRelaxationStep import HCSITSRelaxationStep
from .HeatConduction import HeatConduction
@@ -18,12 +17,10 @@ from .SpringDashpot import SpringDashpot
from .SpringDashpotSpring import SpringDashpotSpring
from .TemperatureIntegration import TemperatureIntegration
from .VelocityVerlet import VelocityVerlet
-from .VelocityVerletWithShape import VelocityVerletWithShape
__all__ = ['DoubleCast',
'DetectAndStoreContacts',
'ExplicitEuler',
- 'ExplicitEulerWithShape',
'ForceLJ',
'HCSITSRelaxationStep',
'HeatConduction',
@@ -38,5 +35,4 @@ __all__ = ['DoubleCast',
'SpringDashpot',
'SpringDashpotSpring',
'TemperatureIntegration',
- 'VelocityVerlet',
- 'VelocityVerletWithShape']
+ 'VelocityVerlet']
diff --git a/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h b/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h
index a81e274f3..e3916291b 100644
--- a/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h
+++ b/python/mesa_pd/templates/kernel/ExplicitEuler.templ.h
@@ -35,7 +35,6 @@ namespace kernel {
/**
* Kernel which explicitly integrates all particles in time.
- * This integrator integrates velocity and position.
*
* This kernel requires the following particle accessor interface
* \code
@@ -53,8 +52,8 @@ namespace kernel {
{%- endfor %}
* \endcode
*
- * \pre All forces acting on the particles have to be set.
- * \post All forces are reset to 0.
+ * \pre All forces and torques acting on the particles have to be set.
+ * \post All forces and torques are reset to 0.
* \ingroup mesa_pd_kernel
*/
class ExplicitEuler
@@ -78,8 +77,27 @@ 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));
+
+ {%- if bIntegrateRotation %}
+ const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
+ ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
+
+ // Calculating the rotation angle
+ const Vec3 phi( ac.getAngularVelocity(idx) * dt_ + wdot * dt_ * dt_);
+
+ // Calculating the new orientation
+ auto rotation = ac.getRotation(idx);
+ rotation.rotate( phi );
+ ac.setRotation(idx, rotation);
+
+ ac.setAngularVelocity(idx, wdot * dt_ + ac.getAngularVelocity(idx));
+ {%- endif %}
}
- ac.setForce (idx, Vec3(real_t(0), real_t(0), real_t(0)));
+
+ ac.setForce (idx, Vec3(real_t(0), real_t(0), real_t(0)));
+ {%- if bIntegrateRotation %}
+ ac.setTorque(idx, Vec3(real_t(0), real_t(0), real_t(0)));
+ {%- endif %}
}
} //namespace kernel
diff --git a/python/mesa_pd/templates/kernel/ExplicitEulerWithShape.templ.h b/python/mesa_pd/templates/kernel/ExplicitEulerWithShape.templ.h
deleted file mode 100644
index 7fb2d73bc..000000000
--- a/python/mesa_pd/templates/kernel/ExplicitEulerWithShape.templ.h
+++ /dev/null
@@ -1,102 +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 ExplicitEulerWithShape.h
-//! \author Sebastian Eibl <sebastian.eibl@fau.de>
-//
-//======================================================================================================================
-
-//======================================================================================================================
-//
-// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
-//
-//======================================================================================================================
-
-#pragma once
-
-#include <mesa_pd/data/DataTypes.h>
-#include <mesa_pd/data/IAccessor.h>
-
-namespace walberla {
-namespace mesa_pd {
-namespace kernel {
-
-/**
- * Kernel which explicitly integrates all particles in time.
- * This integrator integrates velocity and position as well as angular velocity and rotation.
- *
- * This kernel requires the following particle accessor interface
- * \code
- {%- for prop in interface %}
- {%- if 'g' in prop.access %}
- * const {{prop.type}}& get{{prop.name | capFirst}}(const size_t p_idx) const;
- {%- endif %}
- {%- if 's' in prop.access %}
- * void set{{prop.name | capFirst}}(const size_t p_idx, const {{prop.type}}& v);
- {%- endif %}
- {%- if 'r' in prop.access %}
- * {{prop.type}}& get{{prop.name | capFirst}}Ref(const size_t p_idx);
- {%- endif %}
- *
- {%- endfor %}
- * \endcode
- *
- * \pre All forces and torques acting on the particles have to be set.
- * \post All forces and torques are reset to 0.
- * \ingroup mesa_pd_kernel
- */
-class ExplicitEulerWithShape
-{
-public:
- explicit ExplicitEulerWithShape(const real_t dt) : dt_(dt) {}
-
- template <typename Accessor>
- void operator()(const size_t i, Accessor& ac) const;
-private:
- real_t dt_ = real_t(0.0);
-};
-
-template <typename Accessor>
-inline void ExplicitEulerWithShape::operator()(const size_t idx,
- Accessor& ac) const
-{
- static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
-
- if (!data::particle_flags::isSet( ac.getFlags(idx), data::particle_flags::FIXED))
- {
- 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);
-
- // Calculating the rotation angle
- const Vec3 phi( ac.getAngularVelocity(idx) * dt_ + wdot * dt_ * dt_);
-
- // Calculating the new orientation
- auto rotation = ac.getRotation(idx);
- rotation.rotate( phi );
- ac.setRotation(idx, rotation);
-
- ac.setAngularVelocity(idx, wdot * dt_ + ac.getAngularVelocity(idx));
- }
-
- ac.setForce (idx, Vec3(real_t(0), real_t(0), real_t(0)));
- ac.setTorque(idx, Vec3(real_t(0), real_t(0), real_t(0)));
-}
-
-} //namespace kernel
-} //namespace mesa_pd
-} //namespace walberla
diff --git a/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h b/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h
index a3df53325..148ac8c8d 100644
--- a/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h
+++ b/python/mesa_pd/templates/kernel/VelocityVerlet.templ.h
@@ -40,7 +40,6 @@ namespace kernel {
* called before the force calculation and postFroceUpdate afterwards. The
* integration is only complete when both functions are called. The integration
* is symplectic.
- * \attention The force calculation has to be independent of velocity.
*
* This kernel requires the following particle accessor interface
* \code
@@ -92,6 +91,19 @@ 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_);
+
+ {%- if bIntegrateRotation %}
+ const Vec3 wdot = math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(),
+ ac.getInvInertiaBF(p_idx)) * ac.getOldTorque(p_idx);
+
+ // Calculating the rotation angle
+ const Vec3 phi( ac.getAngularVelocity(p_idx) * dt_ + real_t(0.5) * wdot * dt_ * dt_);
+
+ // Calculating the new orientation
+ auto rotation = ac.getRotation(p_idx);
+ rotation.rotate( phi );
+ ac.setRotation(p_idx, rotation);
+ {%- endif %}
}
}
@@ -104,9 +116,24 @@ 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_);
+
+ {%- 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;
+
+ ac.setAngularVelocity(p_idx, ac.getAngularVelocity(p_idx) +
+ real_t(0.5) * wdot * dt_ );
+ {%- endif %}
}
+
ac.setOldForce(p_idx, ac.getForce(p_idx));
ac.setForce(p_idx, Vec3(real_t(0), real_t(0), real_t(0)));
+
+ {%- if bIntegrateRotation %}
+ ac.setOldTorque(p_idx, ac.getTorque(p_idx));
+ ac.setTorque(p_idx, Vec3(real_t(0), real_t(0), real_t(0)));
+ {%- endif %}
}
} //namespace kernel
diff --git a/python/mesa_pd/templates/kernel/VelocityVerletWithShape.templ.h b/python/mesa_pd/templates/kernel/VelocityVerletWithShape.templ.h
deleted file mode 100644
index c166e1be5..000000000
--- a/python/mesa_pd/templates/kernel/VelocityVerletWithShape.templ.h
+++ /dev/null
@@ -1,137 +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 VelocityVerletWithShape.h
-//! \author Sebastian Eibl <sebastian.eibl@fau.de>
-//
-//======================================================================================================================
-
-//======================================================================================================================
-//
-// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
-//
-//======================================================================================================================
-
-#pragma once
-
-#include <mesa_pd/data/DataTypes.h>
-#include <mesa_pd/data/IAccessor.h>
-
-namespace walberla {
-namespace mesa_pd {
-namespace kernel {
-
-/**
- * Velocity verlet integration for all particles.
- *
- * Velocit verlet integration is a two part kernel. preForceUpdate has to be
- * called before the force calculation and postFroceUpdate afterwards. The
- * integration is only complete when both functions are called. The integration
- * is symplectic.
- * \attention The force calculation has to be independent of velocity.
- *
- * This kernel requires the following particle accessor interface
- * \code
- {%- for prop in interface %}
- {%- if 'g' in prop.access %}
- * const {{prop.type}}& get{{prop.name | capFirst}}(const size_t p_idx) const;
- {%- endif %}
- {%- if 's' in prop.access %}
- * void set{{prop.name | capFirst}}(const size_t p_idx, const {{prop.type}}& v);
- {%- endif %}
- {%- if 'r' in prop.access %}
- * {{prop.type}}& get{{prop.name | capFirst}}Ref(const size_t p_idx);
- {%- endif %}
- *
- {%- endfor %}
- * \endcode
- * \ingroup mesa_pd_kernel
- */
-class VelocityVerletWithShapePreForceUpdate
-{
-public:
- VelocityVerletWithShapePreForceUpdate(const real_t dt) : dt_(dt) {}
-
- template <typename Accessor>
- void operator()(const size_t i, Accessor& ac) const;
-
- real_t dt_;
-};
-
-/// \see VelocityVerletPreForceUpdate
-class VelocityVerletWithShapePostForceUpdate
-{
-public:
- VelocityVerletWithShapePostForceUpdate(const real_t dt) : dt_(dt) {}
-
- template <typename Accessor>
- void operator()(const size_t i, Accessor& ac) const;
-
- real_t dt_;
-};
-
-template <typename Accessor>
-inline void VelocityVerletWithShapePreForceUpdate::operator()(const size_t p_idx, Accessor& ac) const
-{
- static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
-
- if (!data::particle_flags::isSet( ac.getFlags(p_idx), data::particle_flags::FIXED))
- {
- 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);
-
- // Calculating the rotation angle
- const Vec3 phi( ac.getAngularVelocity(p_idx) * dt_ + real_t(0.5) * wdot * dt_ * dt_);
-
- // Calculating the new orientation
- auto rotation = ac.getRotation(p_idx);
- rotation.rotate( phi );
- ac.setRotation(p_idx, rotation);
- }
-}
-
-template <typename Accessor>
-inline void VelocityVerletWithShapePostForceUpdate::operator()(const size_t p_idx, Accessor& ac) const
-{
- static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
-
- if (!data::particle_flags::isSet( ac.getFlags(p_idx), data::particle_flags::FIXED))
- {
- 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;
-
- ac.setAngularVelocity(p_idx, ac.getAngularVelocity(p_idx) +
- real_t(0.5) * wdot * dt_ );
- }
-
- ac.setOldForce(p_idx, ac.getForce(p_idx));
- ac.setForce(p_idx, Vec3(real_t(0), real_t(0), real_t(0)));
-
- ac.setOldTorque(p_idx, ac.getTorque(p_idx));
- ac.setTorque(p_idx, Vec3(real_t(0), real_t(0), real_t(0)));
-}
-
-} //namespace kernel
-} //namespace mesa_pd
-} //namespace walberla
diff --git a/src/mesa_pd/kernel/ExplicitEuler.h b/src/mesa_pd/kernel/ExplicitEuler.h
index 47e93bc62..8887219b1 100644
--- a/src/mesa_pd/kernel/ExplicitEuler.h
+++ b/src/mesa_pd/kernel/ExplicitEuler.h
@@ -35,7 +35,6 @@ namespace kernel {
/**
* Kernel which explicitly integrates all particles in time.
- * This integrator integrates velocity and position.
*
* This kernel requires the following particle accessor interface
* \code
@@ -52,10 +51,21 @@ namespace kernel {
*
* const walberla::mesa_pd::data::particle_flags::FlagT& getFlags(const size_t p_idx) const;
*
+ * const walberla::mesa_pd::Rot3& getRotation(const size_t p_idx) const;
+ * void setRotation(const size_t p_idx, const walberla::mesa_pd::Rot3& v);
+ *
+ * const walberla::mesa_pd::Vec3& getAngularVelocity(const size_t p_idx) const;
+ * void setAngularVelocity(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
+ *
+ * const walberla::mesa_pd::Mat3& getInvInertiaBF(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);
+ *
* \endcode
*
- * \pre All forces acting on the particles have to be set.
- * \post All forces are reset to 0.
+ * \pre All forces and torques acting on the particles have to be set.
+ * \post All forces and torques are reset to 0.
* \ingroup mesa_pd_kernel
*/
class ExplicitEuler
@@ -79,8 +89,22 @@ 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);
+
+ // Calculating the rotation angle
+ const Vec3 phi( ac.getAngularVelocity(idx) * dt_ + wdot * dt_ * dt_);
+
+ // Calculating the new orientation
+ auto rotation = ac.getRotation(idx);
+ rotation.rotate( phi );
+ ac.setRotation(idx, rotation);
+
+ ac.setAngularVelocity(idx, wdot * dt_ + ac.getAngularVelocity(idx));
}
- ac.setForce (idx, Vec3(real_t(0), real_t(0), real_t(0)));
+
+ ac.setForce (idx, Vec3(real_t(0), real_t(0), real_t(0)));
+ ac.setTorque(idx, Vec3(real_t(0), real_t(0), real_t(0)));
}
} //namespace kernel
diff --git a/src/mesa_pd/kernel/ExplicitEulerWithShape.h b/src/mesa_pd/kernel/ExplicitEulerWithShape.h
index 6746d156e..d66d9a8aa 100644
--- a/src/mesa_pd/kernel/ExplicitEulerWithShape.h
+++ b/src/mesa_pd/kernel/ExplicitEulerWithShape.h
@@ -13,7 +13,7 @@
// 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 ExplicitEulerWithShape.h
+//! \file ExplicitEuler.h
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//
//======================================================================================================================
@@ -35,7 +35,6 @@ namespace kernel {
/**
* Kernel which explicitly integrates all particles in time.
- * This integrator integrates velocity and position as well as angular velocity and rotation.
*
* This kernel requires the following particle accessor interface
* \code
@@ -50,6 +49,8 @@ namespace kernel {
* const walberla::mesa_pd::Vec3& getForce(const size_t p_idx) const;
* void setForce(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
*
+ * const walberla::mesa_pd::data::particle_flags::FlagT& getFlags(const size_t p_idx) const;
+ *
* const walberla::mesa_pd::Rot3& getRotation(const size_t p_idx) const;
* void setRotation(const size_t p_idx, const walberla::mesa_pd::Rot3& v);
*
@@ -61,18 +62,16 @@ namespace kernel {
* 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);
*
- * const walberla::mesa_pd::data::particle_flags::FlagT& getFlags(const size_t p_idx) const;
- *
* \endcode
*
* \pre All forces and torques acting on the particles have to be set.
* \post All forces and torques are reset to 0.
* \ingroup mesa_pd_kernel
*/
-class ExplicitEulerWithShape
+class ExplicitEuler
{
public:
- explicit ExplicitEulerWithShape(const real_t dt) : dt_(dt) {}
+ explicit ExplicitEuler(const real_t dt) : dt_(dt) {}
template <typename Accessor>
void operator()(const size_t i, Accessor& ac) const;
@@ -81,7 +80,7 @@ private:
};
template <typename Accessor>
-inline void ExplicitEulerWithShape::operator()(const size_t idx,
+inline void ExplicitEuler::operator()(const size_t idx,
Accessor& ac) const
{
static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
@@ -90,7 +89,6 @@ inline void ExplicitEulerWithShape::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);
diff --git a/src/mesa_pd/kernel/VelocityVerlet.h b/src/mesa_pd/kernel/VelocityVerlet.h
index f4874bff0..f029652c2 100644
--- a/src/mesa_pd/kernel/VelocityVerlet.h
+++ b/src/mesa_pd/kernel/VelocityVerlet.h
@@ -40,7 +40,6 @@ namespace kernel {
* called before the force calculation and postFroceUpdate afterwards. The
* integration is only complete when both functions are called. The integration
* is symplectic.
- * \attention The force calculation has to be independent of velocity.
*
* This kernel requires the following particle accessor interface
* \code
@@ -58,6 +57,20 @@ namespace kernel {
* const walberla::mesa_pd::Vec3& getOldForce(const size_t p_idx) const;
* void setOldForce(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
*
+ * const walberla::mesa_pd::Rot3& getRotation(const size_t p_idx) const;
+ * void setRotation(const size_t p_idx, const walberla::mesa_pd::Rot3& v);
+ *
+ * const walberla::mesa_pd::Vec3& getAngularVelocity(const size_t p_idx) const;
+ * void setAngularVelocity(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
+ *
+ * const walberla::mesa_pd::Mat3& getInvInertiaBF(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);
+ *
+ * const walberla::mesa_pd::Vec3& getOldTorque(const size_t p_idx) const;
+ * void setOldTorque(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
+ *
* const walberla::mesa_pd::data::particle_flags::FlagT& getFlags(const size_t p_idx) const;
*
* \endcode
@@ -96,6 +109,16 @@ 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);
+
+ // Calculating the rotation angle
+ const Vec3 phi( ac.getAngularVelocity(p_idx) * dt_ + real_t(0.5) * wdot * dt_ * dt_);
+
+ // Calculating the new orientation
+ auto rotation = ac.getRotation(p_idx);
+ rotation.rotate( phi );
+ ac.setRotation(p_idx, rotation);
}
}
@@ -108,9 +131,18 @@ 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;
+
+ ac.setAngularVelocity(p_idx, ac.getAngularVelocity(p_idx) +
+ real_t(0.5) * wdot * dt_ );
}
+
ac.setOldForce(p_idx, ac.getForce(p_idx));
ac.setForce(p_idx, Vec3(real_t(0), real_t(0), real_t(0)));
+ ac.setOldTorque(p_idx, ac.getTorque(p_idx));
+ ac.setTorque(p_idx, Vec3(real_t(0), real_t(0), real_t(0)));
}
} //namespace kernel
diff --git a/src/mesa_pd/kernel/VelocityVerletWithShape.h b/src/mesa_pd/kernel/VelocityVerletWithShape.h
deleted file mode 100644
index 36c5b2f77..000000000
--- a/src/mesa_pd/kernel/VelocityVerletWithShape.h
+++ /dev/null
@@ -1,155 +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 VelocityVerletWithShape.h
-//! \author Sebastian Eibl <sebastian.eibl@fau.de>
-//
-//======================================================================================================================
-
-//======================================================================================================================
-//
-// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
-//
-//======================================================================================================================
-
-#pragma once
-
-#include <mesa_pd/data/DataTypes.h>
-#include <mesa_pd/data/IAccessor.h>
-
-namespace walberla {
-namespace mesa_pd {
-namespace kernel {
-
-/**
- * Velocity verlet integration for all particles.
- *
- * Velocit verlet integration is a two part kernel. preForceUpdate has to be
- * called before the force calculation and postFroceUpdate afterwards. The
- * integration is only complete when both functions are called. The integration
- * is symplectic.
- * \attention The force calculation has to be independent of velocity.
- *
- * This kernel requires the following particle accessor interface
- * \code
- * const walberla::mesa_pd::Vec3& getPosition(const size_t p_idx) const;
- * void setPosition(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
- *
- * const walberla::mesa_pd::Vec3& getLinearVelocity(const size_t p_idx) const;
- * void setLinearVelocity(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
- *
- * const walberla::real_t& getInvMass(const size_t p_idx) const;
- *
- * const walberla::mesa_pd::Vec3& getForce(const size_t p_idx) const;
- * void setForce(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
- *
- * const walberla::mesa_pd::Vec3& getOldForce(const size_t p_idx) const;
- * void setOldForce(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
- *
- * const walberla::mesa_pd::Rot3& getRotation(const size_t p_idx) const;
- * void setRotation(const size_t p_idx, const walberla::mesa_pd::Rot3& v);
- *
- * const walberla::mesa_pd::Vec3& getAngularVelocity(const size_t p_idx) const;
- * void setAngularVelocity(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
- *
- * const walberla::mesa_pd::Mat3& getInvInertiaBF(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);
- *
- * const walberla::mesa_pd::Vec3& getOldTorque(const size_t p_idx) const;
- * void setOldTorque(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
- *
- * const walberla::mesa_pd::data::particle_flags::FlagT& getFlags(const size_t p_idx) const;
- *
- * \endcode
- * \ingroup mesa_pd_kernel
- */
-class VelocityVerletWithShapePreForceUpdate
-{
-public:
- VelocityVerletWithShapePreForceUpdate(const real_t dt) : dt_(dt) {}
-
- template <typename Accessor>
- void operator()(const size_t i, Accessor& ac) const;
-
- real_t dt_;
-};
-
-/// \see VelocityVerletPreForceUpdate
-class VelocityVerletWithShapePostForceUpdate
-{
-public:
- VelocityVerletWithShapePostForceUpdate(const real_t dt) : dt_(dt) {}
-
- template <typename Accessor>
- void operator()(const size_t i, Accessor& ac) const;
-
- real_t dt_;
-};
-
-template <typename Accessor>
-inline void VelocityVerletWithShapePreForceUpdate::operator()(const size_t p_idx, Accessor& ac) const
-{
- static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
-
- if (!data::particle_flags::isSet( ac.getFlags(p_idx), data::particle_flags::FIXED))
- {
- 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);
-
- // Calculating the rotation angle
- const Vec3 phi( ac.getAngularVelocity(p_idx) * dt_ + real_t(0.5) * wdot * dt_ * dt_);
-
- // Calculating the new orientation
- auto rotation = ac.getRotation(p_idx);
- rotation.rotate( phi );
- ac.setRotation(p_idx, rotation);
- }
-}
-
-template <typename Accessor>
-inline void VelocityVerletWithShapePostForceUpdate::operator()(const size_t p_idx, Accessor& ac) const
-{
- static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
-
- if (!data::particle_flags::isSet( ac.getFlags(p_idx), data::particle_flags::FIXED))
- {
- 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;
-
- ac.setAngularVelocity(p_idx, ac.getAngularVelocity(p_idx) +
- real_t(0.5) * wdot * dt_ );
- }
-
- ac.setOldForce(p_idx, ac.getForce(p_idx));
- ac.setForce(p_idx, Vec3(real_t(0), real_t(0), real_t(0)));
-
- ac.setOldTorque(p_idx, ac.getTorque(p_idx));
- ac.setTorque(p_idx, Vec3(real_t(0), real_t(0), real_t(0)));
-}
-
-} //namespace kernel
-} //namespace mesa_pd
-} //namespace walberla
\ No newline at end of file
diff --git a/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp b/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
index 5aa9a51c2..2cc78a22f 100644
--- a/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
+++ b/tests/lbm_mesapd_coupling/momentum_exchange_method/SettlingSphere.cpp
@@ -71,10 +71,10 @@
#include "mesa_pd/data/shape/Sphere.h"
#include "mesa_pd/domain/BlockForestDomain.h"
#include "mesa_pd/kernel/DoubleCast.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/ParticleSelector.h"
#include "mesa_pd/kernel/SpringDashpot.h"
-#include "mesa_pd/kernel/VelocityVerletWithShape.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/mpi/ContactFilter.h"
#include "mesa_pd/mpi/SyncNextNeighbors.h"
#include "mesa_pd/mpi/ReduceProperty.h"
@@ -588,9 +588,9 @@ int main( int argc, char **argv )
syncCall();
- mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
- mesa_pd::kernel::VelocityVerletWithShapePreForceUpdate vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles));
- mesa_pd::kernel::VelocityVerletWithShapePostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles));
+ mesa_pd::kernel::ExplicitEuler explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
+ mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles));
+ mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles));
mesa_pd::kernel::SpringDashpot collisionResponse(1);
mesa_pd::mpi::ReduceProperty reduceProperty;
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index ffc728390..6397eb6f3 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -99,9 +99,6 @@ waLBerla_execute_test( NAME MESA_PD_Kernel_DoubleCast )
waLBerla_compile_test( NAME MESA_PD_Kernel_ExplicitEuler FILES kernel/ExplicitEuler.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_ExplicitEuler )
-waLBerla_compile_test( NAME MESA_PD_Kernel_ExplicitEulerWithShape FILES kernel/ExplicitEulerWithShape.cpp DEPENDS core )
-waLBerla_execute_test( NAME MESA_PD_Kernel_ExplicitEulerWithShape )
-
waLBerla_compile_test( NAME MESA_PD_Kernel_ForceLJ FILES kernel/ForceLJ.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_ForceLJ )
@@ -158,9 +155,6 @@ waLBerla_execute_test( NAME MESA_PD_Kernel_TemperatureIntegration )
waLBerla_compile_test( NAME MESA_PD_Kernel_VelocityVerlet FILES kernel/VelocityVerlet.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_VelocityVerlet )
-waLBerla_compile_test( NAME MESA_PD_Kernel_VelocityVerletWithShape FILES kernel/VelocityVerletWithShape.cpp DEPENDS core )
-waLBerla_execute_test( NAME MESA_PD_Kernel_VelocityVerletWithShape )
-
waLBerla_compile_test( NAME MESA_PD_MPI_BroadcastProperty FILES mpi/BroadcastProperty.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_MPI_BroadcastProperty PROCESSES 8 )
diff --git a/tests/mesa_pd/DropTestAnalytic.cpp b/tests/mesa_pd/DropTestAnalytic.cpp
index 6cdd51190..7c3c89ed0 100644
--- a/tests/mesa_pd/DropTestAnalytic.cpp
+++ b/tests/mesa_pd/DropTestAnalytic.cpp
@@ -25,7 +25,7 @@
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
+#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
@@ -102,7 +102,7 @@ int main( int argc, char ** argv )
real_t dt = real_t(0.00001);
// Init kernels
- kernel::ExplicitEulerWithShape explicitEulerWithShape( dt );
+ kernel::ExplicitEuler explicitEuler( dt );
kernel::SpringDashpot dem(1);
auto meff = real_t(1.0) / ss->shapes[sp->getShapeID()]->getInvMass();
dem.setParametersFromCOR(0,0,real_t(0.9), dt * real_t(20), meff);
@@ -152,7 +152,7 @@ int main( int argc, char ** argv )
ps->forEachParticle(false,
kernel::SelectLocal(),
accessor,
- explicitEulerWithShape,
+ explicitEuler,
accessor);
// if(i%1 == 0)
diff --git a/tests/mesa_pd/DropTestGeneral.cpp b/tests/mesa_pd/DropTestGeneral.cpp
index e729b1221..3d668694c 100644
--- a/tests/mesa_pd/DropTestGeneral.cpp
+++ b/tests/mesa_pd/DropTestGeneral.cpp
@@ -25,7 +25,7 @@
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
+#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
@@ -122,7 +122,7 @@ int main( int argc, char ** argv )
real_t dt = real_t(0.00001);
// Init kernels
- kernel::ExplicitEulerWithShape explicitEulerWithShape( dt );
+ kernel::ExplicitEuler explicitEuler( dt );
kernel::SpringDashpot dem(1);
auto meff = real_t(1.0) / ss->shapes[sp->getShapeID()]->getInvMass();
dem.setParametersFromCOR(0,0,real_t(0.9), dt * real_t(20), meff);
@@ -172,7 +172,7 @@ int main( int argc, char ** argv )
ps->forEachParticle(false,
kernel::SelectLocal(),
accessor,
- explicitEulerWithShape,
+ explicitEuler,
accessor);
// if(i%1 == 0)
diff --git a/tests/mesa_pd/IntegratorWithShapeEvaluation.ipynb b/tests/mesa_pd/IntegratorWithShapeEvaluation.ipynb
index 9427965a0..887e7d1ef 100644
--- a/tests/mesa_pd/IntegratorWithShapeEvaluation.ipynb
+++ b/tests/mesa_pd/IntegratorWithShapeEvaluation.ipynb
@@ -189,7 +189,7 @@
"metadata": {},
"outputs": [],
"source": [
- "data = np.loadtxt(\"VelocityVerletWithShape.txt\").transpose()\n",
+ "data = np.loadtxt(\"VelocityVerlet.txt\").transpose()\n",
"x = np.arange(len(data[0]))"
]
},
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
index 894cc47ee..7e316264f 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionLSD.cpp
@@ -26,7 +26,7 @@
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "core/Environment.h"
@@ -136,7 +136,7 @@ int main( int argc, char** argv )
collision_detection::AnalyticContactDetection acd;
kernel::DoubleCast double_cast;
- kernel::ExplicitEulerWithShape explEuler(dt);
+ kernel::ExplicitEuler explEuler(dt);
kernel::VelocityVerletPreForceUpdate vvPreForce( dt );
kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
index 0fa4a40cd..c422a70b4 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionNLSD.cpp
@@ -26,7 +26,7 @@
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/NonLinearSpringDashpot.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
@@ -137,7 +137,7 @@ int main( int argc, char** argv )
collision_detection::AnalyticContactDetection acd;
kernel::DoubleCast double_cast;
- kernel::ExplicitEulerWithShape explEuler(dt);
+ kernel::ExplicitEuler explEuler(dt);
kernel::VelocityVerletPreForceUpdate vvPreForce( dt );
kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
diff --git a/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp b/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
index ab3e92e67..300e0eb91 100644
--- a/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
+++ b/tests/mesa_pd/kernel/CoefficientOfRestitutionSD.cpp
@@ -26,7 +26,7 @@
#include "mesa_pd/kernel/DoubleCast.h"
#include "mesa_pd/kernel/VelocityVerlet.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/SpringDashpot.h"
#include "core/Environment.h"
@@ -136,7 +136,7 @@ int main( int argc, char** argv )
collision_detection::AnalyticContactDetection acd;
kernel::DoubleCast double_cast;
- kernel::ExplicitEulerWithShape explEuler(dt);
+ kernel::ExplicitEuler explEuler(dt);
kernel::VelocityVerletPreForceUpdate vvPreForce( dt );
kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
diff --git a/tests/mesa_pd/kernel/ExplicitEuler.cpp b/tests/mesa_pd/kernel/ExplicitEuler.cpp
index 34ac15885..1685ff1ae 100644
--- a/tests/mesa_pd/kernel/ExplicitEuler.cpp
+++ b/tests/mesa_pd/kernel/ExplicitEuler.cpp
@@ -75,14 +75,20 @@ int main( int argc, char ** argv )
integrator(0, accessor);
+ const auto& R = accessor.getRotation(0).getMatrix();
+ const auto wdot = R * accessor.getInvInertiaBF(0) * R.getTranspose() * torque;
+
//check force
WALBERLA_CHECK_FLOAT_EQUAL(accessor.getForce(0), Vec3(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getTorque(0), Vec3(0));
//check velocity
WALBERLA_CHECK_FLOAT_EQUAL(accessor.getLinearVelocity(0), force * accessor.getInvMass(0) * dt + linVel);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getAngularVelocity(0), wdot * dt + angVel);
//check position
WALBERLA_CHECK_FLOAT_EQUAL(accessor.getPosition(0), linVel * dt + force * accessor.getInvMass(0) * dt * dt);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getRotation(0).getQuaternion(), Quat( (wdot * dt + angVel).getNormalized(), (wdot * dt + angVel).length() * dt ));
accessor.setPosition( 0, Vec3(0,0,0));
accessor.setRotation( 0, Rot3(Quat()));
@@ -98,12 +104,15 @@ int main( int argc, char ** argv )
//check force
WALBERLA_CHECK_FLOAT_EQUAL(accessor.getForce(0), Vec3(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getTorque(0), Vec3(0));
//check velocity
WALBERLA_CHECK_FLOAT_EQUAL(accessor.getLinearVelocity(0), linVel);
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getAngularVelocity(0), angVel);
//check position
WALBERLA_CHECK_FLOAT_EQUAL(accessor.getPosition(0), Vec3(0,0,0));
+ WALBERLA_CHECK_FLOAT_EQUAL(accessor.getRotation(0).getQuaternion(), Quat());
return EXIT_SUCCESS;
}
diff --git a/tests/mesa_pd/kernel/ExplicitEulerWithShape.cpp b/tests/mesa_pd/kernel/ExplicitEulerWithShape.cpp
deleted file mode 100644
index b04e420b0..000000000
--- a/tests/mesa_pd/kernel/ExplicitEulerWithShape.cpp
+++ /dev/null
@@ -1,125 +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 ExplicitEulerWithShape.cpp
-//! \author Sebastian Eibl <sebastian.eibl@fau.de>
-//
-//======================================================================================================================
-
-#include <mesa_pd/data/DataTypes.h>
-#include <mesa_pd/data/ParticleAccessor.h>
-
-#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
-
-#include <core/Environment.h>
-#include <core/logging/Logging.h>
-
-#include <iostream>
-
-namespace walberla {
-
-using namespace walberla::mesa_pd;
-
-class SingleParticleAccessor : public data::SingleParticleAccessor
-{
-public:
- const walberla::real_t& getInvMass(const size_t /*p_idx*/) const {return invMass_;}
- 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::real_t invMass_;
- walberla::mesa_pd::Mat3 invInertiaBF_;
-};
-
-int main( int argc, char ** argv )
-{
- Environment env(argc, argv);
- WALBERLA_UNUSED(env);
- mpi::MPIManager::instance()->useWorldComm();
-
- //init data structures
- SingleParticleAccessor accessor;
-
- //initialize particle
- const auto linVel = Vec3(1,2,3);
- const auto angVel = Vec3(1,2,3);
-
- const auto force = Vec3(1,2,3);
- const auto torque = Vec3(1,2,3);
-
- accessor.setPosition( 0, Vec3(0,0,0));
- accessor.setRotation( 0, Rot3(Quat()));
- accessor.setLinearVelocity( 0, linVel);
- accessor.setAngularVelocity( 0, angVel);
- accessor.setForce( 0, force);
- accessor.setTorque( 0, torque);
- accessor.setInvMass( 0, real_t(1.23456));
- accessor.setInvInertiaBF( 0, Mat3(real_t(1.23456), real_t(0), real_t(0), real_t(0), real_t(1.23456), real_t(0), real_t(0), real_t(0), real_t(1.23456)));
-
- //init kernels
- const real_t dt = real_t(1);
- kernel::ExplicitEulerWithShape integrator( dt );
-
- integrator(0, accessor);
-
- const auto& R = accessor.getRotation(0).getMatrix();
- const auto wdot = R * accessor.getInvInertiaBF(0) * R.getTranspose() * torque;
-
- //check force
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getForce(0), Vec3(0));
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getTorque(0), Vec3(0));
-
- //check velocity
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getLinearVelocity(0), force * accessor.getInvMass(0) * dt + linVel);
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getAngularVelocity(0), wdot * dt + angVel);
-
- //check position
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getPosition(0), linVel * dt + force * accessor.getInvMass(0) * dt * dt);
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getRotation(0).getQuaternion(), Quat( (wdot * dt + angVel).getNormalized(), (wdot * dt + angVel).length() * dt ));
-
- accessor.setPosition( 0, Vec3(0,0,0));
- accessor.setRotation( 0, Rot3(Quat()));
- accessor.setLinearVelocity( 0, linVel);
- accessor.setAngularVelocity( 0, angVel);
- accessor.setForce( 0, force);
- accessor.setTorque( 0, torque);
- accessor.setInvMass( 0, real_t(1.23456));
- accessor.setInvInertiaBF( 0, Mat3(real_t(1.23456), real_t(0), real_t(0), real_t(0), real_t(1.23456), real_t(0), real_t(0), real_t(0), real_t(1.23456)));
- data::particle_flags::set( accessor.getFlagsRef(0), data::particle_flags::FIXED );
-
- integrator(0, accessor);
-
- //check force
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getForce(0), Vec3(0));
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getTorque(0), Vec3(0));
-
- //check velocity
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getLinearVelocity(0), linVel);
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getAngularVelocity(0), angVel);
-
- //check position
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getPosition(0), Vec3(0,0,0));
- WALBERLA_CHECK_FLOAT_EQUAL(accessor.getRotation(0).getQuaternion(), Quat());
-
- return EXIT_SUCCESS;
-}
-
-} //namespace walberla
-
-int main( int argc, char ** argv )
-{
- return walberla::main(argc, argv);
-}
diff --git a/tests/mesa_pd/kernel/GenerateLinkedCells.cpp b/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
index c54473cb0..e7a60f36e 100644
--- a/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
+++ b/tests/mesa_pd/kernel/GenerateLinkedCells.cpp
@@ -24,8 +24,9 @@
#include <mesa_pd/kernel/ParticleSelector.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 <core/Environment.h>
#include <core/grid_generator/SCIterator.h>
@@ -51,9 +52,12 @@ int main( int argc, char ** argv )
//init data structures
auto storage = std::make_shared<data::ParticleStorage>(100);
+ auto ss = std::make_shared<data::ShapeStorage>();
+ auto smallSphere = ss->create<data::Sphere>( real_t(1) );
+ ss->shapes[smallSphere]->updateMassAndInertia(real_t(2707));
data::LinkedCells linkedCells(domain.getScaled(2), real_t(1));
- data::ParticleAccessor accessor(storage);
+ data::ParticleAccessorWithShape accessor(storage, ss);
//initialize particles
for (auto it = grid_generator::SCIterator(domain, Vec3(spacing, spacing, spacing) * real_c(0.5), spacing);
@@ -62,6 +66,7 @@ int main( int argc, char ** argv )
{
data::Particle&& p = *storage->create();
p.getPositionRef() = (*it);
+ p.setShapeID(smallSphere);
}
//init kernels
diff --git a/tests/mesa_pd/kernel/IntegratorAccuracy.cpp b/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
index 099267523..695bae5b9 100644
--- a/tests/mesa_pd/kernel/IntegratorAccuracy.cpp
+++ b/tests/mesa_pd/kernel/IntegratorAccuracy.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/kernel/ParticleSelector.h"
@@ -105,7 +105,10 @@ int main( int argc, char ** argv )
//init data structures
auto ps = std::make_shared<data::ParticleStorage>(1);
- data::ParticleAccessor accessor(ps);
+ auto ss = std::make_shared<data::ShapeStorage>();
+ auto smallSphere = ss->create<data::Sphere>( real_t(1) );
+ ss->shapes[smallSphere]->updateMassAndInertia(real_t(2707));
+ data::ParticleAccessorWithShape accessor(ps, ss);
data::Particle&& p = *ps->create();
p.setPosition(pos);
@@ -113,6 +116,7 @@ int main( int argc, char ** argv )
p.setForce(getForce(pos, k));
p.setOldForce(getForce(pos, k));
p.setInvMass(real_t(1) / mass);
+ p.setShapeID(smallSphere);
// velocity verlet
kernel::VelocityVerletPreForceUpdate preForce( dt );
diff --git a/tests/mesa_pd/kernel/LinearSpringDashpot.cpp b/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
index 0513e8086..18beb0732 100644
--- a/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
+++ b/tests/mesa_pd/kernel/LinearSpringDashpot.cpp
@@ -26,8 +26,8 @@
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/kernel/DoubleCast.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
-#include "mesa_pd/kernel/VelocityVerletWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
#include "mesa_pd/kernel/LinearSpringDashpot.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
@@ -148,11 +148,11 @@ int main( int argc, char ** argv )
data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::FIXED);
// velocity verlet
- kernel::VelocityVerletWithShapePreForceUpdate vvPreForce( dt );
- kernel::VelocityVerletWithShapePostForceUpdate vvPostForce( dt );
+ kernel::VelocityVerletPreForceUpdate vvPreForce( dt );
+ kernel::VelocityVerletPostForceUpdate vvPostForce( dt );
// explicit euler
- kernel::ExplicitEulerWithShape explEuler( dt );
+ kernel::ExplicitEuler explEuler( dt );
// collision response
collision_detection::AnalyticContactDetection acd;
diff --git a/tests/mesa_pd/kernel/SpherePile.cpp b/tests/mesa_pd/kernel/SpherePile.cpp
index e326a1110..67e1087a7 100644
--- a/tests/mesa_pd/kernel/SpherePile.cpp
+++ b/tests/mesa_pd/kernel/SpherePile.cpp
@@ -26,7 +26,7 @@
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/kernel/DoubleCast.h"
-#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
#include "mesa_pd/kernel/SpringDashpot.h"
#include "mesa_pd/kernel/SpringDashpotSpring.h"
#include "mesa_pd/mpi/ReduceContactHistory.h"
@@ -98,7 +98,7 @@ int main(int argc, char **argv)
false);
// explicit euler
- kernel::ExplicitEulerWithShape explEuler(dt);
+ kernel::ExplicitEuler explEuler(dt);
collision_detection::AnalyticContactDetection acd;
kernel::DoubleCast double_cast;
kernel::SpringDashpot sd(1);
diff --git a/tests/mesa_pd/kernel/VelocityVerlet.cpp b/tests/mesa_pd/kernel/VelocityVerlet.cpp
index e322d78c1..1dba566fb 100644
--- a/tests/mesa_pd/kernel/VelocityVerlet.cpp
+++ b/tests/mesa_pd/kernel/VelocityVerlet.cpp
@@ -20,6 +20,7 @@
#include <mesa_pd/data/ParticleAccessor.h>
#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/data/shape/Sphere.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/VelocityVerlet.h>
@@ -27,43 +28,68 @@
#include <core/Environment.h>
#include <core/logging/Logging.h>
+#include <cmath>
+#include <fstream>
#include <iostream>
namespace walberla {
using namespace walberla::mesa_pd;
-int main( int argc, char ** argv )
+class ParticleAccessorWithShape : public data::ParticleAccessor
{
- Environment env(argc, argv);
- WALBERLA_UNUSED(env);
- mpi::MPIManager::instance()->useWorldComm();
+public:
+ ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps)
+ : ParticleAccessor(ps)
+ {
+ sp.updateMassAndInertia(real_t(1234));
+ }
- const real_t k = real_t(0.1);
+ const auto& getInvMass(const size_t /*p_idx*/) const {return sp.getInvMass();}
+ const auto& getInvInertiaBF(const size_t /*p_idx*/) const {return sp.getInvInertiaBF();}
+
+ data::BaseShape* getShape(const size_t /*p_idx*/) {return &sp;}
+private:
+ data::Sphere sp = data::Sphere(real_t(0.6));
+};
+
+int main()
+{
//init data structures
auto ps = std::make_shared<data::ParticleStorage>(100);
- data::ParticleAccessor accessor(ps);
+ ParticleAccessorWithShape accessor(ps);
const Vec3 startingVelocity(real_t(1),real_t(2),real_t(3));
//initialize particle
data::Particle&& p = *ps->create();
p.getPositionRef() = Vec3(0,0,0);
- p.getInvMassRef() = real_t(1.1);
p.getLinearVelocityRef() = startingVelocity;
+ p.getAngularVelocityRef() = startingVelocity;
p.getForceRef() = Vec3(0,0,0);
p.getOldForceRef() = Vec3(0,0,0);
- const real_t dt = real_t(0.1);
+ const real_t dt = real_t(0.001);
+
+ const real_t k = real_t(10000);
const real_t w = std::sqrt(k*accessor.getInvMass(0));
- const Vec3 A = accessor.getLinearVelocity(0) / w;
+ const Vec3 A = accessor.getLinearVelocity(0) / w;
WALBERLA_LOG_DEVEL_VAR(w);
WALBERLA_LOG_DEVEL_VAR(A);
auto analyticPosition = [A, w](const real_t t){return A * std::sin(w*t);};
auto analyticVelocity = [A, w](const real_t t){return A * std::cos(w*t) * w;};
+ const real_t kappa = real_t(1000);
+ const real_t omega = std::sqrt(kappa*accessor.getInvInertiaBF(0)[0]);
+ const Vec3 Alpha = accessor.getAngularVelocity(0) / omega;
+ WALBERLA_LOG_DEVEL_VAR(omega);
+ WALBERLA_LOG_DEVEL_VAR(Alpha);
+
+ auto analyticRotation = [Alpha, omega](const real_t t){return Alpha * std::sin(omega*t);};
+ auto analyticAngVel = [Alpha, omega](const real_t t){return Alpha * std::cos(omega*t) * omega;};
+
//init kernels
kernel::VelocityVerletPreForceUpdate preForce( dt );
kernel::VelocityVerletPostForceUpdate postForce( dt );
@@ -71,41 +97,84 @@ int main( int argc, char ** argv )
p.getPositionRef() = analyticPosition(-dt);
p.getLinearVelocityRef() = analyticVelocity(-dt);
p.getOldForceRef() = - k * p.getPosition();
+
+ p.getRotationRef() = Rot3(Quat(analyticRotation(-dt), analyticRotation(-dt).length()));
+ p.getAngularVelocityRef() = analyticAngVel(-dt);
+ p.getOldTorqueRef() = - kappa * p.getRotation().getQuaternion().getAxis() * p.getRotation().getQuaternion().getAngle();
+
ps->forEachParticle(false, kernel::SelectAll(), accessor, preForce, accessor);
+
p.getLinearVelocityRef() = analyticVelocity(real_t(0));
p.getForceRef() = - k * p.getPosition();
+
+ p.getAngularVelocityRef() = analyticAngVel(real_t(0));
+ p.getTorqueRef() = - kappa * p.getRotation().getQuaternion().getAxis() * p.getRotation().getQuaternion().getAngle();
ps->forEachParticle(false, kernel::SelectAll(), accessor, postForce, accessor);
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getPosition(), Vec3(0), real_t(1e-2));
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getLinearVelocity(), startingVelocity, real_t(1e-2));
- for (auto i = 1; i < 500; ++i)
+ std::fstream fout("VelocityVerlet.txt", std::ofstream::out);
+ for (auto i = 1; i < 3000; ++i)
{
ps->forEachParticle(false, kernel::SelectAll(), accessor, preForce, accessor);
p.getForceRef() = - k * p.getPosition();
+ p.getTorqueRef() = - kappa * p.getRotation().getQuaternion().getAxis() * p.getRotation().getQuaternion().getAngle();
ps->forEachParticle(false, kernel::SelectAll(), accessor, postForce, accessor);
- //check force
+ //check force&torque
WALBERLA_CHECK_FLOAT_EQUAL(p.getForce(), Vec3(0), p);
+ WALBERLA_CHECK_FLOAT_EQUAL(p.getTorque(), Vec3(0), p);
//check velocity
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getLinearVelocity(),
analyticVelocity(real_c(i) * dt),
- real_t(1e-2),
+ real_t(1e-4),
"iteration: " << i << "\n" <<
"t: " << real_c(i)*dt << "\n" <<
p);
+ //check angular velocity
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getAngularVelocity(),
+ analyticAngVel(real_c(i) * dt),
+ real_t(1e-4),
+ "iteration: " << i << "\n" <<
+ "t: " << real_c(i)*dt << "\n" <<
+ p);
+ // WALBERLA_LOG_DEVEL_VAR(p.getAngularVelocity());
+ // WALBERLA_LOG_DEVEL_VAR(analyticAngVel(real_c(i) * dt));
+
//check position
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getPosition(),
analyticPosition(real_c(i) * dt),
- real_t(1e-2),
+ real_t(1e-4),
"iteration: " << i << "\n" <<
"t: " << real_c(i)*dt << "\n" <<
p);
-// WALBERLA_LOG_DEVEL_VAR(p.getPosition());
-// WALBERLA_LOG_DEVEL_VAR(analyticPosition(real_c(i) * dt));
+
+ //check rotation
+ auto angSol = Rot3(Quat(analyticRotation(real_c(i) * dt), analyticRotation(real_c(i) * dt).length()));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getRotation().getQuaternion().getAxis() * p.getRotation().getQuaternion().getAngle(),
+ angSol.getQuaternion().getAxis() * angSol.getQuaternion().getAngle(),
+ real_t(1e-4),
+ "iteration: " << i << "\n" <<
+ "t: " << real_c(i)*dt << "\n" <<
+ p);
+ // WALBERLA_LOG_DEVEL_VAR(p.getPosition());
+ // WALBERLA_LOG_DEVEL_VAR(analyticPosition(real_c(i) * dt));
+
+ // WALBERLA_LOG_DEVEL_VAR(/*p.getRotation().getQuaternion().getAxis() * */p.getRotation().getQuaternion().getAngle());
+ // WALBERLA_LOG_DEVEL_VAR(/*angSol.getQuaternion().getAxis() * */angSol.getQuaternion().getAngle());
+ fout << p.getPosition().length() << " " <<
+ analyticPosition(real_c(i) * dt).length() << " " <<
+ p.getLinearVelocity().length() << " " <<
+ analyticVelocity(real_c(i) * dt).length() << " " <<
+ p.getRotation().getQuaternion().getAngle() << " " <<
+ angSol.getQuaternion().getAngle() << " " <<
+ p.getAngularVelocity().length() << " " <<
+ analyticAngVel(real_c(i) * dt).length() << std::endl;
}
+ fout.close();
return EXIT_SUCCESS;
}
@@ -114,5 +183,15 @@ int main( int argc, char ** argv )
int main( int argc, char ** argv )
{
- return walberla::main(argc, argv);
+ walberla::Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+ walberla::mpi::MPIManager::instance()->useWorldComm();
+
+ if (std::is_same<walberla::real_t, float>::value)
+ {
+ WALBERLA_LOG_WARNING("waLBerla build in sp mode: skipping test due to low precision");
+ return EXIT_SUCCESS;
+ }
+
+ return walberla::main();
}
diff --git a/tests/mesa_pd/kernel/VelocityVerletWithShape.cpp b/tests/mesa_pd/kernel/VelocityVerletWithShape.cpp
deleted file mode 100644
index 70d778067..000000000
--- a/tests/mesa_pd/kernel/VelocityVerletWithShape.cpp
+++ /dev/null
@@ -1,197 +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 VelocityVerletWithShape.cpp
-//! \author Sebastian Eibl <sebastian.eibl@fau.de>
-//
-//======================================================================================================================
-
-#include <mesa_pd/data/ParticleAccessor.h>
-#include <mesa_pd/data/ParticleStorage.h>
-#include <mesa_pd/data/shape/Sphere.h>
-
-#include <mesa_pd/kernel/ParticleSelector.h>
-#include <mesa_pd/kernel/VelocityVerletWithShape.h>
-
-#include <core/Environment.h>
-#include <core/logging/Logging.h>
-
-#include <cmath>
-#include <fstream>
-#include <iostream>
-
-namespace walberla {
-
-using namespace walberla::mesa_pd;
-
-class ParticleAccessorWithShape : public data::ParticleAccessor
-{
-public:
- ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps)
- : ParticleAccessor(ps)
- {
- sp.updateMassAndInertia(real_t(1234));
- }
-
- const auto& getInvMass(const size_t /*p_idx*/) const {return sp.getInvMass();}
-
- const auto& getInvInertiaBF(const size_t /*p_idx*/) const {return sp.getInvInertiaBF();}
-
- data::BaseShape* getShape(const size_t /*p_idx*/) {return &sp;}
-private:
- data::Sphere sp = data::Sphere(real_t(0.6));
-};
-
-int main()
-{
- //init data structures
- auto ps = std::make_shared<data::ParticleStorage>(100);
- ParticleAccessorWithShape accessor(ps);
-
- const Vec3 startingVelocity(real_t(1),real_t(2),real_t(3));
-
- //initialize particle
- data::Particle&& p = *ps->create();
- p.getPositionRef() = Vec3(0,0,0);
- p.getLinearVelocityRef() = startingVelocity;
- p.getAngularVelocityRef() = startingVelocity;
- p.getForceRef() = Vec3(0,0,0);
- p.getOldForceRef() = Vec3(0,0,0);
-
- const real_t dt = real_t(0.001);
-
- const real_t k = real_t(10000);
- const real_t w = std::sqrt(k*accessor.getInvMass(0));
- const Vec3 A = accessor.getLinearVelocity(0) / w;
- WALBERLA_LOG_DEVEL_VAR(w);
- WALBERLA_LOG_DEVEL_VAR(A);
-
- auto analyticPosition = [A, w](const real_t t){return A * std::sin(w*t);};
- auto analyticVelocity = [A, w](const real_t t){return A * std::cos(w*t) * w;};
-
- const real_t kappa = real_t(1000);
- const real_t omega = std::sqrt(kappa*accessor.getInvInertiaBF(0)[0]);
- const Vec3 Alpha = accessor.getAngularVelocity(0) / omega;
- WALBERLA_LOG_DEVEL_VAR(omega);
- WALBERLA_LOG_DEVEL_VAR(Alpha);
-
- auto analyticRotation = [Alpha, omega](const real_t t){return Alpha * std::sin(omega*t);};
- auto analyticAngVel = [Alpha, omega](const real_t t){return Alpha * std::cos(omega*t) * omega;};
-
- //init kernels
- kernel::VelocityVerletWithShapePreForceUpdate preForce( dt );
- kernel::VelocityVerletWithShapePostForceUpdate postForce( dt );
-
- p.getPositionRef() = analyticPosition(-dt);
- p.getLinearVelocityRef() = analyticVelocity(-dt);
- p.getOldForceRef() = - k * p.getPosition();
-
- p.getRotationRef() = Rot3(Quat(analyticRotation(-dt), analyticRotation(-dt).length()));
- p.getAngularVelocityRef() = analyticAngVel(-dt);
- p.getOldTorqueRef() = - kappa * p.getRotation().getQuaternion().getAxis() * p.getRotation().getQuaternion().getAngle();
-
- ps->forEachParticle(false, kernel::SelectAll(), accessor, preForce, accessor);
-
- p.getLinearVelocityRef() = analyticVelocity(real_t(0));
- p.getForceRef() = - k * p.getPosition();
-
- p.getAngularVelocityRef() = analyticAngVel(real_t(0));
- p.getTorqueRef() = - kappa * p.getRotation().getQuaternion().getAxis() * p.getRotation().getQuaternion().getAngle();
- ps->forEachParticle(false, kernel::SelectAll(), accessor, postForce, accessor);
-
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getPosition(), Vec3(0), real_t(1e-2));
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getLinearVelocity(), startingVelocity, real_t(1e-2));
-
- std::fstream fout("VelocityVerletWithShape.txt", std::ofstream::out);
- for (auto i = 1; i < 3000; ++i)
- {
- ps->forEachParticle(false, kernel::SelectAll(), accessor, preForce, accessor);
- p.getForceRef() = - k * p.getPosition();
- p.getTorqueRef() = - kappa * p.getRotation().getQuaternion().getAxis() * p.getRotation().getQuaternion().getAngle();
- ps->forEachParticle(false, kernel::SelectAll(), accessor, postForce, accessor);
-
- //check force&torque
- WALBERLA_CHECK_FLOAT_EQUAL(p.getForce(), Vec3(0), p);
- WALBERLA_CHECK_FLOAT_EQUAL(p.getTorque(), Vec3(0), p);
-
- //check velocity
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getLinearVelocity(),
- analyticVelocity(real_c(i) * dt),
- real_t(1e-4),
- "iteration: " << i << "\n" <<
- "t: " << real_c(i)*dt << "\n" <<
- p);
-
- //check angular velocity
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getAngularVelocity(),
- analyticAngVel(real_c(i) * dt),
- real_t(1e-4),
- "iteration: " << i << "\n" <<
- "t: " << real_c(i)*dt << "\n" <<
- p);
- // WALBERLA_LOG_DEVEL_VAR(p.getAngularVelocity());
- // WALBERLA_LOG_DEVEL_VAR(analyticAngVel(real_c(i) * dt));
-
- //check position
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getPosition(),
- analyticPosition(real_c(i) * dt),
- real_t(1e-4),
- "iteration: " << i << "\n" <<
- "t: " << real_c(i)*dt << "\n" <<
- p);
-
- //check rotation
- auto angSol = Rot3(Quat(analyticRotation(real_c(i) * dt), analyticRotation(real_c(i) * dt).length()));
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(p.getRotation().getQuaternion().getAxis() * p.getRotation().getQuaternion().getAngle(),
- angSol.getQuaternion().getAxis() * angSol.getQuaternion().getAngle(),
- real_t(1e-4),
- "iteration: " << i << "\n" <<
- "t: " << real_c(i)*dt << "\n" <<
- p);
- // WALBERLA_LOG_DEVEL_VAR(p.getPosition());
- // WALBERLA_LOG_DEVEL_VAR(analyticPosition(real_c(i) * dt));
-
- // WALBERLA_LOG_DEVEL_VAR(/*p.getRotation().getQuaternion().getAxis() * */p.getRotation().getQuaternion().getAngle());
- // WALBERLA_LOG_DEVEL_VAR(/*angSol.getQuaternion().getAxis() * */angSol.getQuaternion().getAngle());
- fout << p.getPosition().length() << " " <<
- analyticPosition(real_c(i) * dt).length() << " " <<
- p.getLinearVelocity().length() << " " <<
- analyticVelocity(real_c(i) * dt).length() << " " <<
- p.getRotation().getQuaternion().getAngle() << " " <<
- angSol.getQuaternion().getAngle() << " " <<
- p.getAngularVelocity().length() << " " <<
- analyticAngVel(real_c(i) * dt).length() << std::endl;
- }
- fout.close();
-
- return EXIT_SUCCESS;
-}
-
-} //namespace walberla
-
-int main( int argc, char ** argv )
-{
- walberla::Environment env(argc, argv);
- WALBERLA_UNUSED(env);
- walberla::mpi::MPIManager::instance()->useWorldComm();
-
- if (std::is_same<walberla::real_t, float>::value)
- {
- WALBERLA_LOG_WARNING("waLBerla build in sp mode: skipping test due to low precision");
- return EXIT_SUCCESS;
- }
-
- return walberla::main();
-}
diff --git a/tests/mesa_pd/kernel/interfaces/ExplicitEulerInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/ExplicitEulerInterfaceCheck.cpp
index 792965081..4d180c321 100644
--- a/tests/mesa_pd/kernel/interfaces/ExplicitEulerInterfaceCheck.cpp
+++ b/tests/mesa_pd/kernel/interfaces/ExplicitEulerInterfaceCheck.cpp
@@ -51,6 +51,17 @@ public:
const walberla::mesa_pd::data::particle_flags::FlagT& getFlags(const size_t /*p_idx*/) const {return flags_;}
+ 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;}
+
id_t getInvalidUid() const {return UniqueID<int>::invalidID();}
size_t getInvalidIdx() const {return std::numeric_limits<size_t>::max();}
@@ -67,6 +78,10 @@ private:
walberla::real_t invMass_;
walberla::mesa_pd::Vec3 force_;
walberla::mesa_pd::data::particle_flags::FlagT flags_;
+ walberla::mesa_pd::Rot3 rotation_;
+ walberla::mesa_pd::Vec3 angularVelocity_;
+ walberla::mesa_pd::Mat3 invInertiaBF_;
+ walberla::mesa_pd::Vec3 torque_;
};
template void kernel::ExplicitEuler::operator()(const size_t p_idx1, Accessor& ac) const;
diff --git a/tests/mesa_pd/kernel/interfaces/ExplicitEulerWithShapeInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/ExplicitEulerWithShapeInterfaceCheck.cpp
index 3737bd292..a56e29914 100644
--- a/tests/mesa_pd/kernel/interfaces/ExplicitEulerWithShapeInterfaceCheck.cpp
+++ b/tests/mesa_pd/kernel/interfaces/ExplicitEulerWithShapeInterfaceCheck.cpp
@@ -13,7 +13,7 @@
// 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 ExplicitEulerWithShapeInterfaceCheck.cpp
+//! \file ExplicitEulerInterfaceCheck.cpp
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//
//======================================================================================================================
@@ -25,7 +25,7 @@
//======================================================================================================================
#include <mesa_pd/data/IAccessor.h>
-#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
+#include <mesa_pd/kernel/ExplicitEuler.h>
#include <core/UniqueID.h>
@@ -84,7 +84,7 @@ private:
walberla::mesa_pd::data::particle_flags::FlagT flags_;
};
-template void kernel::ExplicitEulerWithShape::operator()(const size_t p_idx1, Accessor& ac) const;
+template void kernel::ExplicitEuler::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/kernel/interfaces/VelocityVerletInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/VelocityVerletInterfaceCheck.cpp
index e0cca71cc..06d8507c7 100644
--- a/tests/mesa_pd/kernel/interfaces/VelocityVerletInterfaceCheck.cpp
+++ b/tests/mesa_pd/kernel/interfaces/VelocityVerletInterfaceCheck.cpp
@@ -52,6 +52,20 @@ public:
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_;}
@@ -70,6 +84,11 @@ private:
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_;
};
diff --git a/tests/mesa_pd/kernel/interfaces/VelocityVerletWithShapeInterfaceCheck.cpp b/tests/mesa_pd/kernel/interfaces/VelocityVerletWithShapeInterfaceCheck.cpp
index d000cb66e..06d8507c7 100644
--- a/tests/mesa_pd/kernel/interfaces/VelocityVerletWithShapeInterfaceCheck.cpp
+++ b/tests/mesa_pd/kernel/interfaces/VelocityVerletWithShapeInterfaceCheck.cpp
@@ -13,7 +13,7 @@
// 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 VelocityVerletWithShapeInterfaceCheck.cpp
+//! \file VelocityVerletInterfaceCheck.cpp
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//
//======================================================================================================================
@@ -25,7 +25,7 @@
//======================================================================================================================
#include <mesa_pd/data/IAccessor.h>
-#include <mesa_pd/kernel/VelocityVerletWithShape.h>
+#include <mesa_pd/kernel/VelocityVerlet.h>
#include <core/UniqueID.h>
@@ -92,8 +92,8 @@ private:
walberla::mesa_pd::data::particle_flags::FlagT flags_;
};
-template void kernel::VelocityVerletWithShapePreForceUpdate::operator()(const size_t p_idx1, Accessor& ac) const;
-template void kernel::VelocityVerletWithShapePostForceUpdate::operator()(const size_t p_idx1, Accessor& ac) const;
+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
--
GitLab