Skip to content
GitLab
Commit b67091ed authored by Christoph Rettinger's avatar Christoph Rettinger Committed by Christoph Schwarzmeier
Browse files

Fix of integration of angular velocity of non-spherical particles

parent 27e2bc06
Hide whitespace changes
Inline Side-by-side
apps/benchmarks/GranularGas/MESA_PD_GranularGas.cpp
View file @ b67091ed
......@@ -18,7 +18,6 @@
//
//======================================================================================================================
#include "Accessor.h"
#include "check.h"
#include "Contact.h"
#include "CreateParticles.h"
......@@ -32,7 +31,7 @@
#include <mesa_pd/collision_detection/AnalyticContactDetection.h>
#include <mesa_pd/data/LinkedCells.h>
#include <mesa_pd/data/ParticleAccessor.h>
#include <mesa_pd/data/ParticleAccessorWithShape.h>
#include <mesa_pd/data/ParticleStorage.h>
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/domain/BlockForestDomain.h>
......@@ -90,7 +89,7 @@ public:
}
inline
void operator()(const size_t idx1, const size_t idx2, ParticleAccessorWithShape& ac)
void operator()(const size_t idx1, const size_t idx2, data::ParticleAccessorWithShape& ac)
{
++contactsChecked_;
......@@ -194,7 +193,7 @@ int main( int argc, char ** argv )
//init data structures
auto ps = std::make_shared<data::ParticleStorage>(100);
auto ss = std::make_shared<data::ShapeStorage>();
ParticleAccessorWithShape accessor(ps, ss);
data::ParticleAccessorWithShape accessor(ps, ss);
data::LinkedCells lc(localDomain.getExtended(params.spacing), params.spacing );
auto smallSphere = ss->create<data::Sphere>( params.radius );
......
apps/benchmarks/GranularGas/MESA_PD_KernelBenchmark.cpp
View file @ b67091ed
......@@ -18,7 +18,6 @@
//
//======================================================================================================================
#include "Accessor.h"
#include "check.h"
#include "Contact.h"
#include "CreateParticles.h"
......@@ -32,7 +31,7 @@
#include <mesa_pd/collision_detection/AnalyticContactDetection.h>
#include <mesa_pd/data/LinkedCells.h>
#include <mesa_pd/data/ParticleAccessor.h>
#include <mesa_pd/data/ParticleAccessorWithShape.h>
#include <mesa_pd/data/ParticleStorage.h>
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/domain/BlockForestDomain.h>
......@@ -209,7 +208,7 @@ int main( int argc, char ** argv )
//init data structures
auto ss = std::make_shared<data::ShapeStorage>();
auto ps = std::make_shared<data::ParticleStorage>(100);
ParticleAccessorWithShape accessor(ps, ss);
data::ParticleAccessorWithShape accessor(ps, ss);
data::LinkedCells lc(localDomain.getExtended(params.spacing), params.spacing );
auto smallSphere = ss->create<data::Sphere>( params.radius );
......
apps/benchmarks/GranularGas/MESA_PD_KernelLoadBalancing.cpp
View file @ b67091ed
......@@ -18,7 +18,6 @@
//
//======================================================================================================================
#include "Accessor.h"
#include "check.h"
#include "Contact.h"
#include "CreateParticles.h"
......@@ -32,7 +31,7 @@
#include <mesa_pd/collision_detection/AnalyticContactDetection.h>
#include <mesa_pd/data/LinkedCells.h>
#include <mesa_pd/data/ParticleAccessor.h>
#include <mesa_pd/data/ParticleAccessorWithShape.h>
#include <mesa_pd/data/ParticleStorage.h>
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/data/SparseLinkedCells.h>
......@@ -198,7 +197,7 @@ int main( int argc, char ** argv )
//init data structures
auto ps = std::make_shared<data::ParticleStorage>(100);
auto ss = std::make_shared<data::ShapeStorage>();
ParticleAccessorWithShape accessor(ps, ss);
data::ParticleAccessorWithShape accessor(ps, ss);
auto lc = std::make_shared<data::LinkedCells>(domain->getUnionOfLocalAABBs().getExtended(params.spacing), params.spacing );
forest->addBlockData(domain::createBlockForestDataHandling(ps), "Storage");
......
apps/benchmarks/GranularGas/MESA_PD_LoadBalancing.cpp
View file @ b67091ed
......@@ -18,7 +18,6 @@
//
//======================================================================================================================
#include "Accessor.h"
#include "check.h"
#include "Contact.h"
#include "CreateParticles.h"
......@@ -32,7 +31,7 @@
#include <mesa_pd/collision_detection/AnalyticContactDetection.h>
#include <mesa_pd/data/SparseLinkedCells.h>
#include <mesa_pd/data/ParticleAccessor.h>
#include <mesa_pd/data/ParticleAccessorWithShape.h>
#include <mesa_pd/data/ParticleStorage.h>
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/data/STLOverloads.h>
......@@ -236,7 +235,7 @@ int main( int argc, char ** argv )
//init data structures
auto ps = std::make_shared<data::ParticleStorage>(100);
auto ss = std::make_shared<data::ShapeStorage>();
ParticleAccessorWithShape accessor(ps, ss);
data::ParticleAccessorWithShape accessor(ps, ss);
auto lc = std::make_shared<data::SparseLinkedCells>(domain->getUnionOfLocalAABBs().getExtended(params.spacing), params.spacing );
forest->addBlockData(domain::createBlockForestDataHandling(ps), "Storage");
......
apps/benchmarks/IntegratorAccuracy/IntegratorAccuracy.cpp
View file @ b67091ed
......@@ -43,15 +43,15 @@ public:
const auto &getInvMass(const size_t /*p_idx*/) const
{ return invMass_; }
void setInvInertiaBF(const size_t /*p_idx*/, const Mat3 &val)
{ invInertiaBF_ = val; }
const auto &getInvInertiaBF(const size_t /*p_idx*/) const // dummy
{ return dummyI_; }
const auto &getInvInertiaBF(const size_t /*p_idx*/) const
{ return invInertiaBF_; }
const auto &getInertiaBF(const size_t /*p_idx*/) const // dummy
{ return dummyI_; }
private:
real_t invMass_;
Mat3 invInertiaBF_;
Mat3 dummyI_{0_r};
};
struct Oscillator
......
apps/showcases/ParticlePacking/Evaluation.h
View file @ b67091ed
......@@ -20,6 +20,7 @@
#pragma once
#include "mesa_pd/data/ContactAccessor.h"
#include "mesa_pd/data/ParticleStorage.h"
#include "Utility.h"
......
apps/showcases/ParticlePacking/ParticlePacking.cpp
View file @ b67091ed
......@@ -35,6 +35,7 @@
#include "mesa_pd/data/ContactStorage.h"
#include "mesa_pd/data/ContactAccessor.h"
#include "mesa_pd/data/ParticleAccessorWithBaseShape.h"
#include "mesa_pd/data/ParticleStorage.h"
#include "mesa_pd/data/HashGrids.h"
......@@ -418,7 +419,7 @@ int main(int argc, char **argv) {
/// MESAPD Data
auto particleStorage = std::make_shared<data::ParticleStorage>(1);
auto contactStorage = std::make_shared<data::ContactStorage>(1);
ParticleAccessorWithShape particleAccessor(particleStorage);
data::ParticleAccessorWithBaseShape particleAccessor(particleStorage);
data::ContactAccessor contactAccessor(contactStorage);
// configure shape creation
......@@ -626,7 +627,7 @@ int main(int argc, char **argv) {
math::DistributedSample diameterSample;
particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
[&diameterSample](const size_t idx, ParticleAccessorWithShape& ac){diameterSample.insert(real_c(2)*ac.getInteractionRadius(idx));}, particleAccessor);
[&diameterSample](const size_t idx, data::ParticleAccessorWithBaseShape& ac){diameterSample.insert(real_c(2)*ac.getInteractionRadius(idx));}, particleAccessor);
diameterSample.mpiAllGather();
WALBERLA_LOG_INFO_ON_ROOT("Statistics of initially created particles' interaction diameters: " << diameterSample.format());
......@@ -713,7 +714,7 @@ int main(int argc, char **argv) {
meshParticleVTK.addFaceOutput< data::SelectParticleLinearVelocity >("LinearVelocity");
meshParticleVTK.addVertexOutput< data::SelectParticlePosition >("Position");
meshParticleVTK.addVertexOutput< data::SelectParticleNumContacts >("numContacts");
auto surfaceVelDataSource = make_shared<mesa_pd::SurfaceVelocityVertexDataSource< mesh::PolyMesh, ParticleAccessorWithShape > >("SurfaceVelocity", particleAccessor);
auto surfaceVelDataSource = make_shared<mesa_pd::SurfaceVelocityVertexDataSource< mesh::PolyMesh, data::ParticleAccessorWithBaseShape > >("SurfaceVelocity", particleAccessor);
meshParticleVTK.setParticleSelector(vtkParticleSelector);
meshParticleVTK.addVertexDataSource(surfaceVelDataSource);
......@@ -850,7 +851,7 @@ int main(int argc, char **argv) {
timing.start("Contact detection");
hashGrids.forEachParticlePairHalf(useOpenMP, kernel::ExcludeInfiniteInfinite(), particleAccessor,
[domain, contactStorage](size_t idx1, size_t idx2, ParticleAccessorWithShape &ac){
[domain, contactStorage](size_t idx1, size_t idx2, data::ParticleAccessorWithBaseShape &ac){
kernel::DoubleCast double_cast;
mpi::ContactFilter contact_filter;
collision_detection::GeneralContactDetection contactDetection;
......@@ -890,7 +891,7 @@ int main(int argc, char **argv) {
} else
{
linkedCells.forEachParticlePairHalf(useOpenMP, kernel::ExcludeInfiniteInfinite(), particleAccessor,
[domain, contactStorage](size_t idx1, size_t idx2, ParticleAccessorWithShape &ac){
[domain, contactStorage](size_t idx1, size_t idx2, data::ParticleAccessorWithBaseShape &ac){
collision_detection::GeneralContactDetection contactDetection;
//Attention: does not use contact threshold in general case (GJK)
......@@ -926,9 +927,9 @@ int main(int argc, char **argv) {
timing.start("Contact eval");
particleStorage->forEachParticle(useOpenMP, kernel::SelectAll(), particleAccessor,
[](size_t p_idx, ParticleAccessorWithShape& ac){ac.setNumContacts(p_idx,0);}, particleAccessor);
[](size_t p_idx, data::ParticleAccessorWithBaseShape& ac){ac.setNumContacts(p_idx,0);}, particleAccessor);
contactStorage->forEachContact(useOpenMP, kernel::SelectAll(), contactAccessor,
[](size_t c, data::ContactAccessor &ca, ParticleAccessorWithShape &pa) {
[](size_t c, data::ContactAccessor &ca, data::ParticleAccessorWithBaseShape &pa) {
auto idx1 = ca.getId1(c);
auto idx2 = ca.getId2(c);
pa.getNumContactsRef(idx1)++;
......@@ -996,7 +997,7 @@ int main(int argc, char **argv) {
timing.start("DEM");
timing.start("Collision");
contactStorage->forEachContact(useOpenMP, kernel::SelectAll(), contactAccessor,
[&dem_collision, coefficientOfRestitution, dem_collisionTime, dem_kappa, dt](size_t c, data::ContactAccessor &ca, ParticleAccessorWithShape &pa){
[&dem_collision, coefficientOfRestitution, dem_collisionTime, dem_kappa, dt](size_t c, data::ContactAccessor &ca, data::ParticleAccessorWithBaseShape &pa){
auto idx1 = ca.getId1(c);
auto idx2 = ca.getId2(c);
auto meff = real_t(1) / (pa.getInvMass(idx1) + pa.getInvMass(idx2));
......@@ -1125,7 +1126,7 @@ int main(int argc, char **argv) {
// apply damping
particleStorage->forEachParticle(useOpenMP, kernel::SelectAll(), particleAccessor,
[velocityDampingFactor](size_t idx, ParticleAccessorWithShape &ac){
[velocityDampingFactor](size_t idx, data::ParticleAccessorWithBaseShape &ac){
ac.getLinearVelocityRef(idx) *= velocityDampingFactor;
ac.getAngularVelocityRef(idx) *= velocityDampingFactor;},
particleAccessor);
......
apps/showcases/ParticlePacking/Utility.h
View file @ b67091ed
......@@ -20,7 +20,11 @@
#pragma once
#include "core/mpi/MPITextFile.h"
#include "mesa_pd/data/ParticleStorage.h"
#include "mesa_pd/data/shape/HalfSpace.h"
#include "mesa_pd/data/shape/CylindricalBoundary.h"
#include <iterator>
#include <algorithm>
......@@ -29,22 +33,6 @@
namespace walberla {
namespace mesa_pd {
class ParticleAccessorWithShape : public data::ParticleAccessor
{
public:
ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& particleStorage)
: ParticleAccessor(particleStorage)
{}
walberla::real_t const & getInvMass(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInvMass();}
Mat3 const & getInvInertiaBF(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInvInertiaBF();}
Mat3 const & getInertiaBF(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInertiaBF();}
data::BaseShape* getShape(const size_t p_idx) const {return ps_->getBaseShape(p_idx).get();}
};
template< typename T>
std::vector<T> parseStringToVector(std::string inputStr)
{
......
python/mesa_pd/kernel/ExplicitEuler.py
View file @ b67091ed
......@@ -20,6 +20,7 @@ class ExplicitEuler:
self.context['interface'].append(create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"))
self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"))
self.context['interface'].append(create_access("inertiaBF", "walberla::mesa_pd::Mat3", access="g"))
self.context['interface'].append(create_access("torque", "walberla::mesa_pd::Vec3", access="gs"))
def generate(self, module):
......
python/mesa_pd/kernel/SemiImplicitEuler.py
View file @ b67091ed
......@@ -20,6 +20,7 @@ class SemiImplicitEuler:
self.context['interface'].append(create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"))
self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"))
self.context['interface'].append(create_access("inertiaBF", "walberla::mesa_pd::Mat3", access="g"))
self.context['interface'].append(create_access("torque", "walberla::mesa_pd::Vec3", access="gs"))
def generate(self, module):
......
python/mesa_pd/kernel/VelocityVerlet.py
View file @ b67091ed
......@@ -17,6 +17,7 @@ class VelocityVerlet:
self.context['interface'].append(create_access("rotation", "walberla::mesa_pd::Rot3", access="gs"))
self.context['interface'].append(create_access("angularVelocity", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("invInertiaBF", "walberla::mesa_pd::Mat3", access="g"))
self.context['interface'].append(create_access("inertiaBF", "walberla::mesa_pd::Mat3", access="g"))
self.context['interface'].append(create_access("torque", "walberla::mesa_pd::Vec3", access="gs"))
self.context['interface'].append(create_access("oldTorque", "walberla::mesa_pd::Vec3", access="gs"))
......
python/mesa_pd/templates/common/ParticleFunctions.templ.h
View file @ b67091ed
......@@ -67,6 +67,21 @@ inline Vec3 transformVectorFromBFtoWF(const size_t p_idx, Accessor& ac, const Ve
return ac.getRotation(p_idx).getMatrix() * vectorBF;
}
/**
* Transform (inverse) particle's moment of inertia from body frame coordinates (as stored by shape) to world frame.
*/
template <typename Accessor>
inline Mat3 getInvInertia(const size_t p_idx, Accessor& ac)
{
return math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(), ac.getInvInertiaBF(p_idx));
}
template <typename Accessor>
inline Mat3 getInertia(const size_t p_idx, Accessor& ac)
{
return math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(), ac.getInertiaBF(p_idx));
}
/**
* Force is applied at the center of mass.
*/
......
python/mesa_pd/templates/kernel/ExplicitEuler.templ.h
View file @ b67091ed
......@@ -28,6 +28,9 @@
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/IAccessor.h>
{%- if bIntegrateRotation %}
#include <mesa_pd/common/ParticleFunctions.h>
{%- endif %}
namespace walberla {
namespace mesa_pd {
......@@ -87,8 +90,12 @@ inline void ExplicitEuler::operator()(const size_t idx,
ac.getLinearVelocity(idx));
{%- if bIntegrateRotation %}
const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
// computation done in body frame: d(omega)/ dt = J^-1 ((J*omega) x omega + T), update in world frame
// see Wachs, 2019, doi:10.1007/s00707-019-02389-9, Eq. 27
const auto omegaBF = transformVectorFromWFtoBF(idx, ac, ac.getAngularVelocity(idx));
const auto torqueBF = transformVectorFromWFtoBF(idx, ac, ac.getTorque(idx));
const Vec3 wdotBF = ac.getInvInertiaBF(idx) * ( ( ac.getInertiaBF(idx) * omegaBF ) % omegaBF + torqueBF );
const Vec3 wdot = transformVectorFromBFtoWF(idx, ac, wdotBF);
// Calculating the rotation angle
const Vec3 phi( 0.5_r * wdot * dt_ * dt_ +
......
python/mesa_pd/templates/kernel/InitParticlesForHCSITS.templ.h
View file @ b67091ed
......@@ -87,7 +87,8 @@ inline void InitParticlesForHCSITS::operator()(size_t j, Accessor& ac, real_t dt
initializeVelocityCorrections( ac, j, ac.getDvRef(j), ac.getDwRef(j), dt ); // use applied external forces to calculate starting velocity
if(!isSet(particle_flags, FIXED)){ // Update velocities with global acceleration and angular velocity with euler eqn
ac.getLinearVelocityRef(j) = ac.getLinearVelocity(j) + getGlobalAcceleration() * dt;
ac.getAngularVelocityRef(j) = ac.getAngularVelocity(j) + dt * ( ac.getInvInertiaBF(j) * ( ( ac.getInertiaBF(j) * ac.getAngularVelocity(j) ) % ac.getAngularVelocity(j) ) );
const auto omegaBF = transformVectorFromWFtoBF(j, ac, ac.getAngularVelocity(j));
ac.getAngularVelocityRef(j) = ac.getAngularVelocity(j) + dt * transformVectorFromBFtoWF(j, ac, ( ac.getInvInertiaBF(j) * ( ( ac.getInertiaBF(j) * omegaBF ) % omegaBF ) ) );
}
}
}
......@@ -110,7 +111,7 @@ template <typename Accessor>
inline void InitParticlesForHCSITS::initializeVelocityCorrections(Accessor& ac, size_t body, Vec3& dv, Vec3& dw, real_t dt ) const
{
dv = ( ac.getInvMass(body) * dt ) * ac.getForce(body);
dw = dt * ( ac.getInvInertiaBF(body) * ac.getTorque(body) );
dw = dt * ( getInvInertia(body, ac) * ac.getTorque(body) );
ac.getForceRef(body) = Vec3();
ac.getTorqueRef(body) = Vec3();
......
python/mesa_pd/templates/kernel/SemiImplicitEuler.templ.h
View file @ b67091ed
......@@ -28,6 +28,9 @@
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/IAccessor.h>
{%- if bIntegrateRotation %}
#include <mesa_pd/common/ParticleFunctions.h>
{%- endif %}
namespace walberla {
namespace mesa_pd {
......@@ -81,8 +84,12 @@ inline void SemiImplicitEuler::operator()(const size_t idx,
ac.getPosition(idx));
{%- if bIntegrateRotation %}
const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
// computation done in body frame: d(omega)/ dt = J^-1 ((J*omega) x omega + T), update in world frame
// see Wachs, 2019, doi:10.1007/s00707-019-02389-9, Eq. 27
const auto omegaBF = transformVectorFromWFtoBF(idx, ac, ac.getAngularVelocity(idx));
const auto torqueBF = transformVectorFromWFtoBF(idx, ac, ac.getTorque(idx));
const Vec3 wdotBF = ac.getInvInertiaBF(idx) * ( ( ac.getInertiaBF(idx) * omegaBF ) % omegaBF + torqueBF );
const Vec3 wdot = transformVectorFromBFtoWF(idx, ac, wdotBF);
ac.setAngularVelocity(idx, wdot * dt_ +
ac.getAngularVelocity(idx));
......
python/mesa_pd/templates/kernel/VelocityVerlet.templ.h
View file @ b67091ed
......@@ -28,6 +28,9 @@
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/IAccessor.h>
{%- if bIntegrateRotation %}
#include <mesa_pd/common/ParticleFunctions.h>
{%- endif %}
namespace walberla {
namespace mesa_pd {
......@@ -95,8 +98,14 @@ inline void VelocityVerletPreForceUpdate::operator()(const size_t p_idx, Accesso
real_t(0.5) * ac.getInvMass(p_idx) * ac.getOldForce(p_idx) * dt_ * dt_);
{%- if bIntegrateRotation %}
const Vec3 wdot = math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(),
ac.getInvInertiaBF(p_idx)) * ac.getOldTorque(p_idx);
// computation done in body frame: d(omega)/ dt = J^-1 ((J*omega) x omega + T), update in world frame
// see Wachs, 2019, doi:10.1007/s00707-019-02389-9, Eq. 27
// note: this implementation (pre and post) is experimental as it is in principle unclear in which order
// angular velocities and rotations (affect again the transformations WF - BF) have to be carried out
const auto omegaBF = transformVectorFromWFtoBF(p_idx, ac, ac.getAngularVelocity(p_idx));
const auto torqueBF = transformVectorFromWFtoBF(p_idx, ac, ac.getOldTorque(p_idx));
const Vec3 wdotBF = ac.getInvInertiaBF(p_idx) * ( ( ac.getInertiaBF(p_idx) * omegaBF ) % omegaBF + torqueBF );
const Vec3 wdot = transformVectorFromBFtoWF(p_idx, ac, wdotBF);
// Calculating the rotation angle
const Vec3 phi( ac.getAngularVelocity(p_idx) * dt_ + real_t(0.5) * wdot * dt_ * dt_);
......@@ -120,12 +129,13 @@ inline void VelocityVerletPostForceUpdate::operator()(const size_t p_idx, Access
real_t(0.5) * ac.getInvMass(p_idx) * (ac.getOldForce(p_idx) + ac.getForce(p_idx)) * dt_);
{%- if bIntegrateRotation %}
const auto torque = ac.getOldTorque(p_idx) + ac.getTorque(p_idx);
const Vec3 wdot = math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(),
ac.getInvInertiaBF(p_idx)) * torque;
const auto omegaBF = transformVectorFromWFtoBF(p_idx, ac, ac.getAngularVelocity(p_idx));
const auto torqueBF = transformVectorFromWFtoBF(p_idx, ac, 0.5_r * (ac.getOldTorque(p_idx) + ac.getTorque(p_idx)));
const Vec3 wdotBF = ac.getInvInertiaBF(p_idx) * ( ( ac.getInertiaBF(p_idx) * omegaBF ) % omegaBF + torqueBF );
const Vec3 wdot = transformVectorFromBFtoWF(p_idx, ac, wdotBF);
ac.setAngularVelocity(p_idx, ac.getAngularVelocity(p_idx) +
real_t(0.5) * wdot * dt_ );
wdot * dt_ );
{%- endif %}
}
......
src/mesa_pd/common/ParticleFunctions.h
View file @ b67091ed
......@@ -67,6 +67,21 @@ inline Vec3 transformVectorFromBFtoWF(const size_t p_idx, Accessor& ac, const Ve
return ac.getRotation(p_idx).getMatrix() * vectorBF;
}
/**
* Transform (inverse) particle's moment of inertia from body frame coordinates (as stored by shape) to world frame.
*/
template <typename Accessor>
inline Mat3 getInvInertia(const size_t p_idx, Accessor& ac)
{
return math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(), ac.getInvInertiaBF(p_idx));
}
template <typename Accessor>
inline Mat3 getInertia(const size_t p_idx, Accessor& ac)
{
return math::transformMatrixRART(ac.getRotation(p_idx).getMatrix(), ac.getInertiaBF(p_idx));
}
/**
* Force is applied at the center of mass.
*/
......
apps/benchmarks/GranularGas/Accessor.h → src/mesa_pd/data/ParticleAccessorWithBaseShape.h
View file @ b67091ed
......@@ -13,33 +13,38 @@
// You should have received a copy of the GNU General Public License along
// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
//
//! \file Accessor.h
//! \author Sebastian Eibl <sebastian.eibl@fau.de>
//! \file
//! \author Christoph Rettinger <sebastian.eibl@fau.de>
//
//======================================================================================================================
#pragma once
#include <mesa_pd/data/ParticleAccessor.h>
#include <mesa_pd/data/ParticleStorage.h>
#include <mesa_pd/data/ShapeStorage.h>
namespace walberla {
namespace mesa_pd {
namespace data {
class ParticleAccessorWithShape : public data::ParticleAccessor
class ParticleAccessorWithBaseShape : public data::ParticleAccessor
{
public:
ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
ParticleAccessorWithBaseShape(std::shared_ptr<data::ParticleStorage>& ps)
: ParticleAccessor(ps)
, ss_(ss)
{}
const auto& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvMass();}
const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
const auto& getInvMass(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInvMass();}
const auto& getMass(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getMass();}
const auto& getInvInertiaBF(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInvInertiaBF();}
const auto& getInertiaBF(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getInertiaBF();}
auto getVolume(const size_t p_idx) const {return ps_->getBaseShapeRef(p_idx)->getVolume();}
data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)].get();}
private:
std::shared_ptr<data::ShapeStorage> ss_;
data::BaseShape* getShape(const size_t p_idx) const {return ps_->getBaseShape(p_idx).get();}
};
} // namespace mesa_pd
} // namespace walberla
} //namespace data
} //namespace mesa_pd
} //namespace walberla
src/mesa_pd/data/ParticleAccessorWithShape.h
View file @ b67091ed
......@@ -40,7 +40,7 @@ public:
const auto& getMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getMass();}
const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInvInertiaBF();}
const auto& getInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getInertiaBF();}
const auto getVolume(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getVolume();}
auto getVolume(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)]->getVolume();}
data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeIDRef(p_idx)].get();}
private:
......
src/mesa_pd/kernel/ExplicitEuler.h
View file @ b67091ed
......@@ -28,6 +28,7 @@
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/IAccessor.h>
#include <mesa_pd/common/ParticleFunctions.h>
namespace walberla {
namespace mesa_pd {
......@@ -64,6 +65,8 @@ namespace kernel {
*
* const walberla::mesa_pd::Mat3& getInvInertiaBF(const size_t p_idx) const;
*
* const walberla::mesa_pd::Mat3& getInertiaBF(const size_t p_idx) const;
*
* const walberla::mesa_pd::Vec3& getTorque(const size_t p_idx) const;
* void setTorque(const size_t p_idx, const walberla::mesa_pd::Vec3& v);
*
......@@ -97,8 +100,12 @@ inline void ExplicitEuler::operator()(const size_t idx,
ac.getPosition(idx));
ac.setLinearVelocity(idx, ac.getInvMass(idx) * ac.getForce(idx) * dt_ +
ac.getLinearVelocity(idx));
const Vec3 wdot = math::transformMatrixRART(ac.getRotation(idx).getMatrix(),
ac.getInvInertiaBF(idx)) * ac.getTorque(idx);
// computation done in body frame: d(omega)/ dt = J^-1 ((J*omega) x omega + T), update in world frame
// see Wachs, 2019, doi:10.1007/s00707-019-02389-9, Eq. 27
const auto omegaBF = transformVectorFromWFtoBF(idx, ac, ac.getAngularVelocity(idx));
const auto torqueBF = transformVectorFromWFtoBF(idx, ac, ac.getTorque(idx));
const Vec3 wdotBF = ac.getInvInertiaBF(idx) * ( ( ac.getInertiaBF(idx) * omegaBF ) % omegaBF + torqueBF );
const Vec3 wdot = transformVectorFromBFtoWF(idx, ac, wdotBF);
// Calculating the rotation angle
const Vec3 phi( 0.5_r * wdot * dt_ * dt_ +
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment