From 908018271a37f763d6addefb6ae5f450f3d4ad76 Mon Sep 17 00:00:00 2001
From: Christoph Rettinger <christoph.rettinger@fau.de>
Date: Thu, 10 Mar 2022 12:49:50 +0100
Subject: [PATCH] Showcase for dense particle packing
---
apps/showcases/CMakeLists.txt | 1 +
apps/showcases/ParticlePacking/CMakeLists.txt | 8 +
.../ParticlePacking/DiameterDistribution.h | 149 ++
apps/showcases/ParticlePacking/Evaluation.h | 599 ++++++++
.../ParticlePacking/ParticlePacking.cfg | 202 +++
.../ParticlePacking/ParticlePacking.cpp | 1330 +++++++++++++++++
.../ParticlePacking/ShapeGeneration.h | 566 +++++++
apps/showcases/ParticlePacking/Utility.h | 221 +++
.../example_meshes/sediment_scan_0.obj | 605 ++++++++
.../example_meshes/sediment_scan_1.obj | 545 +++++++
.../PegIntoSphereBed/PegIntoSphereBed.cpp | 4 +-
python/mesa_pd.py | 10 +
.../templates/mpi/ShapePackUnpack.templ.h | 4 +
src/mesa_pd/data/ParticleAccessor.h | 18 +
src/mesa_pd/data/ParticleStorage.h | 59 +
src/mesa_pd/data/shape/BaseShape.h | 3 +-
src/mesa_pd/kernel/InitParticlesForHCSITS.h | 2 +-
src/mesa_pd/kernel/IntegrateParticlesHCSITS.h | 1 -
src/mesa_pd/mpi/ShapePackUnpack.h | 4 +
.../notifications/NumContactNotification.h | 114 ++
.../notifications/ParticleCopyNotification.h | 4 +
.../ParticleGhostCopyNotification.h | 4 +
.../ConvexPolyhedron/MeshParticleVTKOutput.h | 28 +-
src/mesh_common/MeshOperations.h | 17 +
src/vtk/VTKTrait.h | 8 +
.../mesa_pd/vtk/ConvexPolyhedronVTKOutput.cpp | 4 +-
26 files changed, 4493 insertions(+), 17 deletions(-)
create mode 100644 apps/showcases/ParticlePacking/CMakeLists.txt
create mode 100644 apps/showcases/ParticlePacking/DiameterDistribution.h
create mode 100644 apps/showcases/ParticlePacking/Evaluation.h
create mode 100644 apps/showcases/ParticlePacking/ParticlePacking.cfg
create mode 100644 apps/showcases/ParticlePacking/ParticlePacking.cpp
create mode 100644 apps/showcases/ParticlePacking/ShapeGeneration.h
create mode 100644 apps/showcases/ParticlePacking/Utility.h
create mode 100644 apps/showcases/ParticlePacking/example_meshes/sediment_scan_0.obj
create mode 100644 apps/showcases/ParticlePacking/example_meshes/sediment_scan_1.obj
create mode 100644 src/mesa_pd/mpi/notifications/NumContactNotification.h
diff --git a/apps/showcases/CMakeLists.txt b/apps/showcases/CMakeLists.txt
index 6b34ee252..6330a83bd 100644
--- a/apps/showcases/CMakeLists.txt
+++ b/apps/showcases/CMakeLists.txt
@@ -4,6 +4,7 @@ add_subdirectory( FluidizedBed )
add_subdirectory( HeatConduction )
add_subdirectory( LightRisingParticleInFluidAMR )
add_subdirectory( Mixer )
+add_subdirectory( ParticlePacking )
add_subdirectory( PegIntoSphereBed )
if ( WALBERLA_BUILD_WITH_CODEGEN AND WALBERLA_BUILD_WITH_PYTHON )
add_subdirectory( PhaseFieldAllenCahn )
diff --git a/apps/showcases/ParticlePacking/CMakeLists.txt b/apps/showcases/ParticlePacking/CMakeLists.txt
new file mode 100644
index 000000000..abe7636fa
--- /dev/null
+++ b/apps/showcases/ParticlePacking/CMakeLists.txt
@@ -0,0 +1,8 @@
+waLBerla_link_files_to_builddir( *.cfg )
+
+if (WALBERLA_MESAPD_CONVEX_POLYHEDRON_AVAILABLE)
+ waLBerla_add_executable ( NAME ParticlePacking
+ FILES ParticlePacking.cpp
+ DEPENDS blockforest core mesa_pd postprocessing sqlite vtk )
+endif()
+
diff --git a/apps/showcases/ParticlePacking/DiameterDistribution.h b/apps/showcases/ParticlePacking/DiameterDistribution.h
new file mode 100644
index 000000000..27b4699e1
--- /dev/null
+++ b/apps/showcases/ParticlePacking/DiameterDistribution.h
@@ -0,0 +1,149 @@
+//======================================================================================================================
+//
+// 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 DiameterDistribution.h
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <random>
+
+#include "Utility.h"
+
+namespace walberla {
+namespace mesa_pd {
+
+
+class DiameterGenerator
+{
+public:
+ virtual real_t get() = 0;
+ virtual ~DiameterGenerator() = default;
+};
+
+class LogNormal: public DiameterGenerator
+{
+public:
+ LogNormal(real_t mu, real_t var, uint_t seed = 42)
+ : gen_(seed)
+ {
+ auto m = std::log(mu*mu/std::sqrt(var+mu*mu));
+ auto s = std::sqrt(log(var/(mu*mu) + 1_r));
+ dist_ = std::lognormal_distribution<> (m, s);
+ }
+
+ real_t get() override
+ {
+ return real_c(dist_(gen_));
+ }
+private:
+ std::mt19937 gen_;
+ std::lognormal_distribution<> dist_;
+};
+
+class Uniform: public DiameterGenerator
+{
+public:
+ Uniform(real_t diameter)
+ : diameter_(diameter)
+ { }
+
+ real_t get() override
+ {
+ return diameter_;
+ }
+private:
+ real_t diameter_;
+};
+
+class DiscreteSieving: public DiameterGenerator
+{
+public:
+ DiscreteSieving(std::vector<real_t> diameters, std::vector<real_t> massFractions, uint_t seed,
+ real_t normalParticleVolume, real_t totalParticleMass = 1_r, real_t particleDensity = 1_r)
+ : diameters_(diameters), gen_(seed)
+ {
+ WALBERLA_CHECK_EQUAL(diameters.size(), massFractions.size(), "Number of entries in diameter and mass-fraction array has to be the same!");
+ WALBERLA_CHECK_FLOAT_EQUAL(real_t(1), std::accumulate(massFractions.begin(), massFractions.end(), real_t(0)), "Sum of mass fractions has to be 1!");
+ auto particleNumbers = transferMassFractionsToParticleNumbers(massFractions, diameters, normalParticleVolume, totalParticleMass, particleDensity);
+ std::string outString = "Discrete Sieving: Expected particle numbers per diameter: | ";
+ bool issueWarning = false;
+ for(const auto & p : particleNumbers){
+ issueWarning |= (p > 0_r && p < 1_r);
+ outString += std::to_string(int(p)) + " | ";
+ }
+ auto totalParticles = std::accumulate(particleNumbers.begin(), particleNumbers.end(), real_t(0));
+ outString += " -> total = " + std::to_string(int(totalParticles));
+ WALBERLA_LOG_INFO_ON_ROOT(outString);
+ if(issueWarning) WALBERLA_LOG_WARNING_ON_ROOT("Attention: The simulated particle mass is not enough to have at least a single particle of all size classes!");
+ dist_ = std::discrete_distribution<uint_t>(particleNumbers.begin(), particleNumbers.end());
+ }
+
+ real_t get() override
+ {
+ return diameters_[dist_(gen_)];
+ }
+private:
+ std::vector<real_t> diameters_;
+ std::mt19937 gen_;
+ std::discrete_distribution<uint_t> dist_;
+};
+
+
+class ContinuousSieving: public DiameterGenerator
+{
+public:
+ ContinuousSieving(std::vector<real_t> sieveSizes, std::vector<real_t> massFractions, uint_t seed,
+ real_t normalParticleVolume, real_t totalParticleMass = 1_r, real_t particleDensity = 1_r)
+ : gen_(seed)
+ {
+ WALBERLA_CHECK_EQUAL(sieveSizes.size(), massFractions.size()+1, "Number of entries in sieves has to be one larger than the mass-fraction array!");
+ WALBERLA_CHECK_FLOAT_EQUAL(real_t(1), std::accumulate(massFractions.begin(), massFractions.end(), real_t(0)), "Sum of mass fractions has to be 1!");
+
+ auto meanDiameters = getMeanDiametersFromSieveSizes(sieveSizes);
+ auto particleNumbers = transferMassFractionsToParticleNumbers(massFractions, meanDiameters, normalParticleVolume, totalParticleMass, particleDensity);
+ std::string outString = "Continuous Sieving: Expected particle numbers per diameter: | ";
+ bool issueWarning = false;
+ for(const auto & p : particleNumbers){
+ issueWarning |= (p > 0_r && p < 1_r);
+ outString += std::to_string(int(p)) + " | ";
+ }
+ auto totalParticles = std::accumulate(particleNumbers.begin(), particleNumbers.end(), real_t(0));
+ outString += " -> total = " + std::to_string(int(totalParticles));
+ WALBERLA_LOG_INFO_ON_ROOT(outString);
+ if(issueWarning) WALBERLA_LOG_WARNING_ON_ROOT("Attention: The simulated particle mass is not enough to have at least a single particle of all size classes!");
+ std::vector<real_t> logSieveSizes(sieveSizes.size(),0_r); // = phi-scale -> uniformly distributed
+ for(uint_t i = 0; i < logSieveSizes.size(); ++i)
+ {
+ logSieveSizes[i] = - std::log2(sieveSizes[i]);
+ }
+
+ dist_ = std::piecewise_constant_distribution<real_t>(logSieveSizes.begin(), logSieveSizes.end(), particleNumbers.begin());
+ }
+
+ real_t get() override
+ {
+ return std::pow(2_r,-dist_(gen_)); //re-transfer from phi-scale to actual diameter
+ }
+private:
+ std::mt19937 gen_;
+ std::piecewise_constant_distribution<real_t> dist_;
+};
+
+
+} // namespace mesa_pd
+} // namespace walberla
diff --git a/apps/showcases/ParticlePacking/Evaluation.h b/apps/showcases/ParticlePacking/Evaluation.h
new file mode 100644
index 000000000..930fec6ac
--- /dev/null
+++ b/apps/showcases/ParticlePacking/Evaluation.h
@@ -0,0 +1,599 @@
+//======================================================================================================================
+//
+// 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 ParticlePacking.cpp
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "mesa_pd/data/ParticleStorage.h"
+
+#include "Utility.h"
+#include "ShapeGeneration.h"
+
+namespace walberla {
+namespace mesa_pd {
+
+
+struct ParticleInfo
+{
+ real_t averageVelocity = 0_r;
+ real_t maximumVelocity = 0_r;
+ uint_t numParticles = 0;
+ real_t maximumHeight = 0_r;
+ real_t particleVolume = 0_r;
+ real_t heightOfMass = 0_r;
+ real_t kinEnergy = 0_r; // = m * v^2/2
+ real_t meanCoordinationNumber = 0_r;
+
+ void allReduce()
+ {
+ walberla::mpi::allReduceInplace(numParticles, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(averageVelocity, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(maximumVelocity, walberla::mpi::MAX);
+ walberla::mpi::allReduceInplace(maximumHeight, walberla::mpi::MAX);
+ walberla::mpi::allReduceInplace(particleVolume, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(heightOfMass, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(kinEnergy, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(meanCoordinationNumber, walberla::mpi::SUM);
+
+ averageVelocity /= real_c(numParticles);
+ heightOfMass /= particleVolume;
+ meanCoordinationNumber /= real_c(numParticles);
+ }
+};
+
+std::ostream &operator<<(std::ostream &os, ParticleInfo const &m) {
+ return os << "Particle Info: uAvg = " << m.averageVelocity << ", uMax = " << m.maximumVelocity
+ << ", numParticles = " << m.numParticles << ", zMax = " << m.maximumHeight << ", Vp = "
+ << m.particleVolume << ", zMass = " << m.heightOfMass << ", kin. energy = " << m.kinEnergy << ", MCN = " << m.meanCoordinationNumber;
+}
+
+template< typename Accessor_T>
+ParticleInfo evaluateParticleInfo(const Accessor_T & ac)
+{
+
+ ParticleInfo info;
+ for(uint_t i = 0; i < ac.size(); ++i)
+ {
+ if (isSet(ac.getFlags(i), data::particle_flags::GHOST)) continue;
+ if (isSet(ac.getFlags(i), data::particle_flags::GLOBAL)) continue;
+
+ ++info.numParticles;
+ real_t velMagnitude = ac.getLinearVelocity(i).length();
+ real_t particleVolume = ac.getBaseShape(i)->getVolume();
+ real_t particleMass = ac.getBaseShape(i)->getMass();
+ real_t height = ac.getPosition(i)[2];
+ info.averageVelocity += velMagnitude;
+ info.maximumVelocity = std::max(info.maximumVelocity, velMagnitude);
+ info.maximumHeight = std::max(info.maximumHeight, height);
+ info.particleVolume += particleVolume;
+ info.heightOfMass += particleVolume*height;
+ info.kinEnergy += 0.5_r * particleMass * velMagnitude * velMagnitude;
+ info.meanCoordinationNumber += real_c(ac.getNumContacts(i));
+ }
+
+ info.allReduce();
+
+ return info;
+}
+
+struct ContactInfo
+{
+ real_t averagePenetrationDepth = 0_r;
+ real_t maximumPenetrationDepth = 0_r;
+ uint_t numContacts = 0;
+
+ void allReduce()
+ {
+ walberla::mpi::allReduceInplace(numContacts, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(averagePenetrationDepth, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(maximumPenetrationDepth, walberla::mpi::MAX);
+
+ if(numContacts > 0) averagePenetrationDepth /= real_t(numContacts);
+ }
+};
+
+std::ostream &operator<<(std::ostream &os, ContactInfo const &m) {
+ return os << "Contact Info: numContacts = " << m.numContacts << ", deltaAvg = " << m.averagePenetrationDepth << ", deltaMax = " << m.maximumPenetrationDepth;
+}
+
+ContactInfo evaluateContactInfo(const data::ContactAccessor & ca)
+{
+ ContactInfo info;
+ for(uint_t i = 0; i < ca.size(); ++i)
+ {
+ real_t penetrationDepth = -ca.getDistance(i);
+ info.maximumPenetrationDepth = std::max(info.maximumPenetrationDepth, penetrationDepth);
+ if(penetrationDepth > 0_r) {
+ info.averagePenetrationDepth += penetrationDepth;
+ ++info.numContacts;
+ }
+ }
+ info.allReduce();
+ return info;
+}
+
+class SizeEvaluator
+{
+public:
+ explicit SizeEvaluator(ScaleMode scaleMode) : scaleMode_(scaleMode) { }
+
+ struct SizeInfo
+ {
+ real_t size;
+ Vec3 shapeSemiAxes;
+ real_t volume;
+ };
+
+ SizeInfo get(data::BaseShape & bs)
+ {
+ SizeInfo sizeInfo;
+ sizeInfo.volume = bs.getVolume();
+
+ auto shapeType = bs.getShapeType();
+ if(shapeType == data::ConvexPolyhedron::SHAPE_TYPE) {
+ auto convexPolyhedron = static_cast<data::ConvexPolyhedron*>(&bs);
+ sizeInfo.shapeSemiAxes = semiAxesFromInertiaTensor(convexPolyhedron->getInertiaBF(), convexPolyhedron->getMass());
+ } else if(shapeType == data::Ellipsoid::SHAPE_TYPE) {
+ auto ellipsoid = static_cast<data::Ellipsoid*>(&bs);
+ sizeInfo.shapeSemiAxes = ellipsoid->getSemiAxes();
+ } else if(shapeType == data::Sphere::SHAPE_TYPE) {
+ auto sphere = static_cast<data::Sphere*>(&bs);
+ sizeInfo.shapeSemiAxes = Vec3(sphere->getRadius());
+ } else {
+ WALBERLA_ABORT("Shape handling not implemented!");
+ }
+
+ switch(scaleMode_)
+ {
+ case ScaleMode::sphereEquivalent:
+ {
+ sizeInfo.size = diameterFromSphereVolume(sizeInfo.volume);
+ break;
+ }
+ case ScaleMode::sieveLike:
+ {
+ sizeInfo.size = sizeFromSemiAxes(sizeInfo.shapeSemiAxes);
+ break;
+ }
+ default: WALBERLA_ABORT("Unknown shape scale mode");
+ }
+ return sizeInfo;
+ }
+
+private:
+ ScaleMode scaleMode_;
+};
+
+// shape evaluation functions
+real_t getFlatnessFromSemiAxes(Vec3 semiAxes)
+{
+ sortVector(semiAxes);
+ return semiAxes[0] / semiAxes[1];
+}
+
+real_t getElongationFromSemiAxes(Vec3 semiAxes)
+{
+ sortVector(semiAxes);
+ return semiAxes[1] / semiAxes[2];
+}
+
+real_t getEquancyFromSemiAxes(Vec3 semiAxes)
+{
+ sortVector(semiAxes);
+ return semiAxes[0] / semiAxes[2];
+}
+
+class ParticleHistogram
+{
+public:
+ ParticleHistogram(std::vector<real_t> sizeBins,
+ SizeEvaluator sizeEvaluator,
+ std::vector<std::vector<real_t>> shapeBins,
+ std::vector<std::tuple<std::string,std::function<real_t(Vec3)>>> shapeEvaluators) :
+ sizeBins_(sizeBins), sizeEvaluator_(sizeEvaluator),
+ shapeBins_(shapeBins), shapeEvaluators_(shapeEvaluators) {
+ bool areBinsAscending = (sizeBins[1] - sizeBins[0]) > real_t(0);
+ if(!areBinsAscending) std::reverse(sizeBins_.begin(), sizeBins_.end());
+
+ massFractionHistogram_ = std::vector<real_t>(sizeBins.size() + 1,real_t(0));
+ numberHistogram_ = std::vector<uint_t>(sizeBins.size() + 1,0);
+
+ WALBERLA_CHECK_EQUAL(shapeBins.size(), shapeEvaluators.size(), "Different number of shape evaluators and bins!");
+
+ shapeHistograms_ = std::vector<std::vector<uint_t>>(shapeBins.size());
+ for(uint_t i = 0; i < shapeBins.size(); ++i)
+ {
+ shapeHistograms_[i] = std::vector<uint_t>(shapeBins[i].size() + 1,0);
+ }
+ }
+
+ template<typename Accessor_T>
+ void operator()(size_t p_idx, Accessor_T& ac)
+ {
+ if (isSet(ac.getFlags(p_idx), data::particle_flags::GLOBAL) || isSet(ac.getFlags(p_idx), data::particle_flags::GHOST)) return;
+
+ auto sizeInfo = sizeEvaluator_.get(*ac.getBaseShape(p_idx));
+
+ auto size = sizeInfo.size;
+ auto volume = sizeInfo.volume;
+ auto semiAxes = sizeInfo.shapeSemiAxes;
+
+ auto sizeHistogramIdx = getHistogramIdx(sizeBins_, size);
+ massFractionHistogram_[sizeHistogramIdx] += volume;
+ numberHistogram_[sizeHistogramIdx]++;
+
+ for(uint_t i = 0; i < shapeBins_.size(); ++i)
+ {
+ auto shapeParam = std::get<1>(shapeEvaluators_[i])(semiAxes);
+ auto shapeHistogramIdx = getHistogramIdx(shapeBins_[i], shapeParam);
+ shapeHistograms_[i][shapeHistogramIdx]++;
+ }
+ }
+
+ void evaluate()
+ {
+ walberla::mpi::reduceInplace(numberHistogram_, walberla::mpi::SUM);
+ walberla::mpi::reduceInplace(massFractionHistogram_, walberla::mpi::SUM);
+ for(uint_t i = 0; i < shapeHistograms_.size(); ++i) walberla::mpi::reduceInplace(shapeHistograms_[i], walberla::mpi::SUM);
+ WALBERLA_ROOT_SECTION()
+ {
+ auto total = std::accumulate(massFractionHistogram_.begin(), massFractionHistogram_.end(), real_t(0));
+ for(auto& h: massFractionHistogram_) h /= total;
+ }
+ }
+
+ void clear()
+ {
+ for(auto& h: massFractionHistogram_) h = real_t(0);
+ for(auto& h: numberHistogram_) h = 0;
+ for(uint_t i = 0; i < shapeHistograms_.size(); ++i)
+ {
+ for(auto& h: shapeHistograms_[i]) h = 0;
+ }
+ }
+
+ std::vector<real_t> getSizeBins() const {return sizeBins_;}
+ std::vector<real_t> getMassFractionHistogram() const {return massFractionHistogram_;}
+ std::vector<uint_t> getNumberHistogram() const {return numberHistogram_;}
+
+ uint_t getNumberOfShapeEvaluators() const {return shapeEvaluators_.size();}
+ std::vector<real_t> getShapeBins(uint_t idx) const {return shapeBins_[idx];}
+ std::vector<uint_t> getShapeHistogram(uint_t idx) const {return shapeHistograms_[idx];}
+ std::tuple<std::string, std::function<real_t(Vec3)>> getShapeEvaluator(uint_t idx) const {return shapeEvaluators_[idx];}
+
+private:
+
+ uint_t getHistogramIdx(const std::vector<real_t> & bins, real_t value) const
+ {
+ auto numBins = bins.size();
+ if(value >= bins[numBins-1]) {
+ return numBins;
+ }
+ else {
+ for(uint_t b = 0; b < numBins; ++b){
+ if(value < bins[b])
+ {
+ return b;
+ }
+ }
+ }
+ WALBERLA_ABORT("No valid histogram bin found for value " << value);
+ return 0;
+ }
+
+
+ std::vector<real_t> sizeBins_;
+ SizeEvaluator sizeEvaluator_;
+ std::vector<std::vector<real_t>> shapeBins_;
+ std::vector<std::tuple<std::string,std::function<real_t(Vec3)>>> shapeEvaluators_; // vector of parameter label and function pairs
+
+ std::vector<real_t> massFractionHistogram_;
+ std::vector<uint_t> numberHistogram_;
+ std::vector<std::vector<uint_t>> shapeHistograms_;
+};
+
+std::ostream &operator<<(std::ostream &os, ParticleHistogram const &ph) {
+ auto bins = ph.getSizeBins();
+ os << "Size bins: | 0 - ";
+ for(auto b: bins) os << b*real_t(1e3) << " | " << b*real_t(1e3) << " - ";
+ os << "infty | \n";
+
+ auto hist = ph.getMassFractionHistogram();
+ os << "Mass Fraction Hist: | ";
+ for(auto h: hist) os << h << " | ";
+
+ os << "\n";
+
+ auto numHist = ph.getNumberHistogram();
+ os << "Number Hist: | ";
+ for(auto h: numHist) os << h << " | ";
+
+ for(uint_t i = 0; i < ph.getNumberOfShapeEvaluators(); ++i)
+ {
+ os << "\n\nShape parameter " << i << " ("<< std::get<0>(ph.getShapeEvaluator(i)) << ") bins: -infty - ";
+ for(auto b: ph.getShapeBins(i)) os << b << " | " << b << " - ";
+ os << "infty \n";
+ os << "Number Hist: | ";
+ for(auto h: ph.getShapeHistogram(i)) os << h << " | ";
+ }
+
+ return os;
+}
+
+class PorosityPerHorizontalLayerEvaluator
+{
+public:
+ PorosityPerHorizontalLayerEvaluator(real_t layerHeight, const AABB & simulationDomain, std::string domainSetup)
+ : layerHeight_(layerHeight)
+ {
+ if(domainSetup == "container") {
+ layerVolume_ = simulationDomain.xSize() * simulationDomain.xSize() * math::pi / real_t(4) * layerHeight;
+ }
+ else {
+ layerVolume_ = simulationDomain.xSize() * simulationDomain.ySize() * layerHeight;
+ }
+ porosityPerLayer_ = std::vector<real_t>(uint_c(simulationDomain.zSize() / layerHeight), real_t(0));
+ radiusPerLayer_ = std::vector<real_t>(uint_c(simulationDomain.zSize() / layerHeight), real_t(0));
+ }
+
+ template<typename Accessor_T>
+ void operator()(size_t p_idx, Accessor_T& ac)
+ {
+ if (isSet(ac.getFlags(p_idx), data::particle_flags::GLOBAL)) return;
+
+ real_t volume = ac.getBaseShapeRef(p_idx)->getVolume();
+ real_t radius = radiusFromSphereVolume(volume);
+ real_t zPos = ac.getPosition(p_idx)[2];
+ int heightIdxBegin = int(std::floor((zPos - radius)/layerHeight_));
+ int heightIdxEnd = int((zPos + radius)/layerHeight_) + 1;
+ for(int i = heightIdxBegin; i < heightIdxEnd; ++i)
+ {
+ real_t sphereSegmentVolume = calculateSphericalSegmentVolume(real_c(i) * layerHeight_ - zPos,
+ real_c(i+1) * layerHeight_ - zPos,
+ radius);
+ uint_t writeIdx = uint_c(std::min(std::max(i, 0), int(porosityPerLayer_.size() - 1) ) );
+ // i could be negative if overlap with bottom plane -> add this to idx 0 to not lose volume
+ // i could also be too large if evaluated during the simulation so cap its max value to avoid seg faults
+ porosityPerLayer_[writeIdx] += sphereSegmentVolume;
+ radiusPerLayer_[writeIdx] += sphereSegmentVolume * radius;
+
+ }
+ }
+
+ void evaluate()
+ {
+ walberla::mpi::reduceInplace(porosityPerLayer_, walberla::mpi::SUM);
+ walberla::mpi::reduceInplace(radiusPerLayer_, walberla::mpi::SUM);
+ WALBERLA_ROOT_SECTION()
+ {
+ for(uint_t i = 0; i < radiusPerLayer_.size(); ++i) radiusPerLayer_[i] /= porosityPerLayer_[i]; // = volume-averaged radius per layer
+ for(auto& p: porosityPerLayer_) p /= layerVolume_;
+ }
+ }
+
+ void clear()
+ {
+ for(uint_t i = 0; i < porosityPerLayer_.size(); ++i)
+ {
+ porosityPerLayer_[i] = 0_r;
+ radiusPerLayer_[i] = 0_r;
+ }
+ }
+
+ void printToFile(std::string fileName)
+ {
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ofstream file;
+ file.open( fileName.c_str() );
+ file << "#\t z\t porosity\t avgRadius\n";
+
+ for(uint_t i = 0; i < porosityPerLayer_.size(); ++i)
+ {
+ file << layerHeight_ * (real_t(0.5)+real_c(i)) << " " << porosityPerLayer_[i] << " " << radiusPerLayer_[i] << "\n";
+ }
+ file.close();
+ }
+ }
+
+ real_t estimateTotalPorosity()
+ {
+ /*
+ uint_t zMaxIdx = 0;
+ // find first empty index
+ while(porosityPerLayer_[zMaxIdx] > 0_r)
+ {
+ ++zMaxIdx;
+ }
+ zMaxIdx = uint_c(0.95_r * real_c(zMaxIdx)); // remove top 5% according to Schruff, 2018
+ return 1_r - std::accumulate(porosityPerLayer_.begin(), std::next(porosityPerLayer_.begin(), zMaxIdx), 0_r) / real_c(zMaxIdx+1);
+ */
+
+ const real_t cutOffPorosity = 0.5_r; // some value
+ const real_t cutOffPhi = 1_r-cutOffPorosity; // solid volume fraction
+
+ uint_t endEvalIdx = 0;
+ uint_t numLayers = porosityPerLayer_.size();
+ for(uint_t i = numLayers-1; i > 0; --i)
+ {
+ if(porosityPerLayer_[i] <= cutOffPhi && porosityPerLayer_[i-1] >= cutOffPhi)
+ {
+ endEvalIdx = i;
+ break;
+ }
+ }
+ if(endEvalIdx > 0) return 1_r - std::accumulate(porosityPerLayer_.begin(), std::next(porosityPerLayer_.begin(), static_cast<long>(endEvalIdx)), 0_r) / real_c(endEvalIdx);
+ else return 1_r;
+
+ }
+
+private:
+
+ real_t calculateSphericalSegmentVolume(real_t lowerLayerHeight, real_t upperLayerHeight, real_t radius) const
+ {
+ real_t r1 = std::max(std::min(lowerLayerHeight,radius),-radius);
+ real_t r2 = std::max(std::min(upperLayerHeight,radius),-radius);
+ real_t a1 = std::sqrt(std::max(radius*radius - r1*r1, real_t(0)));
+ real_t a2 = std::sqrt(std::max(radius*radius - r2*r2, real_t(0)));
+ real_t h = r2-r1;
+ // wikipedia Kugelschicht
+ return math::pi * h / real_t(6) * (real_t(3)*a1*a1 + real_t(3)*a2*a2 + h*h);
+ }
+
+ real_t layerHeight_;
+ real_t layerVolume_;
+ std::vector<real_t> porosityPerLayer_;
+ std::vector<real_t> radiusPerLayer_;
+};
+
+
+class ContactInfoPerHorizontalLayerEvaluator
+{
+public:
+ ContactInfoPerHorizontalLayerEvaluator(real_t layerHeight, const AABB & simulationDomain)
+ : layerHeight_(layerHeight)
+ {
+ contactsPerLayer_ = std::vector<uint_t>(uint_c(simulationDomain.zSize() / layerHeight), 0);
+ avgPenetrationPerLayer_ = std::vector<real_t>(uint_c(simulationDomain.zSize() / layerHeight), 0_r);
+ maxPenetrationPerLayer_ = std::vector<real_t>(uint_c(simulationDomain.zSize() / layerHeight), 0_r);
+ }
+
+ template<typename ContactAccessor_T>
+ void operator()(size_t c_idx, ContactAccessor_T& cac)
+ {
+ real_t zPos = cac.getPosition(c_idx)[2];
+ real_t penetrationDepth = -cac.getDistance(c_idx);
+ uint_t heightIdx= uint_c(std::max(0_r,zPos / layerHeight_));
+ contactsPerLayer_[heightIdx]++;
+ avgPenetrationPerLayer_[heightIdx] += penetrationDepth;
+ maxPenetrationPerLayer_[heightIdx] = std::max(maxPenetrationPerLayer_[heightIdx], penetrationDepth);
+ }
+
+ void evaluate()
+ {
+ walberla::mpi::reduceInplace(contactsPerLayer_, walberla::mpi::SUM);
+ walberla::mpi::reduceInplace(avgPenetrationPerLayer_, walberla::mpi::SUM);
+ walberla::mpi::reduceInplace(maxPenetrationPerLayer_, walberla::mpi::MAX);
+ WALBERLA_ROOT_SECTION()
+ {
+ for(uint_t i = 0; i < contactsPerLayer_.size(); ++i)
+ {
+ if(contactsPerLayer_[i] > 0)
+ {
+ avgPenetrationPerLayer_[i] /= real_c(contactsPerLayer_[i]);
+ }
+ }
+ }
+ }
+
+ void clear()
+ {
+ for(uint_t i = 0; i < contactsPerLayer_.size(); ++i)
+ {
+ contactsPerLayer_[i] = 0;
+ avgPenetrationPerLayer_[i] = 0_r;
+ maxPenetrationPerLayer_[i] = 0_r;
+ }
+ }
+
+ void printToFile(std::string fileName)
+ {
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ofstream file;
+ file.open( fileName.c_str() );
+ file << "#\t z\t numContacts\t avgPenetration\t maxPenetration\n";
+
+ for(uint_t i = 0; i < contactsPerLayer_.size(); ++i)
+ {
+ file << layerHeight_ * (real_t(0.5)+real_c(i)) << " " << contactsPerLayer_[i]
+ << " " << avgPenetrationPerLayer_[i] << " " << maxPenetrationPerLayer_[i] << "\n";
+ }
+ file.close();
+ }
+ }
+
+private:
+
+ real_t layerHeight_;
+ std::vector<uint_t> contactsPerLayer_;
+ std::vector<real_t> avgPenetrationPerLayer_;
+ std::vector<real_t> maxPenetrationPerLayer_;
+};
+
+
+class LoggingWriter
+{
+public:
+ explicit LoggingWriter(std::string fileName) : fileName_(fileName){
+ WALBERLA_ROOT_SECTION() {
+ std::ofstream file;
+ file.open(fileName_.c_str());
+ file << "# t numParticles maxVel avgVel maxHeight massHeight kinEnergy numContacts maxPenetration avgPenetration porosity MCN\n";
+ file.close();
+ }
+ }
+
+ void operator()(real_t t, const ParticleInfo & particleInfo, const ContactInfo & contactInfo, real_t porosity)
+ {
+ WALBERLA_ROOT_SECTION() {
+ std::ofstream file;
+ file.open(fileName_.c_str(), std::ios_base::app);
+ file << t << " " << particleInfo.numParticles << " " << particleInfo.maximumVelocity << " "
+ << particleInfo.averageVelocity << " " << particleInfo.maximumHeight << " " << particleInfo.heightOfMass << " "
+ << particleInfo.kinEnergy << " " << contactInfo.numContacts << " " << contactInfo.maximumPenetrationDepth << " "
+ << contactInfo.averagePenetrationDepth << " " << porosity << " " << particleInfo.meanCoordinationNumber << "\n";
+ file.close();
+ }
+ }
+
+private:
+ std::string fileName_;
+};
+
+
+std::string assembleParticleInformation(data::ParticleStorage& ps, SizeEvaluator sizeEvaluator, int precision)
+{
+ std::ostringstream ossData;
+
+ for (auto pIt : ps)
+ {
+
+ using namespace data::particle_flags;
+ if( isSet(pIt->getFlags(), GLOBAL) || isSet(pIt->getFlags(), GHOST)) continue;
+
+ auto position = pIt->getPosition();
+ auto sizeInfo = sizeEvaluator.get(*pIt->getBaseShape());
+
+ ossData << std::setprecision( precision );
+ ossData << position[0] << " " << position[1] << " " << position[2] << " "
+ << sizeInfo.size << " " << sizeInfo.volume << " "
+ << sizeInfo.shapeSemiAxes[0] << " " << sizeInfo.shapeSemiAxes[1] << " " << sizeInfo.shapeSemiAxes[2] << " "
+ << pIt->getNumContacts()
+ << "\n";
+ }
+
+ return ossData.str();
+}
+
+
+
+} // namespace msa_pd
+} // namespace walberla
+
diff --git a/apps/showcases/ParticlePacking/ParticlePacking.cfg b/apps/showcases/ParticlePacking/ParticlePacking.cfg
new file mode 100644
index 000000000..e1d4d9e59
--- /dev/null
+++ b/apps/showcases/ParticlePacking/ParticlePacking.cfg
@@ -0,0 +1,202 @@
+ParticlePacking
+{
+ domainSetup periodic; // container, periodic
+ domainWidth 0.104; // m
+ domainHeight 1; // m
+
+ particleDensity 2650; // kg/m^3, 2500 (spheres), 2650 (real sediments);
+ ambientDensity 1000; // kg/m^3
+ gravitationalAcceleration 9.81; // m/s^2
+
+ particleDistribution SievingCurve; // size distribution, see 'Distribution' block for options
+ particleShape Mesh; // see 'Shape' block
+
+ limitVelocity -1; // m/s, negative switches limiting off
+ initialVelocity 1; // m/s
+ initialGenerationHeightRatioStart 0.1; // -
+ initialGenerationHeightRatioEnd 1; // -
+ generationSpacing 0.010; // m
+ scaleGenerationSpacingWithForm true;
+ generationHeightRatioStart 0.6; // -
+ generationHeightRatioEnd 1; // -
+
+ totalParticleMass 4; // kg
+
+ numBlocksPerDirection <3,3,1>;
+
+ visSpacing 0.01; // s
+ infoSpacing 0.01; // s
+ loggingSpacing 0.01; // s
+
+ terminalVelocity 1e-3; // m/s
+ terminalRelativeHeightChange 1e-5; // -
+ terminationCheckingSpacing 0.01; // s
+ minimalTerminalRunTime .0; // s, minimal runtime after generation & shaking
+
+ solver DEM; // particle simulation approach, see 'Solver' block for options
+
+ shaking true; // see 'Shaking' block
+
+ velocityDampingCoefficient 1e-3; // continuous reduction of velocity in last simulation phase
+
+ useHashGrids false;
+ particleSortingSpacing 1000; // time steps, non-positive values switch sorting off, performance optimization
+}
+
+Shaking
+{
+ amplitude 3e-4; // m
+ period 0.025; // s
+ duration 6.0; // s, duration of shaking AFTER creation of all particles
+ activeFromBeginning true;
+}
+
+Solver
+{
+ coefficientOfRestitution 0.1; // -
+ frictionCoefficient 0.5; // -
+ dt 5e-6; // s, time step size
+ DEM
+ {
+ collisionTime 20e-5; // s
+ poissonsRatio 0.22; // -
+ }
+
+ HCSITS
+ {
+ errorReductionParameter 0.8;
+ relaxationParameter 0.75;
+ numberOfIterations 10;
+ relaxationModel InelasticGeneralizedMaximumDissipationContact;
+ }
+}
+
+Distribution
+{
+ randomSeed 41; // if negative, seed is randomly chosen
+
+ Uniform{
+ diameter 11e-3; // m
+ }
+
+ LogNormal{
+ mu 1e-3; // m
+ variance 1e-7;
+ }
+
+ DiameterMassFractions
+ {
+ // from Schruff, 2016
+ //diameters 1.5e-3 2e-3 3e-3 4e-3 6e-3 8e-3 11e-3 16e-3 22e-3;
+ //massFractions .00 .00 1.00 .00 .00 .00 .00 .00 .00; // I, U
+ //massFractions .00 .00 .00 1.00 .00 .00 .00 .00 .00; // II, U
+ //massFractions .00 .00 .00 .00 1.00 .00 .00 .00 .00; // III, U
+ //massFractions .00 .00 .50 .00 .50 .00 .00 .00 .00; // IV, B
+ //massFractions .00 .00 .00 .50 .00 .50 .00 .00 .00; // V, B
+ //massFractions .00 .00 .00 .00 .00 .50 .00 .50 .00; // VI, B
+ //massFractions .00 .00 .25 .50 .25 .00 .00 .00 .00; // VII, T
+ //massFractions .00 .00 .00 .25 .50 .25 .00 .00 .00; // VIII, T
+ //massFractions .00 .00 .00 .00 .25 .50 .25 .00 .00; // IX, T
+ //massFractions .025 .050 .100 .200 .250 .200 .100 .050 .025; // X, LN
+
+ // Liang, 2015, discrete spheres (3,4,6,8,11,16,22 mm)
+ diameters 3e-3 4e-3 6e-3 8e-3 11e-3 16e-3 22e-3;
+ //massFractions 0 0 1 0 0 0 0; // A
+ //massFractions .0 .0 .21 .58 .21 0 0; // B
+ //massFractions .0 .06 .24 .4 .24 .06 .0; // C
+ //massFractions .04 .11 .22 .26 .22 .11 .04; // D
+ //massFractions .08 .13 .08 .06 .18 .29 .18; // E
+ massFractions .13 .21 .13 .06 .13 .21 .13; // F
+ //massFractions .18 .29 .18 .06 .08 .13 .08; // G
+ }
+
+ SievingCurve
+ {
+ // Frings, 2011, Verification, slightly truncated, table 2
+ //sieveSizes 0.063e-3 0.09e-3 0.125e-3 0.18e-3 0.25e-3 0.355e-3 0.5e-3 0.71e-3 1e-3 1.4e-3 2e-3 2.8e-3 4e-3 5.6e-3 8e-3 11.2e-3 16e-3 22.4e-3 31.5e-3 45e-3 63e-3;
+ //massFractions .00 .00 .01 .01 .05 .15 .12 .06 .03 .03 .03 .05 .06 .07 .07 .07 .06 .06 .04 .03; // case 45
+
+ // Lian, 2015, DigiPack, values from porosity report (Versuch3), in cylinder with diam 104mm
+ sieveSizes 2.8e-3 4e-3 5.6e-3 8e-3 11.2e-3 16e-3 22.4e-3 31.5e-3;
+ //massFractions 0 0 1 0 0 0 0; // A
+ //massFractions .0 .0 .21 .58 .21 0 0; // B
+ //massFractions .0 .06 .24 .4 .24 .06 .0; // C
+ //massFractions .04 .11 .22 .26 .22 .11 .04; // D
+ //massFractions .08 .13 .08 .06 .18 .29 .18; // E
+ massFractions .13 .21 .13 .06 .13 .21 .13; // F
+ //massFractions .18 .29 .18 .06 .08 .13 .08; // G
+
+ //massFractions 0 0 0 0 0 0 1; // D
+
+ useDiscreteForm false; // only use average sieving diameters
+ }
+}
+
+Shape
+{
+ scaleMode sieveLike; // sphereEquivalent, sieveLike
+
+ Sphere
+ {
+ }
+
+ Ellipsoid
+ {
+ semiAxes <50,5,1>; // will be scaled to obtain desired size
+ }
+
+ EquivalentEllipsoid
+ {
+ // specify either a mesh file or a folder containing mesh files ( have to be .obj files)
+ //path example_meshes/sediment_scan_0.obj;
+ path example_meshes;
+ }
+
+ EllipsoidFormDistribution
+ {
+ elongationMean 0.4;
+ elongationStdDev 0.2;
+ flatnessMean 0.4;
+ flatnessStdDev 0.2;
+ }
+
+ Mesh
+ {
+ // specify either a mesh file or a folder containing mesh files ( have to be .obj files)
+ //path example_meshes/sediment_scan_0.obj;
+ path example_meshes;
+ }
+
+ MeshFormDistribution
+ {
+ path example_meshes; // meshes used as basis, are then scaled to match size and form
+
+ elongationMean 0.4; //0.682;
+ elongationStdDev 0.2; //0.139;
+ flatnessMean 0.4; //0.65;
+ flatnessStdDev 0.2; //0.175;
+ }
+
+ UnscaledMeshesPerFraction
+ {
+ // expects subfolders 0, 1,... in given folder that contain a set of meshes
+ // those meshes are taken as is and not scaled during creation
+ // the mass fraction from Distribution/SievingCurve is used to determine the generation probabilities per fraction
+ folder meshes_collection_unscaled;
+ }
+}
+
+Evaluation
+{
+ vtkFolder vtk_out_container;
+ vtkFinalFolder vtk_files;
+
+ //Rhein, Frings, 2011, Verification, Table 2
+ histogramBins 200e-3 125e-3 90e-3 63e-3 45e-3 31.5e-3 22.4e-3 16e-3 11.2e-3 8e-3 5.6e-3 4e-3 2.8e-3 2e-3 1.4e-3 1e-3 0.71e-3 0.5e-3 0.355e-3 0.25e-3 0.18e-3 0.125e-3 0.09e-3 0.063e-3;
+
+ //for horizontal layer evaluation (only spheres)
+ porosityProfileFolder porosity_profiles;
+ layerHeight 1e-3; // m
+
+ sqlDBFileName db_ParticlePacking.sqlite;
+}
diff --git a/apps/showcases/ParticlePacking/ParticlePacking.cpp b/apps/showcases/ParticlePacking/ParticlePacking.cpp
new file mode 100644
index 000000000..41f97f8b9
--- /dev/null
+++ b/apps/showcases/ParticlePacking/ParticlePacking.cpp
@@ -0,0 +1,1330 @@
+//======================================================================================================================
+//
+// 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 ParticlePacking.cpp
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+
+#include "blockforest/Initialization.h"
+#include "blockforest/BlockForest.h"
+#include "core/Environment.h"
+#include "core/grid_generator/SCIterator.h"
+#include "core/grid_generator/HCPIterator.h"
+#include "core/math/all.h"
+#include "core/timing/TimingTree.h"
+#include "vtk/VTKOutput.h"
+
+#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
+#include "mesa_pd/collision_detection/GeneralContactDetection.h"
+
+#include "mesa_pd/common/ParticleFunctions.h"
+
+#include "mesa_pd/data/ContactStorage.h"
+#include "mesa_pd/data/ContactAccessor.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/HashGrids.h"
+
+#include "mesa_pd/domain/BlockForestDomain.h"
+
+#include "mesa_pd/kernel/AssocToBlock.h"
+#include "mesa_pd/kernel/InsertParticleIntoLinkedCells.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/kernel/DetectAndStoreContacts.h"
+#include "mesa_pd/kernel/InitContactsForHCSITS.h"
+#include "mesa_pd/kernel/InitParticlesForHCSITS.h"
+#include "mesa_pd/kernel/IntegrateParticlesHCSITS.h"
+#include "mesa_pd/kernel/HCSITSRelaxationStep.h"
+#include "mesa_pd/kernel/SemiImplicitEuler.h"
+#include "mesa_pd/kernel/LinearSpringDashpot.h"
+
+#include "mesa_pd/mpi/ContactFilter.h"
+#include "mesa_pd/mpi/ReduceProperty.h"
+#include "mesa_pd/mpi/ReduceContactHistory.h"
+#include "mesa_pd/mpi/BroadcastProperty.h"
+#include "mesa_pd/mpi/SyncNextNeighborsBlockForest.h"
+#include "mesa_pd/mpi/SyncGhostOwners.h"
+#include "mesa_pd/mpi/notifications/VelocityUpdateNotification.h"
+#include "mesa_pd/mpi/notifications/VelocityCorrectionNotification.h"
+#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
+#include "mesa_pd/mpi/notifications/NumContactNotification.h"
+
+#include "mesa_pd/sorting/LinearizedCompareFunctor.h"
+
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+#include "mesa_pd/vtk/TensorGlyph.h"
+#include "mesa_pd/vtk/ConvexPolyhedron/MeshParticleVTKOutput.h"
+#include "mesa_pd/vtk/ConvexPolyhedron/data_sources/SurfaceVelocityVertexDataSource.h"
+
+#include "sqlite/SQLite.h"
+
+#include <iostream>
+#include <random>
+#include <chrono>
+
+#include "Utility.h"
+#include "Evaluation.h"
+#include "DiameterDistribution.h"
+#include "ShapeGeneration.h"
+
+namespace walberla {
+namespace mesa_pd {
+
+
+kernel::HCSITSRelaxationStep::RelaxationModel relaxationModelFromString(const std::string & model)
+{
+ if(model == "InelasticFrictionlessContact")
+ return kernel::HCSITSRelaxationStep::RelaxationModel::InelasticFrictionlessContact;
+ if(model == "ApproximateInelasticCoulombContactByDecoupling")
+ return kernel::HCSITSRelaxationStep::RelaxationModel::ApproximateInelasticCoulombContactByDecoupling;
+ if(model == "ApproximateInelasticCoulombContactByOrthogonalProjections")
+ return kernel::HCSITSRelaxationStep::RelaxationModel::ApproximateInelasticCoulombContactByOrthogonalProjections;
+ if(model == "InelasticCoulombContactByDecoupling")
+ return kernel::HCSITSRelaxationStep::RelaxationModel::InelasticCoulombContactByDecoupling;
+ if(model == "InelasticCoulombContactByOrthogonalProjections")
+ return kernel::HCSITSRelaxationStep::RelaxationModel::InelasticCoulombContactByOrthogonalProjections;
+ if(model == "InelasticGeneralizedMaximumDissipationContact")
+ return kernel::HCSITSRelaxationStep::RelaxationModel::InelasticGeneralizedMaximumDissipationContact;
+ if(model == "InelasticProjectedGaussSeidel")
+ return kernel::HCSITSRelaxationStep::RelaxationModel::InelasticProjectedGaussSeidel;
+ WALBERLA_ABORT("Unknown relaxation model " << model);
+}
+
+class ParticleCreator
+{
+public:
+ ParticleCreator(const std::shared_ptr<data::ParticleStorage> & particleStorage, const std::shared_ptr<domain::IDomain> & particleDomain,
+ const AABB & simulationDomain, const std::string & domainSetup, real_t particleDensity, bool scaleGenerationSpacingWithForm ) :
+ particleStorage_(particleStorage), particleDomain_(particleDomain), simulationDomain_(simulationDomain),
+ domainSetup_(domainSetup), particleDensity_(particleDensity), scaleGenerationSpacingWithForm_(scaleGenerationSpacingWithForm),
+ gen_(static_cast<unsigned long>(walberla::mpi::MPIManager::instance()->rank()))
+ { }
+
+ void createParticles( real_t zMin, real_t zMax, real_t spacing,
+ const shared_ptr<DiameterGenerator> & diameterGenerator, const shared_ptr<ShapeGenerator> & shapeGenerator,
+ real_t initialVelocity, real_t maximumAllowedInteractionRadius )
+ {
+ // this scaling is done to flexibly change the generation scaling in x,y, and z direction, based on the average form
+ auto spacingScaling = (scaleGenerationSpacingWithForm_) ? shapeGenerator->getNormalFormParameters() / shapeGenerator->getNormalFormParameters()[1] // divide by I for normalization
+ : Vec3(1_r); // no scaling (= equal spacing) in all directions
+ sortVector(spacingScaling); // S, I, L
+ Vec3 invScaling(1_r/spacingScaling[0], 1_r/spacingScaling[1], 1_r/spacingScaling[2]);
+
+ AABB creationDomain(simulationDomain_.xMin()*invScaling[0], simulationDomain_.yMin()*invScaling[1], zMin*invScaling[2],
+ simulationDomain_.xMax()*invScaling[0], simulationDomain_.yMax()*invScaling[1], zMax*invScaling[2]);
+ Vec3 pointOfReference(0,0,(zMax+zMin)*0.5_r*invScaling[2]);
+
+ WALBERLA_LOG_INFO_ON_ROOT("Creating particles between z = " << zMin << " and " << zMax);
+ for (auto ptUnscaled : grid_generator::HCPGrid(creationDomain.getExtended(Vec3(-0.5_r * spacing, -0.5_r * spacing, 0_r)), pointOfReference, spacing))
+ {
+ Vec3 pt(ptUnscaled[0]*spacingScaling[0], ptUnscaled[1]*spacingScaling[1], ptUnscaled[2]*spacingScaling[2]); // scale back
+ auto diameter = diameterGenerator->get();
+ if(!particleDomain_->isContainedInLocalSubdomain(pt,real_t(0))) continue;
+ if(domainSetup_ == "container")
+ {
+ auto domainCenter = simulationDomain_.center();
+ auto distanceFromDomainCenter = pt - domainCenter;
+ distanceFromDomainCenter[2] = real_t(0);
+ auto distance = distanceFromDomainCenter.length();
+ real_t containerRadius = real_t(0.5)*simulationDomain_.xSize();
+ if(distance > containerRadius - real_t(0.5) * spacing) continue;
+ }
+
+ // create particle
+ auto p = particleStorage_->create();
+ p->getPositionRef() = pt;
+
+ shapeGenerator->setShape(diameter, maximumAllowedInteractionRadius, p->getBaseShapeRef(), p->getInteractionRadiusRef());
+
+ p->getBaseShapeRef()->updateMassAndInertia(particleDensity_);
+
+ p->setLinearVelocity( Vec3(0.1_r*math::realRandom(-initialVelocity, initialVelocity, gen_),
+ 0.1_r*math::realRandom(-initialVelocity, initialVelocity, gen_),
+ -initialVelocity) );
+
+ p->setAngularVelocity( 0.1_r*Vec3(math::realRandom(-initialVelocity, initialVelocity, gen_),
+ math::realRandom(-initialVelocity, initialVelocity, gen_),
+ math::realRandom(-initialVelocity, initialVelocity, gen_)) / diameter );
+
+ p->getOwnerRef() = walberla::mpi::MPIManager::instance()->rank();
+ p->getTypeRef() = 0;
+ }
+ }
+
+private:
+ std::shared_ptr<data::ParticleStorage> particleStorage_;
+ std::shared_ptr<domain::IDomain> particleDomain_;
+ AABB simulationDomain_;
+ std::string domainSetup_;
+ real_t particleDensity_;
+ bool scaleGenerationSpacingWithForm_;
+ std::mt19937 gen_;
+};
+
+void addConfigToDatabase(Config & config,
+ std::map< std::string, walberla::int64_t > & integerProperties,
+ std::map< std::string, double > & realProperties,
+ std::map< std::string, std::string > & stringProperties)
+{
+
+ const Config::BlockHandle mainConf = config.getBlock("ParticlePacking");
+
+ Vector3<uint_t> numBlocksPerDirection = mainConf.getParameter< Vector3<uint_t> >("numBlocksPerDirection");
+ integerProperties["numBlocksX"] = int64_c(numBlocksPerDirection[0]);
+ integerProperties["numBlocksY"] = int64_c(numBlocksPerDirection[1]);
+ integerProperties["numBlocksZ"] = int64_c(numBlocksPerDirection[2]);
+ integerProperties["useHashGrids"] = (mainConf.getParameter<bool>("useHashGrids")) ? 1 : 0;
+ integerProperties["scaleGenerationSpacingWithForm"] = (mainConf.getParameter<bool>("scaleGenerationSpacingWithForm")) ? 1 : 0;
+ stringProperties["domainSetup"] = mainConf.getParameter<std::string>("domainSetup");
+ stringProperties["particleDistribution"] = mainConf.getParameter<std::string>("particleDistribution");
+ stringProperties["particleShape"] = mainConf.getParameter<std::string>("particleShape");
+ stringProperties["solver"] = mainConf.getParameter<std::string>("solver");
+ realProperties["domainWidth"] = mainConf.getParameter<double>("domainWidth");
+ realProperties["domainHeight"] = mainConf.getParameter<double>("domainHeight");
+ realProperties["particleDensity"] = mainConf.getParameter<double>("particleDensity");
+ realProperties["ambientDensity"] = mainConf.getParameter<double>("ambientDensity");
+ realProperties["gravitationalAcceleration"] = mainConf.getParameter<double>("gravitationalAcceleration");
+ realProperties["limitVelocity"] = mainConf.getParameter<double>("limitVelocity");
+ realProperties["initialVelocity"] = mainConf.getParameter<double>("initialVelocity");
+ realProperties["initialGenerationHeightRatioStart"] = mainConf.getParameter<double>("initialGenerationHeightRatioStart");
+ realProperties["initialGenerationHeightRatioEnd"] = mainConf.getParameter<double>("initialGenerationHeightRatioEnd");
+ realProperties["generationSpacing"] = mainConf.getParameter<double>("generationSpacing");
+ realProperties["generationHeightRatioStart"] = mainConf.getParameter<double>("generationHeightRatioStart");
+ realProperties["generationHeightRatioEnd"] = mainConf.getParameter<double>("generationHeightRatioEnd");
+ realProperties["totalParticleMass"] = mainConf.getParameter<double>("totalParticleMass");
+ realProperties["terminalVelocity"] = mainConf.getParameter<double>("terminalVelocity");
+ realProperties["terminalRelativeHeightChange"] = mainConf.getParameter<double>("terminalRelativeHeightChange");
+ realProperties["minimalTerminalRunTime"] = mainConf.getParameter<double>("minimalTerminalRunTime");
+ realProperties["terminationCheckingSpacing"] = mainConf.getParameter<double>("terminationCheckingSpacing");
+ realProperties["velocityDampingCoefficient"] = mainConf.getParameter<double>("velocityDampingCoefficient");
+
+ const Config::BlockHandle solverConf = config.getBlock("Solver");
+ realProperties["dt"] = solverConf.getParameter<double>("dt");
+ realProperties["frictionCoefficient"] = solverConf.getParameter<double>("frictionCoefficient");
+ realProperties["coefficientOfRestitution"] = solverConf.getParameter<double>("coefficientOfRestitution");
+ const Config::BlockHandle solverHCSITSConf = solverConf.getBlock("HCSITS");
+ integerProperties["hcsits_numberOfIterations"] = solverHCSITSConf.getParameter<int64_t>("numberOfIterations");
+ stringProperties["hcsits_relaxationModel"] = solverHCSITSConf.getParameter<std::string>("relaxationModel");
+ realProperties["hcsits_errorReductionParameter"] = solverHCSITSConf.getParameter<double>("errorReductionParameter");
+ realProperties["hcsits_relaxationParameter"] = solverHCSITSConf.getParameter<double>("relaxationParameter");
+ const Config::BlockHandle solverDEMConf = solverConf.getBlock("DEM");
+ realProperties["dem_collisionTimeNonDim"] = solverDEMConf.getParameter<double>("collisionTime") / solverConf.getParameter<double>("dt");
+ realProperties["dem_poissonsRatio"] = solverDEMConf.getParameter<double>("poissonsRatio");
+
+ const Config::BlockHandle distributionConf = config.getBlock("Distribution");
+ integerProperties["distribution_randomSeed"] = distributionConf.getParameter<int64_t>("randomSeed");
+ const Config::BlockHandle uniformConf = distributionConf.getBlock("Uniform");
+ realProperties["distribution_uniform_diameter"] = uniformConf.getParameter<real_t>("diameter");
+ const Config::BlockHandle logNormalConf = distributionConf.getBlock("LogNormal");
+ realProperties["distribution_logNormal_mu"] = logNormalConf.getParameter<real_t>("mu");
+ realProperties["distribution_logNormal_variance"] = logNormalConf.getParameter<real_t>("variance");
+ const Config::BlockHandle diamMassFracsConf = distributionConf.getBlock("DiameterMassFractions");
+ stringProperties["distribution_diamMassFracs_diameters"] = diamMassFracsConf.getParameter<std::string>("diameters");
+ stringProperties["distribution_diamMassFracs_massFractions"] = diamMassFracsConf.getParameter<std::string>("massFractions");
+ const Config::BlockHandle sievingConf = distributionConf.getBlock("SievingCurve");
+ stringProperties["distribution_sievingCurve_sieveSizes"] = sievingConf.getParameter<std::string>("sieveSizes");
+ stringProperties["distribution_sievingCurve_massFractions"] = sievingConf.getParameter<std::string>("massFractions");
+ integerProperties["distribution_sievingCurve_useDiscreteForm"] = (sievingConf.getParameter<bool>("useDiscreteForm")) ? 1 : 0;
+
+ const Config::BlockHandle shapeConf = config.getBlock("Shape");
+ stringProperties["shape_scaleMode"] = shapeConf.getParameter<std::string>("scaleMode");
+
+ const Config::BlockHandle ellipsoidConf = shapeConf.getBlock("Ellipsoid");
+ Vec3 ellipsoid_semiAxes = ellipsoidConf.getParameter< Vec3 >("semiAxes");
+ realProperties["shape_ellipsoid_semiAxis0"] = double(ellipsoid_semiAxes[0]);
+ realProperties["shape_ellipsoid_semiAxis1"] = double(ellipsoid_semiAxes[1]);
+ realProperties["shape_ellipsoid_semiAxis2"] = double(ellipsoid_semiAxes[2]);
+ const Config::BlockHandle equivalentEllipsoidConf = shapeConf.getBlock("EquivalentEllipsoid");
+ stringProperties["shape_equivalentEllipsoid_path"] = equivalentEllipsoidConf.getParameter<std::string>("path");
+ const Config::BlockHandle ellipsoidDistributionConf = shapeConf.getBlock("EllipsoidFormDistribution");
+ realProperties["shape_ellipsoidFromDistribution_elongationMean"] = ellipsoidDistributionConf.getParameter<real_t>("elongationMean");
+ realProperties["shape_ellipsoidFromDistribution_elongationStdDev"] = ellipsoidDistributionConf.getParameter<real_t>("elongationStdDev");
+ realProperties["shape_ellipsoidFromDistribution_flatnessMean"] = ellipsoidDistributionConf.getParameter<real_t>("flatnessMean");
+ realProperties["shape_ellipsoidFromDistribution_flatnessStdDev"] = ellipsoidDistributionConf.getParameter<real_t>("flatnessStdDev");
+ const Config::BlockHandle meshConf = shapeConf.getBlock("Mesh");
+ stringProperties["shape_mesh_path"] = meshConf.getParameter<std::string>("path");
+ const Config::BlockHandle meshDistributionConf = shapeConf.getBlock("MeshFormDistribution");
+ stringProperties["shape_meshFromDistribution_path"] = meshDistributionConf.getParameter<std::string>("path");
+ realProperties["shape_meshFromDistribution_elongationMean"] = meshDistributionConf.getParameter<real_t>("elongationMean");
+ realProperties["shape_meshFromDistribution_elongationStdDev"] = meshDistributionConf.getParameter<real_t>("elongationStdDev");
+ realProperties["shape_meshFromDistribution_flatnessMean"] = meshDistributionConf.getParameter<real_t>("flatnessMean");
+ realProperties["shape_meshFromDistribution_flatnessStdDev"] = meshDistributionConf.getParameter<real_t>("flatnessStdDev");
+ const Config::BlockHandle meshesUnscaledConf = shapeConf.getBlock("UnscaledMeshesPerFraction");
+ stringProperties["shape_meshesUnscaled_folder"] = meshesUnscaledConf.getParameter<std::string>("folder");
+
+ const Config::BlockHandle evaluationConf = config.getBlock("evaluation");
+ stringProperties["evaluation_histogramBins"] = evaluationConf.getParameter<std::string>("histogramBins");
+ realProperties["evaluation_layerHeight"] = evaluationConf.getParameter<real_t>("layerHeight");
+
+ integerProperties["shaking"] = (mainConf.getParameter<bool>("shaking")) ? 1 : 0;
+ const Config::BlockHandle shakingConf = config.getBlock("Shaking");
+ realProperties["shaking_amplitude"] = shakingConf.getParameter<double>("amplitude");
+ realProperties["shaking_period"] = shakingConf.getParameter<double>("period");
+ realProperties["shaking_duration"] = shakingConf.getParameter<double>("duration");
+ integerProperties["shaking_activeFromBeginning"] = (shakingConf.getParameter<bool>("activeFromBeginning")) ? 1 : 0;
+
+}
+
+class SelectTensorGlyphForEllipsoids
+{
+public:
+ using return_type = vtk::TensorGlyph;
+ return_type operator()(const data::Particle& p) const {
+ WALBERLA_CHECK_EQUAL(p.getBaseShape()->getShapeType(), data::Ellipsoid::SHAPE_TYPE);
+
+ auto ellipsoid = static_cast<data::Ellipsoid*>(p.getBaseShape().get());
+ return vtk::createTensorGlyph(ellipsoid->getSemiAxes(),p.getRotation());
+ }
+};
+
+
+
+/*
+ * Application to generate dense random packings of particles
+ *
+ * Main features:
+ * - two domain setups: cylindrical container, horizontally periodic
+ * - two simulation approaches: discrete element method (DEM), hard-contact semi-implicit timestepping solver (HCSITS)
+ * - different size distributions
+ * - different shapes: spherical, ellipsoidal, polygonal as given by mesh
+ * - evaluation of vertical porosity profile
+ * - VTK visualization
+ * - logging of final result and all properties into SQlite database
+ * - requires OpenMesh
+ *
+ * Simulation process:
+ * - Generation phase: continuous generation in upper part of domain and settling due to gravity
+ * - Shaking phase (optional, can also be active during generation phase): Shaking in a horizontal direction to compactify packing
+ * - Termination phase: Run until converged state is reached
+ *
+ * See corresponding publication by C. Rettinger for more infos
+ *
+ */
+int main(int argc, char **argv) {
+
+ /// Setup
+ Environment env(argc, argv);
+
+ /// Config
+ auto cfg = env.config();
+ WALBERLA_LOG_INFO_ON_ROOT(*cfg);
+ if (cfg == nullptr) WALBERLA_ABORT("No config specified!");
+ const Config::BlockHandle mainConf = cfg->getBlock("ParticlePacking");
+
+ std::string domainSetup = mainConf.getParameter<std::string>("domainSetup");
+ WALBERLA_CHECK(domainSetup == "container" || domainSetup == "periodic");
+ real_t domainWidth = mainConf.getParameter<real_t>("domainWidth");
+ real_t domainHeight = mainConf.getParameter<real_t>("domainHeight");
+ real_t particleDensity = mainConf.getParameter<real_t>("particleDensity");
+ real_t ambientDensity = mainConf.getParameter<real_t>("ambientDensity");
+ real_t gravitationalAcceleration = mainConf.getParameter<real_t>("gravitationalAcceleration");
+ real_t reducedGravitationalAcceleration = (particleDensity - ambientDensity) / particleDensity * gravitationalAcceleration;
+
+ std::string particleDistribution = mainConf.getParameter<std::string>("particleDistribution");
+ std::string particleShape = mainConf.getParameter<std::string>("particleShape");
+ real_t limitVelocity = mainConf.getParameter<real_t>("limitVelocity");
+ real_t initialVelocity = mainConf.getParameter<real_t>("initialVelocity");
+ real_t initialGenerationHeightRatioStart = mainConf.getParameter<real_t>("initialGenerationHeightRatioStart");
+ real_t initialGenerationHeightRatioEnd = mainConf.getParameter<real_t>("initialGenerationHeightRatioEnd");
+ real_t generationSpacing = mainConf.getParameter<real_t>("generationSpacing");
+ WALBERLA_CHECK_GREATER(domainWidth, generationSpacing, "Generation Spacing has to be smaller than domain size");
+ real_t generationHeightRatioStart = mainConf.getParameter<real_t>("generationHeightRatioStart");
+ real_t generationHeightRatioEnd = mainConf.getParameter<real_t>("generationHeightRatioEnd");
+ bool scaleGenerationSpacingWithForm = mainConf.getParameter<bool>("scaleGenerationSpacingWithForm");
+ real_t totalParticleMass = mainConf.getParameter<real_t>("totalParticleMass");
+
+ real_t visSpacingInSeconds = mainConf.getParameter<real_t>("visSpacing");
+ real_t infoSpacingInSeconds = mainConf.getParameter<real_t>("infoSpacing");
+ real_t loggingSpacingInSeconds = mainConf.getParameter<real_t>("loggingSpacing");
+ Vector3<uint_t> numBlocksPerDirection = mainConf.getParameter< Vector3<uint_t> >("numBlocksPerDirection");
+ real_t terminalVelocity = mainConf.getParameter<real_t>("terminalVelocity");
+ real_t terminalRelativeHeightChange = mainConf.getParameter<real_t>("terminalRelativeHeightChange");
+ real_t terminationCheckingSpacing = mainConf.getParameter<real_t>("terminationCheckingSpacing");
+ real_t minimalTerminalRunTime = mainConf.getParameter<real_t>("minimalTerminalRunTime");
+
+ bool useHashGrids = mainConf.getParameter<bool>("useHashGrids");
+
+ std::string solver = mainConf.getParameter<std::string>("solver");
+
+ int particleSortingSpacing = mainConf.getParameter< int >("particleSortingSpacing");
+
+
+ const Config::BlockHandle solverConf = cfg->getBlock("Solver");
+ real_t dt = solverConf.getParameter<real_t>("dt");
+ real_t frictionCoefficient = solverConf.getParameter<real_t>("frictionCoefficient");
+ real_t coefficientOfRestitution = solverConf.getParameter<real_t>("coefficientOfRestitution");
+ real_t velocityDampingCoefficient = solverConf.getParameter<real_t>("velocityDampingCoefficient");
+
+ uint_t visSpacing = uint_c(visSpacingInSeconds / dt);
+ uint_t infoSpacing = uint_c(infoSpacingInSeconds / dt);
+ uint_t loggingSpacing = uint_c(loggingSpacingInSeconds / dt);
+ WALBERLA_LOG_INFO_ON_ROOT("VTK spacing = " << visSpacing << ", info spacing = " << infoSpacing << ", logging spacing = " << loggingSpacing);
+
+ const Config::BlockHandle solverHCSITSConf = solverConf.getBlock("HCSITS");
+ real_t hcsits_errorReductionParameter = solverHCSITSConf.getParameter<real_t>("errorReductionParameter");
+ real_t hcsits_relaxationParameter = solverHCSITSConf.getParameter<real_t>("relaxationParameter");
+ std::string hcsits_relaxationModel = solverHCSITSConf.getParameter<std::string>("relaxationModel");
+ uint_t hcsits_numberOfIterations = solverHCSITSConf.getParameter<uint_t>("numberOfIterations");
+
+ const Config::BlockHandle solverDEMConf = solverConf.getBlock("DEM");
+ real_t dem_collisionTime = solverDEMConf.getParameter<real_t>("collisionTime");
+ //real_t dem_stiffnessNormal = solverDEMConf.getParameter<real_t>("stiffnessNormal");
+ real_t dem_poissonsRatio = solverDEMConf.getParameter<real_t>("poissonsRatio");
+ real_t dem_kappa = real_t(2) * ( real_t(1) - dem_poissonsRatio ) / ( real_t(2) - dem_poissonsRatio ) ; // from Thornton et al
+
+ bool shaking = mainConf.getParameter<bool>("shaking");
+ const Config::BlockHandle shakingConf = cfg->getBlock("Shaking");
+ real_t shaking_amplitude = shakingConf.getParameter<real_t>("amplitude");
+ real_t shaking_period = shakingConf.getParameter<real_t>("period");
+ real_t shaking_duration = shakingConf.getParameter<real_t>("duration");
+ bool shaking_activeFromBeginning = shakingConf.getParameter<bool>("activeFromBeginning");
+
+ const Config::BlockHandle evaluationConf = cfg->getBlock("evaluation");
+ auto evaluationHistogramBins = parseStringToVector<real_t>(evaluationConf.getParameter<std::string>("histogramBins"));
+ //real_t voxelsPerMm = evaluationConf.getParameter<real_t>("voxelsPerMm");
+ std::string porosityProfileFolder = evaluationConf.getParameter<std::string>("porosityProfileFolder");
+ real_t evaluationLayerHeight = evaluationConf.getParameter<real_t>("layerHeight");
+ std::string vtkOutputFolder = evaluationConf.getParameter<std::string>("vtkFolder");
+ std::string vtkFinalFolder = evaluationConf.getParameter<std::string>("vtkFinalFolder");
+ std::string sqlDBFileName = evaluationConf.getParameter<std::string>("sqlDBFileName");
+
+ const Config::BlockHandle shapeConf = cfg->getBlock("Shape");
+ ScaleMode shapeScaleMode = str_to_scaleMode(shapeConf.getParameter<std::string>("scaleMode"));
+ const Config::BlockHandle distributionConf = cfg->getBlock("Distribution");
+
+ /// BlockForest
+ math::AABB simulationDomain(-0.5_r*domainWidth, -0.5_r*domainWidth, 0_r,
+ 0.5_r*domainWidth, 0.5_r*domainWidth, domainHeight);
+ Vector3<bool> isPeriodic = (domainSetup == "container") ? Vector3<bool>(false) : Vector3<bool>(true, true, false);
+
+ WALBERLA_LOG_INFO_ON_ROOT("Creating domain of size " << simulationDomain);
+
+ shared_ptr<BlockForest> forest = blockforest::createBlockForest(simulationDomain, numBlocksPerDirection, isPeriodic);
+ auto domain = std::make_shared<mesa_pd::domain::BlockForestDomain>(forest);
+
+ /// MESAPD Data
+ auto particleStorage = std::make_shared<data::ParticleStorage>(1);
+ auto contactStorage = std::make_shared<data::ContactStorage>(1);
+ ParticleAccessorWithShape particleAccessor(particleStorage);
+ data::ContactAccessor contactAccessor(contactStorage);
+
+ // configure shape creation
+ shared_ptr<ShapeGenerator> shapeGenerator;
+ if(particleShape == "Sphere")
+ {
+ shapeGenerator = make_shared<SphereGenerator>();
+ } else if(particleShape == "Ellipsoid")
+ {
+ auto ellipsoidConfig = shapeConf.getBlock("Ellipsoid");
+ std::vector<Vec3> semiAxes = {ellipsoidConfig.getParameter<Vec3>("semiAxes")};
+
+ shared_ptr<NormalizedFormGenerator> normalizedFormGenerator = make_shared<SampleFormGenerator>(semiAxes, shapeScaleMode);
+ shapeGenerator = make_shared<EllipsoidGenerator>(normalizedFormGenerator);
+ } else if(particleShape == "EquivalentEllipsoid")
+ {
+ auto ellipsoidConfig = shapeConf.getBlock("EquivalentEllipsoid");
+ std::string meshPath = ellipsoidConfig.getParameter<std::string>("path");
+
+ auto meshFileNames = getMeshFilesFromPath(meshPath);
+ auto semiAxes = extractSemiAxesFromMeshFiles(meshFileNames);
+ shared_ptr<NormalizedFormGenerator> normalizedFormGenerator = make_shared<SampleFormGenerator>(semiAxes, shapeScaleMode);
+
+ shapeGenerator = make_shared<EllipsoidGenerator>(normalizedFormGenerator);
+ } else if(particleShape == "EllipsoidFormDistribution")
+ {
+ auto ellipsoidConfig = shapeConf.getBlock("EllipsoidFormDistribution");
+ auto elongationMean = ellipsoidConfig.getParameter<real_t>("elongationMean");
+ auto elongationStdDev = ellipsoidConfig.getParameter<real_t>("elongationStdDev");
+ auto flatnessMean = ellipsoidConfig.getParameter<real_t>("flatnessMean");
+ auto flatnessStdDev = ellipsoidConfig.getParameter<real_t>("flatnessStdDev");
+
+ shared_ptr<NormalizedFormGenerator> normalizedFormGenerator = make_shared<DistributionFormGenerator>(elongationMean, elongationStdDev, flatnessMean, flatnessStdDev, shapeScaleMode);
+
+ shapeGenerator = make_shared<EllipsoidGenerator>(normalizedFormGenerator);
+ }
+ else if(particleShape == "Mesh")
+ {
+ auto meshConfig = shapeConf.getBlock("Mesh");
+ std::string meshPath = meshConfig.getParameter<std::string>("path");
+
+ auto meshFileNames = getMeshFilesFromPath(meshPath);
+ shared_ptr<NormalizedFormGenerator> normalizedFormGenerator = make_shared<ConstFormGenerator>();
+
+ shapeGenerator = make_shared<MeshesGenerator>(meshFileNames, shapeScaleMode, normalizedFormGenerator);
+ } else if(particleShape == "MeshFormDistribution")
+ {
+ auto meshConfig = shapeConf.getBlock("MeshFormDistribution");
+ std::string meshPath = meshConfig.getParameter<std::string>("path");
+ auto elongationMean = meshConfig.getParameter<real_t>("elongationMean");
+ auto elongationStdDev = meshConfig.getParameter<real_t>("elongationStdDev");
+ auto flatnessMean = meshConfig.getParameter<real_t>("flatnessMean");
+ auto flatnessStdDev = meshConfig.getParameter<real_t>("flatnessStdDev");
+
+ auto meshFileNames = getMeshFilesFromPath(meshPath);
+ shared_ptr<NormalizedFormGenerator> normalizedFormGenerator = make_shared<DistributionFormGenerator>(elongationMean, elongationStdDev, flatnessMean, flatnessStdDev, shapeScaleMode);
+
+ shapeGenerator = make_shared<MeshesGenerator>(meshFileNames, shapeScaleMode, normalizedFormGenerator);
+ } else if(particleShape == "UnscaledMeshesPerFraction")
+ {
+ shapeGenerator = make_shared<UnscaledMeshesPerFractionGenerator>(shapeConf, parseStringToVector<real_t>(distributionConf.getBlock("SievingCurve").getParameter<std::string>("massFractions")));
+ } else
+ {
+ WALBERLA_ABORT("Unknown shape " << particleShape);
+ }
+ WALBERLA_LOG_INFO_ON_ROOT("Will create particles with ");
+ WALBERLA_LOG_INFO_ON_ROOT(" - maximum diameter scaling of " << shapeGenerator->getMaxDiameterScalingFactor());
+ WALBERLA_LOG_INFO_ON_ROOT(" - normal volume " << shapeGenerator->getNormalVolume());
+ WALBERLA_LOG_INFO_ON_ROOT(" - " << (shapeGenerator->generatesSingleShape() ? "single shape" : "multiple shapes"));
+
+ // configure size creation
+ int randomSeedFromConfig = distributionConf.getParameter<int>("randomSeed");
+ uint_t randomSeed = (randomSeedFromConfig >= 0) ? uint_c(randomSeedFromConfig) : uint_c(time(nullptr));
+ WALBERLA_LOG_INFO_ON_ROOT("Random seed of " << randomSeed);
+
+ shared_ptr<DiameterGenerator> diameterGenerator;
+ real_t minGenerationParticleDiameter = real_t(0);
+ real_t maxGenerationParticleDiameter = std::numeric_limits<real_t>::max();
+
+ if(particleDistribution == "LogNormal")
+ {
+ const Config::BlockHandle logNormalConf = distributionConf.getBlock("LogNormal");
+ real_t mu = logNormalConf.getParameter<real_t>("mu");
+ real_t variance = logNormalConf.getParameter<real_t>("variance");
+ diameterGenerator = make_shared<LogNormal>(mu, variance, randomSeed);
+ // min and max diameter not determinable
+ WALBERLA_LOG_INFO_ON_ROOT("Using log-normal distribution with mu = " << mu << ", var = " << variance);
+ }
+ else if(particleDistribution == "Uniform")
+ {
+ const Config::BlockHandle uniformConf = distributionConf.getBlock("Uniform");
+ real_t diameter = uniformConf.getParameter<real_t>("diameter");
+ diameterGenerator = make_shared<Uniform>(diameter);
+ minGenerationParticleDiameter = diameter;
+ maxGenerationParticleDiameter = diameter;
+ WALBERLA_LOG_INFO_ON_ROOT("Using uniform distribution");
+ }
+ else if(particleDistribution == "DiameterMassFractions")
+ {
+ const Config::BlockHandle sievingConf = distributionConf.getBlock("DiameterMassFractions");
+ auto diameters = parseStringToVector<real_t>(sievingConf.getParameter<std::string>("diameters"));
+ auto massFractions = parseStringToVector<real_t>(sievingConf.getParameter<std::string>("massFractions"));
+ diameterGenerator = make_shared<DiscreteSieving>(diameters, massFractions, randomSeed, shapeGenerator->getNormalVolume(), totalParticleMass, particleDensity);
+
+ maxGenerationParticleDiameter = real_t(0);
+ minGenerationParticleDiameter = std::numeric_limits<real_t>::max();
+ for(uint_t i = 0; i < diameters.size(); ++i) {
+ if(massFractions[i] > real_t(0)) {
+ maxGenerationParticleDiameter = std::max(maxGenerationParticleDiameter, diameters[i]);
+ minGenerationParticleDiameter = std::min(minGenerationParticleDiameter, diameters[i]);
+ }
+ }
+ WALBERLA_LOG_INFO_ON_ROOT("Using diameter - mass fraction distribution");
+ }
+ else if(particleDistribution == "SievingCurve")
+ {
+ const Config::BlockHandle sievingConf = distributionConf.getBlock("SievingCurve");
+ auto sieveSizes = parseStringToVector<real_t>(sievingConf.getParameter<std::string>("sieveSizes"));
+ auto massFractions = parseStringToVector<real_t>(sievingConf.getParameter<std::string>("massFractions"));
+ bool useDiscreteForm = sievingConf.getParameter<bool>("useDiscreteForm");
+
+ auto diameters = getMeanDiametersFromSieveSizes(sieveSizes);
+ real_t d50 = computePercentileFromSieveDistribution(diameters, massFractions, 50_r);
+ real_t d16 = computePercentileFromSieveDistribution(diameters, massFractions, 16_r);
+ real_t d84 = computePercentileFromSieveDistribution(diameters, massFractions, 84_r);
+ real_t stdDev = std::sqrt(d84/d16);
+ WALBERLA_LOG_INFO_ON_ROOT("Curve properties: D50 = " << d50 << ", D16 = " << d16 << ", D84 = " << d84 << ", estimated std. dev. = " << stdDev);
+
+ maxGenerationParticleDiameter = real_t(0);
+ minGenerationParticleDiameter = std::numeric_limits<real_t>::max();
+ if(useDiscreteForm)
+ {
+ diameterGenerator = make_shared<DiscreteSieving>(diameters, massFractions, randomSeed, shapeGenerator->getNormalVolume(), totalParticleMass, particleDensity);
+ for(uint_t i = 0; i < diameters.size(); ++i) {
+ if(massFractions[i] > real_t(0)) {
+ maxGenerationParticleDiameter = std::max(maxGenerationParticleDiameter, diameters[i]);
+ minGenerationParticleDiameter = std::min(minGenerationParticleDiameter, diameters[i]);
+ }
+ }
+ WALBERLA_LOG_INFO_ON_ROOT("Using discrete sieving curve distribution");
+
+ } else {
+ diameterGenerator = make_shared<ContinuousSieving>(sieveSizes, massFractions, randomSeed, shapeGenerator->getNormalVolume(), totalParticleMass, particleDensity);
+ for(uint_t i = 0; i < sieveSizes.size()-1; ++i) {
+ if(massFractions[i] > real_t(0)) {
+ maxGenerationParticleDiameter = std::max(maxGenerationParticleDiameter, std::max(sieveSizes[i],sieveSizes[i+1]));
+ minGenerationParticleDiameter = std::min(minGenerationParticleDiameter, std::min(sieveSizes[i],sieveSizes[i+1]));
+ }
+ }
+ WALBERLA_LOG_INFO_ON_ROOT("Using piece-wise constant / continuous sieving curve distribution");
+ }
+ }
+ else
+ {
+ WALBERLA_ABORT("Unknown particle distribution specified: " << particleDistribution);
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("Generate with diameters in range [" << minGenerationParticleDiameter << ", " << maxGenerationParticleDiameter << "] and generation spacing = " << generationSpacing);
+
+ bool useOpenMP = false;
+
+ real_t smallestBlockSize = std::min(simulationDomain.xSize() / real_c(numBlocksPerDirection[0]),
+ std::min(simulationDomain.ySize() / real_c(numBlocksPerDirection[1]),
+ simulationDomain.zSize() / real_c(numBlocksPerDirection[2])));
+
+ /*
+ auto geometricMeanDiameter = std::sqrt(minGenerationParticleDiameter * maxGenerationParticleDiameter);
+ if(solver=="DEM")
+ {
+ auto getCollisionTime = [dem_stiffnessNormal, coefficientOfRestitution, particleDensity](real_t d){
+ auto meanEffMass = d * d * d * math::pi * particleDensity/ (6_r*2_r);
+ return std::sqrt((std::pow(std::log(coefficientOfRestitution),2) + math::pi * math::pi ) / (dem_stiffnessNormal / meanEffMass));};
+ WALBERLA_LOG_INFO_ON_ROOT("DEM parameterization with kn = " << dem_stiffnessNormal << " and cor = " << coefficientOfRestitution
+ << " expects collision time in range [" << getCollisionTime(minGenerationParticleDiameter) << ", " << getCollisionTime(maxGenerationParticleDiameter) << "]");
+ }
+ */
+
+ // plane at top and bottom
+ createPlane(particleStorage, Vector3<real_t>(0), Vector3<real_t>(0, 0, 1));
+ createPlane(particleStorage, Vector3<real_t>(0_r,0_r,simulationDomain.zMax()), Vector3<real_t>(0, 0, -1));
+
+ real_t domainVolume = simulationDomain.volume();
+ if(domainSetup == "container")
+ {
+ createCylindricalBoundary(particleStorage, Vector3<real_t>(0), Vector3<real_t>(0, 0, 1), 0.5_r*domainWidth);
+ domainVolume = math::pi * domainWidth * domainWidth * 0.25_r * simulationDomain.zSize();
+ }
+
+ real_t maximumAllowedInteractionRadius = std::numeric_limits<real_t>::infinity();
+ if(domainSetup == "periodic")
+ {
+ // avoid that two large particles are next to each other and would, due to periodic mapping, have 2 different contact points with each other |( p1 () p2 ()| p1 )
+ maximumAllowedInteractionRadius = 0.25_r * domainWidth; // max diameter = domainWidth / 2
+ WALBERLA_LOG_INFO_ON_ROOT("Periodic case: the maximum interaction radius is restricted to " << maximumAllowedInteractionRadius << " to ensure valid periodic interaction" );
+ if(numBlocksPerDirection[0] < 3 || numBlocksPerDirection[1] < 3) WALBERLA_LOG_INFO_ON_ROOT("Warning: At least 3 blocks per periodic direction required for proper simulation!")
+ }
+
+ // fill domain with particles initially
+ real_t maxGenerationHeight = simulationDomain.zMax() - generationSpacing;
+ real_t minGenerationHeight = generationSpacing;
+ ParticleCreator particleCreator(particleStorage, domain, simulationDomain, domainSetup, particleDensity, scaleGenerationSpacingWithForm);
+ particleCreator.createParticles(std::max(minGenerationHeight, initialGenerationHeightRatioStart * simulationDomain.zMax()),
+ std::min(maxGenerationHeight, initialGenerationHeightRatioEnd * simulationDomain.zMax()),
+ generationSpacing, diameterGenerator, shapeGenerator, initialVelocity, maximumAllowedInteractionRadius );
+
+ 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.mpiAllGather();
+ WALBERLA_LOG_INFO_ON_ROOT("Statistics of initially created particles' interaction diameters: " << diameterSample.format());
+
+ real_t maxParticleDiameter = maxGenerationParticleDiameter * shapeGenerator->getMaxDiameterScalingFactor();
+ if(maxParticleDiameter < diameterSample.max())
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Maximum interaction diameter from samples is larger than estimated maximum diameter, will use sampled one instead.")
+ maxParticleDiameter = 1.1_r * diameterSample.max(); // 10% safety margin
+ }
+ if(maxParticleDiameter > 2_r * maximumAllowedInteractionRadius)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Warning: Maximum expected particle interaction diameter ("<< maxParticleDiameter << ") is larger than maximum allowed interaction diameter - check that the generated size & form distributions match the expected ones!");
+ maxParticleDiameter = 2_r * maximumAllowedInteractionRadius;
+ }
+
+ bool useNextNeighborSync = 2_r * smallestBlockSize > maxParticleDiameter;
+
+ WALBERLA_LOG_INFO_ON_ROOT("Sync info: maximum expected interaction diameter = " << maxParticleDiameter << " and smallest block size = " << smallestBlockSize);
+
+ // sync functionality
+ kernel::AssocToBlock associateToBlock(forest);
+ std::function<void(void)> syncCall;
+ if(useNextNeighborSync)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Using next neighbor sync!");
+ syncCall = [&particleStorage,&forest,&domain](){
+ mpi::SyncNextNeighborsBlockForest syncParticles;
+ syncParticles(*particleStorage, forest, domain);
+ };
+ } else {
+ WALBERLA_LOG_INFO_ON_ROOT("Using ghost owner sync!");
+ syncCall = [&particleStorage,&domain](){
+ mpi::SyncGhostOwners syncGhostOwnersFunc;
+ syncGhostOwnersFunc(*particleStorage, *domain);
+ };
+ }
+
+ // initial sync
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor, associateToBlock, particleAccessor);
+ if(useNextNeighborSync){ syncCall(); }
+ else { for(uint_t i = 0; i < uint_c(std::ceil(maxParticleDiameter/smallestBlockSize)); ++i) syncCall(); }
+
+
+ // create linked cells data structure
+ real_t linkedCellWidth = 1.01_r * maxParticleDiameter;
+ WALBERLA_LOG_INFO_ON_ROOT("Using linked cells with cell width = " << linkedCellWidth);
+ data::LinkedCells linkedCells(domain->getUnionOfLocalAABBs().getExtended(linkedCellWidth), linkedCellWidth );
+
+ {
+ auto info = evaluateParticleInfo(particleAccessor);
+ WALBERLA_LOG_INFO_ON_ROOT(info);
+ }
+
+ /// VTK Output
+ if(visSpacing > 0)
+ {
+ auto vtkDomainOutput = walberla::vtk::createVTKOutput_DomainDecomposition( forest, "domain_decomposition", 1, vtkOutputFolder, "simulation_step" );
+ vtkDomainOutput->write();
+ }
+
+ // mesapd particle output
+ auto particleVtkOutput = make_shared<vtk::ParticleVtkOutput>(particleStorage);
+ particleVtkOutput->addOutput<data::SelectParticleUid>("uid");
+ particleVtkOutput->addOutput<data::SelectParticleOwner>("owner");
+ particleVtkOutput->addOutput<data::SelectParticleInteractionRadius>("interactionRadius");
+ if(particleShape.find("Ellipsoid") != std::string::npos)
+ {
+ particleVtkOutput->addOutput<SelectTensorGlyphForEllipsoids>("tensorGlyph");
+ }
+ particleVtkOutput->addOutput<data::SelectParticleLinearVelocity>("velocity");
+ particleVtkOutput->addOutput<data::SelectParticleNumContacts>("numContacts");
+ auto vtkParticleSelector = [](const mesa_pd::data::ParticleStorage::iterator &pIt) {
+ return (pIt->getBaseShape()->getShapeType() != data::HalfSpace::SHAPE_TYPE &&
+ pIt->getBaseShape()->getShapeType() != data::CylindricalBoundary::SHAPE_TYPE &&
+ !isSet(pIt->getFlags(), data::particle_flags::GHOST));
+ };
+ particleVtkOutput->setParticleSelector(vtkParticleSelector);
+ auto particleVtkWriter = walberla::vtk::createVTKOutput_PointData(particleVtkOutput, "particles", visSpacing, vtkOutputFolder, "simulation_step");
+
+ mesa_pd::MeshParticleVTKOutput< mesh::PolyMesh > meshParticleVTK(particleStorage, "mesh", visSpacing, vtkOutputFolder);
+ meshParticleVTK.addFaceOutput< data::SelectParticleUid >("UID");
+ meshParticleVTK.addVertexOutput< data::SelectParticleInteractionRadius >("InteractionRadius");
+ 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);
+ meshParticleVTK.setParticleSelector(vtkParticleSelector);
+ meshParticleVTK.addVertexDataSource(surfaceVelDataSource);
+
+
+ /// MESAPD kernels
+
+ //collision detection
+ data::HashGrids hashGrids;
+ kernel::InsertParticleIntoLinkedCells initializeLinkedCells;
+ kernel::DetectAndStoreContacts detectAndStore(*contactStorage);
+
+ //DEM
+ kernel::SemiImplicitEuler dem_integration( dt );
+ kernel::LinearSpringDashpot dem_collision(1);
+ dem_collision.setFrictionCoefficientStatic(0,0,0_r); // no static friction
+ dem_collision.setFrictionCoefficientDynamic(0,0, frictionCoefficient);
+ // stiffness and damping depend on effective mass -> have to be calculated for each collision individually
+
+ //HCSITS
+ kernel::InitContactsForHCSITS hcsits_initContacts(1);
+ hcsits_initContacts.setFriction(0,0,frictionCoefficient);
+ hcsits_initContacts.setErp(hcsits_errorReductionParameter);
+ kernel::InitParticlesForHCSITS hcsits_initParticles;
+ hcsits_initParticles.setGlobalAcceleration(Vector3<real_t>(0,0,-reducedGravitationalAcceleration));
+ kernel::HCSITSRelaxationStep hcsits_relaxationStep;
+ hcsits_relaxationStep.setRelaxationModel(relaxationModelFromString(hcsits_relaxationModel));
+ hcsits_relaxationStep.setCor(coefficientOfRestitution); // Only effective for PGSM
+ kernel::IntegrateParticlesHCSITS hcsits_integration;
+
+ // sync
+ mpi::ReduceContactHistory reduceAndSwapContactHistory;
+ mpi::BroadcastProperty broadcastKernel;
+ mpi::ReduceProperty reductionKernel;
+
+ uint_t timestep = 0;
+ WALBERLA_LOG_INFO_ON_ROOT("Starting simulation in domain of volume " << domainVolume << " m^3.");
+ WALBERLA_LOG_INFO_ON_ROOT("Will terminate generation when particle mass is above " << totalParticleMass << " kg.");
+
+ real_t velocityDampingFactor = std::pow(velocityDampingCoefficient, dt);
+ WALBERLA_LOG_INFO_ON_ROOT("Once all particles are created, will apply velocity damping of " << velocityDampingFactor << " per time step.");
+ real_t oldAvgParticleHeight = real_t(1);
+ real_t oldMaxParticleHeight = real_t(1);
+ real_t timeLastTerminationCheck = real_t(0);
+ real_t timeLastCreation = real_t(0);
+ real_t maximumTimeBetweenCreation = (generationHeightRatioEnd-generationHeightRatioStart)*simulationDomain.zSize() / initialVelocity; // = time, particles need at max to clear/pass the creation domain
+ WALBERLA_LOG_INFO_ON_ROOT("Maximum time between creation steps: " << maximumTimeBetweenCreation);
+
+ bool isShakingActive = false;
+ real_t timeBeginShaking = real_t(-1);
+ if(shaking && shaking_activeFromBeginning)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Will use shaking from beginning.");
+ isShakingActive = true;
+ timeBeginShaking = real_t(0);
+ }
+ real_t timeEndShaking = real_t(-1);
+ real_t timeBeginDamping = real_t(-1);
+
+ if(limitVelocity > 0_r) WALBERLA_LOG_INFO_ON_ROOT("Will apply limiting of translational particle velocity to maximal magnitude of " << limitVelocity);
+
+ std::string uniqueFileIdentifier = std::to_string(std::chrono::system_clock::now().time_since_epoch().count()); // used as hash to identify this run
+ walberla::mpi::broadcastObject(uniqueFileIdentifier);
+
+ SizeEvaluator particleSizeEvaluator(shapeScaleMode);
+ std::vector<std::tuple<std::string, std::function<real_t(Vec3)>>> particleShapeEvaluators = {std::make_tuple("flatness",getFlatnessFromSemiAxes),
+ std::make_tuple("elongation",getElongationFromSemiAxes),
+ std::make_tuple("equancy",getEquancyFromSemiAxes)};
+ uint_t numShapeBins = 17;
+ std::vector<std::vector<real_t>> particleShapeBins(particleShapeEvaluators.size());
+ for(auto& pSB: particleShapeBins)
+ {
+ pSB = std::vector<real_t>(numShapeBins, real_t(0));
+ real_t binBegin = 0_r;
+ real_t binEnd = 1_r;
+ real_t inc = (binEnd - binBegin) / real_t(numShapeBins-1);
+ real_t val = binBegin;
+ for(uint_t i = 0; i < numShapeBins; ++i, val += inc) pSB[i] = val;
+ }
+
+ ParticleHistogram particleHistogram(evaluationHistogramBins, particleSizeEvaluator, particleShapeBins, particleShapeEvaluators);
+
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ particleHistogram, particleAccessor);
+ particleHistogram.evaluate();
+ WALBERLA_LOG_INFO_ON_ROOT(particleHistogram);
+
+
+ PorosityPerHorizontalLayerEvaluator porosityEvaluator(evaluationLayerHeight, simulationDomain, domainSetup);
+
+ std::string loggingFileName = porosityProfileFolder + "/" + uniqueFileIdentifier + "_logging.txt";
+ WALBERLA_LOG_INFO_ON_ROOT("Writing logging file to " << loggingFileName);
+ LoggingWriter loggingWriter(loggingFileName);
+
+
+
+ //particleStorage->forEachParticle(useOpenMP, kernel::SelectAll(), particleAccessor,
+ // [](const size_t idx, auto &ac){WALBERLA_LOG_INFO(ac.getUid(idx) << " " << ac.getPosition(idx) << " " <<!isSet(ac.getFlags(idx), data::particle_flags::GHOST) << "\n");}, particleAccessor);
+
+ //std::ofstream file;
+ //std::string fileName = "rank" + std::to_string(walberla::mpi::MPIManager::instance()->rank()) + ".txt";
+ //file.open( fileName.c_str() );
+
+ WcTimingTree timing;
+
+ timing.start("Simulation");
+
+ bool terminateSimulation = false;
+ while (!terminateSimulation) {
+
+ real_t currentTime = dt * real_c(timestep);
+
+ timing.start("Sorting");
+ if(particleSortingSpacing > 0 && timestep % uint_c(particleSortingSpacing) == 0 && !useHashGrids)
+ {
+ sorting::LinearizedCompareFunctor linearSorting(linkedCells.domain_, linkedCells.numCellsPerDim_);
+ particleStorage->sort(linearSorting);
+ }
+ timing.stop("Sorting");
+
+ timing.start("VTK");
+ if(particleShape.find("Mesh") != std::string::npos) meshParticleVTK(particleAccessor);
+ else particleVtkWriter->write();
+ timing.stop("VTK");
+
+
+ contactStorage->clear();
+
+ if(useHashGrids)
+ {
+ timing.start("Hash grid");
+ hashGrids.clearAll();
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectAll(), particleAccessor, hashGrids, particleAccessor);
+ timing.stop("Hash grid");
+
+ timing.start("Contact detection");
+ hashGrids.forEachParticlePairHalf(useOpenMP, kernel::ExcludeInfiniteInfinite(), particleAccessor,
+ [domain, contactStorage](size_t idx1, size_t idx2, ParticleAccessorWithShape &ac){
+ kernel::DoubleCast double_cast;
+ mpi::ContactFilter contact_filter;
+ collision_detection::GeneralContactDetection contactDetection;
+ //Attention: does not use contact threshold in general case (GJK)
+
+ if (double_cast(idx1, idx2, ac, contactDetection, ac)) {
+ if (contact_filter(contactDetection.getIdx1(), contactDetection.getIdx2(), ac, contactDetection.getContactPoint(), *domain)) {
+ auto c = contactStorage->create();
+ c->setId1(contactDetection.getIdx1());
+ c->setId2(contactDetection.getIdx2());
+ c->setDistance(contactDetection.getPenetrationDepth());
+ c->setNormal(contactDetection.getContactNormal());
+ c->setPosition(contactDetection.getContactPoint());
+ }
+ }
+ }, particleAccessor);
+ timing.stop("Contact detection");
+
+ } else
+ {
+ // use linked cells
+
+ timing.start("Linked cells");
+ linkedCells.clear();
+ particleStorage->forEachParticle(useOpenMP,kernel::SelectAll(),particleAccessor,
+ initializeLinkedCells, particleAccessor, linkedCells);
+ timing.stop("Linked cells");
+
+ timing.start("Contact detection");
+ if(particleShape == "Sphere")
+ {
+ collision_detection::AnalyticContactDetection contactDetection;
+ //acd.getContactThreshold() = contactThreshold;
+ linkedCells.forEachParticlePairHalf(useOpenMP, kernel::ExcludeInfiniteInfinite(), particleAccessor,
+ detectAndStore, particleAccessor, *domain, contactDetection);
+
+ } else
+ {
+ linkedCells.forEachParticlePairHalf(useOpenMP, kernel::ExcludeInfiniteInfinite(), particleAccessor,
+ [domain, contactStorage](size_t idx1, size_t idx2, ParticleAccessorWithShape &ac){
+
+ collision_detection::GeneralContactDetection contactDetection;
+ //Attention: does not use contact threshold in general case (GJK)
+
+ // coarse collision detection via interaction radii
+ data::Sphere sp1(ac.getInteractionRadius(idx1));
+ data::Sphere sp2(ac.getInteractionRadius(idx2));
+ if(contactDetection(idx1, idx2, sp1, sp2, ac)) {
+ //NOTE: this works also for infinite particles ( plane, cylindrical boundary) since contact_detection return true
+ // and the following contact_filter treats all local-global interactions independently of contact detection result, which would be non-sense for this interaction
+
+ mpi::ContactFilter contact_filter;
+ if (contact_filter(contactDetection.getIdx1(), contactDetection.getIdx2(), ac, contactDetection.getContactPoint(), *domain)) {
+ //NOTE: usually we first do fine collision detection and then the exact contact location determines the process that handles this contact
+ // however, along periodic boundaries, the GJK/EPA for meshes seems to be numerical unstable and yields (sometimes) different contact points for the same interaction pair (but periodically transformed)
+ // as a result, the same contact appears twice and potentially handled by two processes simultaneously.
+ // thus we change the ordering and do the contact filtering according to the result of the coarse collision detection, i.e. the bounding sphere check
+ kernel::DoubleCast double_cast;
+ if (double_cast(idx1, idx2, ac, contactDetection, ac)) {
+ auto c = contactStorage->create();
+ c->setId1(contactDetection.getIdx1());
+ c->setId2(contactDetection.getIdx2());
+ c->setDistance(contactDetection.getPenetrationDepth());
+ c->setNormal(contactDetection.getContactNormal());
+ c->setPosition(contactDetection.getContactPoint());
+ }
+ }
+ }}, particleAccessor);
+ }
+
+ timing.stop("Contact detection");
+ }
+
+ timing.start("Contact eval");
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectAll(), particleAccessor,
+ [](size_t p_idx, ParticleAccessorWithShape& ac){ac.setNumContacts(p_idx,0);}, particleAccessor);
+ contactStorage->forEachContact(useOpenMP, kernel::SelectAll(), contactAccessor,
+ [](size_t c, data::ContactAccessor &ca, ParticleAccessorWithShape &pa) {
+ auto idx1 = ca.getId1(c);
+ auto idx2 = ca.getId2(c);
+ pa.getNumContactsRef(idx1)++;
+ pa.getNumContactsRef(idx2)++;
+ }, contactAccessor, particleAccessor);
+
+ reductionKernel.operator()<NumContactNotification>(*particleStorage);
+ timing.stop("Contact eval");
+
+
+ timing.start("Shaking");
+ if(isShakingActive)
+ {
+ real_t shaking_common_term = 2_r * math::pi / shaking_period;
+ real_t shakingAcceleration = shaking_amplitude * std::sin((currentTime - timeBeginShaking) * shaking_common_term) * shaking_common_term * shaking_common_term;
+ if(solver == "DEM")
+ {
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ [shakingAcceleration](const size_t idx, auto &ac){addForceAtomic(idx, ac, Vec3(shakingAcceleration,0_r,0_r) / ac.getInvMass(idx));}, particleAccessor);
+ } else
+ {
+ hcsits_initParticles.setGlobalAcceleration(Vector3<real_t>(shakingAcceleration,0_r,-reducedGravitationalAcceleration));
+ }
+ }
+ timing.stop("Shaking");
+
+ if(solver == "HCSITS")
+ {
+ timing.start("HCSITS");
+
+ timing.start("Init contacts");
+ contactStorage->forEachContact(useOpenMP, kernel::SelectAll(), contactAccessor,
+ hcsits_initContacts, contactAccessor, particleAccessor);
+ timing.stop("Init contacts");
+ timing.start("Init particles");
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectAll(), particleAccessor,
+ hcsits_initParticles, particleAccessor, dt);
+ timing.stop("Init particles");
+
+ timing.start("Velocity update");
+ VelocityUpdateNotification::Parameters::relaxationParam = real_t(1.0); // must be set to 1.0 such that dv and dw caused by external forces and torques are not falsely altered
+ reductionKernel.operator()<VelocityCorrectionNotification>(*particleStorage);
+ broadcastKernel.operator()<VelocityUpdateNotification>(*particleStorage);
+ timing.stop("Velocity update");
+
+ VelocityUpdateNotification::Parameters::relaxationParam = hcsits_relaxationParameter;
+ for(uint_t i = uint_t(0); i < hcsits_numberOfIterations; i++){
+ timing.start("Relaxation step");
+ contactStorage->forEachContact(useOpenMP, kernel::SelectAll(), contactAccessor,
+ hcsits_relaxationStep, contactAccessor, particleAccessor, dt);
+ timing.stop("Relaxation step");
+ timing.start("Velocity update");
+ reductionKernel.operator()<VelocityCorrectionNotification>(*particleStorage);
+ broadcastKernel.operator()<VelocityUpdateNotification>(*particleStorage);
+ timing.stop("Velocity update");
+ }
+ timing.start("Integration");
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectAll(), particleAccessor,
+ hcsits_integration, particleAccessor, dt);
+ timing.stop("Integration");
+ timing.stop("HCSITS");
+ }
+ else if(solver == "DEM")
+ {
+ 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){
+ auto idx1 = ca.getId1(c);
+ auto idx2 = ca.getId2(c);
+ auto meff = real_t(1) / (pa.getInvMass(idx1) + pa.getInvMass(idx2));
+
+ /*
+ // if given stiffness
+ dem_collision.setStiffnessN(0,0,dem_stiffnessNormal);
+ dem_collision.setStiffnessT(0,0,dem_kappa*dem_stiffnessNormal);
+
+ // Wachs 2019, given stiffness and cor, we can compute damping (but formula in Wachs is probably wrong...)
+ auto log_en = std::log(coefficientOfRestitution);
+ auto dampingN = - 2_r * std::sqrt(dem_stiffnessNormal * meff) * log_en / (log_en*log_en+math::pi*math::pi);
+ dem_collision.setDampingN(0,0,dampingN);
+ dem_collision.setDampingT(0,0,std::sqrt(dem_kappa) * dampingN);
+ */
+
+ dem_collision.setStiffnessAndDamping(0,0,coefficientOfRestitution, dem_collisionTime, dem_kappa, meff);
+
+ dem_collision(idx1, idx2, pa, ca.getPosition(c), ca.getNormal(c), ca.getDistance(c), dt);
+ },
+ contactAccessor, particleAccessor);
+ timing.stop("Collision");
+
+
+ timing.start("Apply gravity");
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ [reducedGravitationalAcceleration](const size_t idx, auto &ac){addForceAtomic(idx, ac, Vec3(0_r,0_r,-reducedGravitationalAcceleration) / ac.getInvMass(idx));}, particleAccessor);
+ timing.stop("Apply gravity");
+
+ timing.start("Reduce");
+ reduceAndSwapContactHistory(*particleStorage);
+ reductionKernel.operator()<ForceTorqueNotification>(*particleStorage);
+ timing.stop("Reduce");
+
+ timing.start("Integration");
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ dem_integration, particleAccessor);
+ timing.stop("Integration");
+
+ timing.stop("DEM");
+
+ }
+
+ if(limitVelocity > 0_r)
+ {
+ timing.start("Velocity limiting");
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ [limitVelocity](const size_t idx, auto &ac){
+ auto velMagnitude = ac.getLinearVelocity(idx).length();
+ if(velMagnitude > limitVelocity) ac.getLinearVelocityRef(idx) *= (limitVelocity / velMagnitude );
+ }, particleAccessor);
+ timing.stop("Velocity limiting");
+ }
+
+
+ timing.start("Sync");
+ syncCall();
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ associateToBlock, particleAccessor);
+ timing.stop("Sync");
+
+ timing.start("Evaluate particles");
+ auto particleInfo = evaluateParticleInfo(particleAccessor);
+
+ if(particleInfo.particleVolume * particleDensity < totalParticleMass)
+ {
+
+ timing.start("Generation");
+ // check if generation
+ if(particleInfo.maximumHeight < generationHeightRatioStart * simulationDomain.zSize() - generationSpacing || currentTime - timeLastCreation > maximumTimeBetweenCreation)
+ {
+ particleCreator.createParticles( std::max(minGenerationHeight, generationHeightRatioStart * simulationDomain.zMax()),
+ std::min(maxGenerationHeight, generationHeightRatioEnd * simulationDomain.zMax()),
+ generationSpacing, diameterGenerator, shapeGenerator, initialVelocity, maximumAllowedInteractionRadius);
+
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ associateToBlock, particleAccessor);
+
+ if(useNextNeighborSync){ syncCall(); }
+ else { for(uint_t i = 0; i < uint_c(std::ceil(maxParticleDiameter/smallestBlockSize)); ++i) syncCall(); }
+
+ timeLastCreation = currentTime;
+
+ // write current particle distribution info
+ particleHistogram.clear();
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ particleHistogram, particleAccessor);
+ particleHistogram.evaluate();
+ WALBERLA_LOG_INFO_ON_ROOT(particleHistogram);
+ }
+ timing.stop("Generation");
+ } else if(shaking)
+ {
+ timing.start("Shaking");
+ // apply shaking
+ if(timeEndShaking < 0_r){
+ if(!isShakingActive)
+ {
+ isShakingActive = true;
+ timeBeginShaking = currentTime;
+ timeEndShaking = currentTime + shaking_duration;
+ WALBERLA_LOG_INFO_ON_ROOT("Beginning of shaking at time " << currentTime << " s for " << shaking_duration << " s.");
+ } else {
+ //timeEndShaking = real_c(std::ceil((currentTime + shaking_duration) / shaking_period)) * shaking_period; // make sure to shake only full periods
+ timeEndShaking = currentTime + shaking_duration; // since its unclear if full periods are really necessary and actually "improve" results, we skip this here
+ WALBERLA_LOG_INFO_ON_ROOT("Continue of shaking at time " << currentTime << " s until time " << timeEndShaking << " s.");
+ }
+ }
+
+ if(currentTime > timeEndShaking)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Ending of shaking at time " << currentTime << " s.");
+ shaking = false;
+ isShakingActive = false;
+ }
+ timing.stop("Shaking");
+
+ } else
+ {
+ timing.start("Damping");
+
+ if(timeBeginDamping < 0_r){
+ timeBeginDamping = currentTime;
+ WALBERLA_LOG_INFO_ON_ROOT("Beginning of damping at time " << timeBeginDamping << " s with damping factor " << velocityDampingFactor << " until convergence");
+ }
+
+ // apply damping
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectAll(), particleAccessor,
+ [velocityDampingFactor](size_t idx, ParticleAccessorWithShape &ac){
+ ac.getLinearVelocityRef(idx) *= velocityDampingFactor;
+ ac.getAngularVelocityRef(idx) *= velocityDampingFactor;},
+ particleAccessor);
+
+ // check if termination
+ if(currentTime - timeBeginDamping > minimalTerminalRunTime)
+ {
+ if(currentTime - timeLastTerminationCheck > terminationCheckingSpacing)
+ {
+ if(particleInfo.maximumVelocity < terminalVelocity)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Reached terminal max velocity - terminating.");
+ terminateSimulation = true;
+ }
+
+ real_t relDiffAvgHeight = std::abs(particleInfo.heightOfMass - oldAvgParticleHeight) / oldAvgParticleHeight;
+ real_t relDiffMaxHeight = std::abs(particleInfo.maximumHeight - oldMaxParticleHeight) / oldMaxParticleHeight;
+ if(relDiffMaxHeight < 10_r * terminalRelativeHeightChange && relDiffAvgHeight < terminalRelativeHeightChange)
+ {
+ // check of max height has to be included to avoid early termination if only little mass is created per generation step
+ WALBERLA_LOG_INFO_ON_ROOT("Reached converged maximum and mass-averaged height - terminating.");
+ terminateSimulation = true;
+ }
+
+ oldAvgParticleHeight = particleInfo.heightOfMass;
+ oldMaxParticleHeight = particleInfo.maximumHeight;
+ timeLastTerminationCheck = currentTime;
+ }
+ }
+ timing.stop("Damping");
+ }
+ timing.stop("Evaluate particles");
+
+ if(( infoSpacing > 0 && timestep % infoSpacing == 0) || (loggingSpacing > 0 && timestep % loggingSpacing == 0))
+ {
+ timing.start("Evaluate infos");
+ auto contactInfo = evaluateContactInfo(contactAccessor);
+
+ porosityEvaluator.clear();
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ porosityEvaluator, particleAccessor);
+ porosityEvaluator.evaluate();
+ real_t estimatedPorosity = porosityEvaluator.estimateTotalPorosity();
+
+
+ if(loggingSpacing > 0 && timestep % loggingSpacing == 0)
+ {
+ loggingWriter(currentTime, particleInfo, contactInfo, estimatedPorosity);
+ }
+
+ if(infoSpacing > 0 && timestep % infoSpacing == 0) {
+ WALBERLA_LOG_INFO_ON_ROOT("t = " << timestep << " = " << currentTime << " s");
+ WALBERLA_LOG_INFO_ON_ROOT(particleInfo << " => " << particleInfo.particleVolume * particleDensity << " kg" << ", current porosity = " << estimatedPorosity);
+ real_t ensembleAverageDiameter = diameterFromSphereVolume(particleInfo.particleVolume / real_c(particleInfo.numParticles));
+ WALBERLA_LOG_INFO_ON_ROOT(contactInfo << " => " << contactInfo.maximumPenetrationDepth / ensembleAverageDiameter * real_t(100) << "% of avg diameter " << ensembleAverageDiameter);
+ }
+
+ timing.stop("Evaluate infos");
+ }
+
+ ++timestep;
+ }
+
+ if(timing.isTimerRunning("Evaluate particles")) timing.stop("Evaluate particles");
+
+ timing.stop("Simulation");
+
+ particleHistogram.clear();
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ particleHistogram, particleAccessor);
+ particleHistogram.evaluate();
+ WALBERLA_LOG_INFO_ON_ROOT(particleHistogram);
+
+ porosityEvaluator.clear();
+ particleStorage->forEachParticle(useOpenMP, kernel::SelectLocal(), particleAccessor,
+ porosityEvaluator, particleAccessor);
+ porosityEvaluator.evaluate();
+ real_t estimatedFinalPorosity = porosityEvaluator.estimateTotalPorosity();
+ WALBERLA_LOG_INFO_ON_ROOT("Estimated total porosity based on layers = " << estimatedFinalPorosity);
+
+ std::string porosityFileName = porosityProfileFolder + "/" + uniqueFileIdentifier + "_layers.txt";
+ WALBERLA_LOG_INFO_ON_ROOT("Writing porosity profile file to " << porosityFileName);
+ porosityEvaluator.printToFile(porosityFileName);
+
+ ContactInfoPerHorizontalLayerEvaluator contactEvaluator(evaluationLayerHeight, simulationDomain);
+ contactStorage->forEachContact(useOpenMP, kernel::SelectAll(), particleAccessor,
+ contactEvaluator, contactAccessor);
+ contactEvaluator.evaluate();
+ std::string contactInfoFileName = porosityProfileFolder + "/" + uniqueFileIdentifier + "_contact_layers.txt";
+ WALBERLA_LOG_INFO_ON_ROOT("Writing contact info profile file to " << contactInfoFileName);
+ contactEvaluator.printToFile(contactInfoFileName);
+
+ auto reducedTT = timing.getReduced();
+ WALBERLA_LOG_INFO_ON_ROOT(reducedTT);
+
+ bool logToProcessLocalFiles = false;
+ std::string particleInfoFileName = porosityProfileFolder + "/" + uniqueFileIdentifier + "_particle_info";
+ if(logToProcessLocalFiles)
+ {
+ particleInfoFileName += "_" + std::to_string(walberla::mpi::MPIManager::instance()->rank()) + ".txt";
+ WALBERLA_LOG_INFO_ON_ROOT("Writing particle info file to process local files like " << particleInfoFileName);
+
+ } else {
+ particleInfoFileName += ".txt";
+ WALBERLA_LOG_INFO_ON_ROOT("Writing particle info file to " << particleInfoFileName);
+ }
+ auto particleInfoString = assembleParticleInformation(*particleStorage, particleSizeEvaluator, 12);
+ writeParticleInformationToFile(particleInfoFileName, particleInfoString, logToProcessLocalFiles);
+
+ // write to sqlite data base
+ auto particleInfo = evaluateParticleInfo(particleAccessor);
+ auto contactInfo = evaluateContactInfo(contactAccessor);
+
+ WALBERLA_ROOT_SECTION() {
+ std::map<std::string, walberla::int64_t> sql_integerProperties;
+ std::map<std::string, double> sql_realProperties;
+ std::map<std::string, std::string> sql_stringProperties;
+ addConfigToDatabase(*cfg, sql_integerProperties, sql_realProperties, sql_stringProperties);
+
+ // store particle info
+ sql_integerProperties["numParticles"] = int64_c(particleInfo.numParticles);
+ sql_realProperties["maxParticlePosition"] = double(particleInfo.maximumHeight);
+ sql_realProperties["particleVolume"] = double(particleInfo.particleVolume);
+
+ // store contact info
+ sql_integerProperties["numContacts"] = int64_c(contactInfo.numContacts);
+ sql_realProperties["maxPenetrationDepth"] = double(contactInfo.maximumPenetrationDepth);
+ sql_realProperties["avgPenetrationDepth"] = double(contactInfo.averagePenetrationDepth);
+
+ // other info
+ sql_realProperties["simulationTime"] = double(reducedTT["Simulation"].total());
+ sql_integerProperties["numProcesses"] = int64_c(walberla::mpi::MPIManager::instance()->numProcesses());
+ sql_integerProperties["timesteps"] = int64_c(timestep);
+ sql_stringProperties["file_identifier"] = uniqueFileIdentifier;
+ sql_realProperties["generationSpacing"] = double(generationSpacing);
+ std::string histogramData = "";
+ for (auto h : particleHistogram.getMassFractionHistogram()) histogramData += std::to_string(h) + " ";
+ sql_stringProperties["evaluation_histogramData"] = histogramData;
+ std::string numberHistogramData = "";
+ for (auto h : particleHistogram.getNumberHistogram()) numberHistogramData += std::to_string(h) + " ";
+ sql_stringProperties["evaluation_numberHistogramData"] = numberHistogramData;
+ sql_integerProperties["singleShape"] = (shapeGenerator->generatesSingleShape()) ? 1 : 0;
+ sql_realProperties["maxAllowedInteractionRadius"] = double(maximumAllowedInteractionRadius);
+
+ for(uint_t i = 0; i < particleHistogram.getNumberOfShapeEvaluators(); ++i)
+ {
+ std::string shapeBins = "";
+ for (auto b : particleHistogram.getShapeBins(i)) shapeBins += std::to_string(b) + " ";
+ sql_stringProperties["evaluation_"+std::get<0>(particleHistogram.getShapeEvaluator(i))+"_bins"] = shapeBins;
+
+ std::string shapeHistogram = "";
+ for (auto h : particleHistogram.getShapeHistogram(i)) shapeHistogram += std::to_string(h) + " ";
+ sql_stringProperties["evaluation_"+std::get<0>(particleHistogram.getShapeEvaluator(i))+"_histogramData"] = shapeHistogram;
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("Storing run and timing data in sql database file " << sqlDBFileName);
+ auto sql_runID = sqlite::storeRunInSqliteDB(sqlDBFileName, sql_integerProperties, sql_stringProperties, sql_realProperties);
+ sqlite::storeTimingTreeInSqliteDB(sqlDBFileName, sql_runID, reducedTT, "Timing");
+ }
+
+ if(!vtkFinalFolder.empty())
+ {
+
+ WALBERLA_LOG_INFO_ON_ROOT("Writing final VTK file to folder " << vtkFinalFolder);
+ if(particleShape.find("Mesh") != std::string::npos)
+ {
+ mesa_pd::MeshParticleVTKOutput< mesh::PolyMesh > finalMeshParticleVTK(particleStorage, uniqueFileIdentifier, uint_t(1), vtkFinalFolder);
+ finalMeshParticleVTK.addFaceOutput< data::SelectParticleUid >("UID");
+ finalMeshParticleVTK.addVertexOutput< data::SelectParticleInteractionRadius >("InteractionRadius");
+ finalMeshParticleVTK.addFaceOutput< data::SelectParticleLinearVelocity >("LinearVelocity");
+ finalMeshParticleVTK.addVertexOutput< data::SelectParticlePosition >("Position");
+ finalMeshParticleVTK.addVertexOutput< data::SelectParticleNumContacts>("numContacts");
+ finalMeshParticleVTK.addVertexDataSource(surfaceVelDataSource);
+ finalMeshParticleVTK.setParticleSelector(vtkParticleSelector);
+ finalMeshParticleVTK(particleAccessor);
+ }
+ else {
+ auto finalParticleVtkWriter = walberla::vtk::createVTKOutput_PointData(particleVtkOutput, uniqueFileIdentifier, uint_t(1), vtkFinalFolder, "final");
+ finalParticleVtkWriter->write();
+ }
+
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ofstream file;
+ std::string configFileCopyName = vtkFinalFolder + "/" + uniqueFileIdentifier + ".cfg";
+ WALBERLA_LOG_INFO_ON_ROOT("Storing config file as " << configFileCopyName);
+ file.open( configFileCopyName.c_str() );
+ file << *cfg;
+ file.close();
+ }
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("Simulation terminated successfully");
+
+ return EXIT_SUCCESS;
+}
+} // namespace msa_pd
+} // namespace walberla
+
+int main( int argc, char* argv[] ) {
+ return walberla::mesa_pd::main( argc, argv );
+}
\ No newline at end of file
diff --git a/apps/showcases/ParticlePacking/ShapeGeneration.h b/apps/showcases/ParticlePacking/ShapeGeneration.h
new file mode 100644
index 000000000..7dc8dd5ab
--- /dev/null
+++ b/apps/showcases/ParticlePacking/ShapeGeneration.h
@@ -0,0 +1,566 @@
+//======================================================================================================================
+//
+// 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 ShapeGeneration.h
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/config/Config.h"
+#include "core/Filesystem.h"
+
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/shape/Sphere.h"
+#include "mesa_pd/data/shape/Ellipsoid.h"
+#include "mesa_pd/data/shape/ConvexPolyhedron.h"
+
+#include "mesh_common/TriangleMeshes.h"
+#include "mesh_common/MeshIO.h"
+#include "mesh_common/MeshOperations.h"
+#include "mesh_common/PolyMeshes.h"
+
+#include "Utility.h"
+
+#include <random>
+
+namespace walberla {
+namespace mesa_pd {
+
+Vec3 getSemiAxesFromMesh(mesh::TriangleMesh& mesh)
+{
+ auto volume = mesh::computeVolume(mesh);
+ mesh::translate(mesh, mesh::toWalberla(-mesh::computeCentroid(mesh)));
+ auto inertiaTensor = mesh::computeInertiaTensor(mesh);
+ auto mass = volume * 1_r; // mesh has unit mass with density 1
+ return semiAxesFromInertiaTensor(inertiaTensor, mass);
+}
+
+std::vector<Vec3> extractSemiAxesFromMeshFiles(const std::vector<std::string> & meshFiles)
+{
+ std::vector<Vec3> semiAxesVector;
+
+ for(const auto& meshFileName : meshFiles)
+ {
+ mesh::TriangleMesh mesh;
+ mesh::readAndBroadcast<mesh::TriangleMesh>(meshFileName, mesh);
+ auto volume = mesh::computeVolume(mesh);
+
+ auto semiAxes = getSemiAxesFromMesh(mesh);
+
+ WALBERLA_LOG_INFO_ON_ROOT("Read mesh file: " << meshFileName << " (" << mesh.n_vertices() << " vertices, "
+ << mesh.n_faces()
+ << " faces, volume = " << volume << ", semi axes = " << semiAxes << ")");
+
+ semiAxesVector.push_back(semiAxes);
+ }
+
+ return semiAxesVector;
+}
+
+
+std::vector<std::string> getMeshFilesFromPath(const std::string & meshPath)
+{
+ std::vector<std::string> meshNames;
+
+ if(filesystem::is_directory(filesystem::path(meshPath)))
+ {
+ // assuming this folder contains the mesh files
+ for(const auto& entry : filesystem::directory_iterator(meshPath)) {
+ std::string meshFileName = entry.path();
+ if(meshFileName.find(".obj") == std::string::npos)
+ {
+ // open mesh seemingly can only read .obj files reliably, so skip all others
+ WALBERLA_LOG_WARNING_ON_ROOT("Ignoring mesh file " << meshFileName << ", since not of .obj file type.");
+ continue;
+ }
+ meshNames.push_back(meshFileName);
+ }
+ } else {
+ // assuming path is a single (mesh) file
+ meshNames.push_back(meshPath);
+ }
+
+ return meshNames;
+}
+
+enum ScaleMode { sphereEquivalent, sieveLike };
+
+ScaleMode str_to_scaleMode(const std::string & str)
+{
+ if(str == "sphereEquivalent") return ScaleMode::sphereEquivalent;
+ else if(str == "sieveLike") return ScaleMode::sieveLike;
+ else WALBERLA_ABORT("Unknown shape scale mode " << str);
+}
+
+using FormParameters = Vector3<real_t>; // contains L, I, S
+
+class NormalizedFormGenerator
+{
+public:
+ virtual FormParameters get() = 0;
+ virtual FormParameters getNormalFormParameters() const = 0;
+ virtual real_t getMaxDiameterScalingFactor() const = 0; // used to estimate maximum particle extent
+ virtual real_t getNormalVolume() const = 0; // volume of a typical particle with diameter 1
+ virtual bool generatesSingleShape() const = 0; // generate a single shape or multiple shapes?
+ virtual ~NormalizedFormGenerator() = default;
+};
+
+class ConstFormGenerator : public NormalizedFormGenerator
+{
+public:
+ explicit ConstFormGenerator(){}
+ FormParameters get() override { return FormParameters(1_r, 1_r, 1_r); }
+ FormParameters getNormalFormParameters() const override { return FormParameters(1_r, 1_r, 1_r); }
+ real_t getMaxDiameterScalingFactor() const override { return 1_r; }
+ real_t getNormalVolume() const override { return 1_r; }
+ bool generatesSingleShape() const override { return true; }
+};
+
+class SampleFormGenerator : public NormalizedFormGenerator
+{
+public:
+ SampleFormGenerator(const std::vector<Vec3> & semiAxesVector,
+ ScaleMode scaleMode) :
+ gen_(static_cast<unsigned long>(walberla::mpi::MPIManager::instance()->rank()))
+ {
+
+ maxDiameterScalingFactor_ = 1_r;
+ normalVolume_ = 0_r;
+ for(const auto& semiAxes : semiAxesVector)
+ {
+ FormParameters normalizedFormParameters;
+ switch(scaleMode)
+ {
+ case ScaleMode::sphereEquivalent:
+ {
+ real_t ellipsoidScalingFactor = std::cbrt( 1_r / ( 8_r * semiAxes[0] * semiAxes[1] * semiAxes[2] ) ); // = V_sphere / V_ellipsoid
+ normalizedFormParameters = 2_r * semiAxes * ellipsoidScalingFactor;
+ break;
+ }
+ case ScaleMode::sieveLike:
+ {
+ auto scalingFactor = sizeFromSemiAxes(semiAxes);
+ normalizedFormParameters = 2_r * semiAxes / scalingFactor;
+ break;
+ }
+ default:
+ WALBERLA_ABORT("Unknown shape scale mode");
+ }
+
+ auto volume = 1_r / 6_r * math::pi * normalizedFormParameters[0] * normalizedFormParameters[1] * normalizedFormParameters[2];
+
+ maxDiameterScalingFactor_ = std::max(maxDiameterScalingFactor_, normalizedFormParameters.max());
+ normalVolume_ += volume;
+
+ normalizedFormParametersVector_.push_back(normalizedFormParameters);
+ }
+
+ WALBERLA_CHECK(!normalizedFormParametersVector_.empty(), "No shape samples provided");
+
+ auto numFormParameters = normalizedFormParametersVector_.size();
+ normalVolume_ /= real_c(numFormParameters); // use average normal volume here as estimation
+
+ normalFormParameters_ = std::accumulate(normalizedFormParametersVector_.begin(),
+ normalizedFormParametersVector_.end(),
+ FormParameters(0_r)) / real_t(numFormParameters); // = average of form parameters
+
+ dist_ = std::uniform_int_distribution<uint_t>(0, numFormParameters-1);
+ }
+
+ FormParameters get() override
+ {
+ return normalizedFormParametersVector_[dist_(gen_)];
+ }
+
+ FormParameters getNormalFormParameters() const override { return normalFormParameters_; }
+ real_t getMaxDiameterScalingFactor() const override { return maxDiameterScalingFactor_; }
+ real_t getNormalVolume() const override { return normalVolume_; }
+ bool generatesSingleShape() const override { return normalizedFormParametersVector_.size() == 1; }
+
+private:
+ std::vector<FormParameters> normalizedFormParametersVector_;
+ std::mt19937 gen_;
+ std::uniform_int_distribution<uint_t> dist_;
+
+ real_t maxDiameterScalingFactor_;
+ real_t normalVolume_;
+ FormParameters normalFormParameters_;
+
+};
+
+// assumes individual normal distribution of form factors elongation (= I/L) and flatness (= S/I)
+// NOTE: since S <= I <= L, the combination elongation & equancy would not work trivially as it follows: equancy <= elongation,
+// these two parameters are not be completely independent
+// limits both values to interval (eps,1]
+// -> best way would be to use a truncated Normal distribution https://en.wikipedia.org/wiki/Truncated_normal_distribution
+// however, formula is more complicated and no library for sampling is available
+class DistributionFormGenerator : public NormalizedFormGenerator
+{
+public:
+ DistributionFormGenerator( real_t elongationMean, real_t elongationStdDev,
+ real_t flatnessMean, real_t flatnessStdDev,
+ ScaleMode scaleMode) :
+ scaleMode_(scaleMode), gen_(static_cast<unsigned long>(walberla::mpi::MPIManager::instance()->rank()))
+ {
+ elongationDist_ = std::normal_distribution<real_t>(elongationMean, elongationStdDev);
+ flatnessDist_ = std::normal_distribution<real_t>(flatnessMean, flatnessStdDev);
+
+ normalFormParameters_ = getNormalizedFormParameters(elongationMean, flatnessMean);
+ normalVolume_ = 1_r / 6_r * math::pi * normalFormParameters_[0] * normalFormParameters_[1] * normalFormParameters_[2];
+
+ WALBERLA_LOG_INFO_ON_ROOT("Shape generation from distribution with mean size-normalized form parameters (L,I,S) = " << normalFormParameters_);
+
+ real_t extremeElongation = std::max(eps_, elongationMean - 5_r * elongationStdDev); // covers almost all values
+ real_t extremeFlatness = std::max(eps_, flatnessMean - 5_r * flatnessStdDev);
+ auto extremeFormParams = getNormalizedFormParameters(extremeElongation, extremeFlatness);
+ maxDiameterScalingFactor_ = extremeFormParams.max();
+
+ }
+
+ FormParameters get() override
+ {
+ // variant: re-roll until within bounds
+ real_t elongation = elongationDist_(gen_);
+ while(elongation < eps_ || elongation > 1_r) elongation = elongationDist_(gen_);
+ real_t flatness = flatnessDist_(gen_);
+ while(flatness < eps_ || flatness > 1_r) flatness = flatnessDist_(gen_);
+
+ // variant: cap values to min/max
+ //real_t elongation = std::min(std::max(eps_, elongationFromDist), 1_r);
+ //real_t flatness = std::min(std::max(eps_, flatnessFromDist), 1_r);
+
+ return getNormalizedFormParameters(elongation, flatness);
+ }
+
+ FormParameters getNormalFormParameters() const override { return normalFormParameters_; }
+ real_t getNormalVolume() const override { return normalVolume_; }
+ real_t getMaxDiameterScalingFactor() const override { return maxDiameterScalingFactor_; }
+ bool generatesSingleShape() const override { return false; }
+
+private:
+
+ FormParameters getNormalizedFormParameters(real_t elongation, real_t flatness)
+ {
+ real_t I = 1_r;
+ switch(scaleMode_)
+ {
+ case ScaleMode::sieveLike:
+ I = std::sqrt(2_r) / std::sqrt(1_r + flatness * flatness); break;
+ case ScaleMode::sphereEquivalent:
+ I = std::cbrt(elongation / flatness); break;
+ }
+ return FormParameters( I / elongation, I, I * flatness);
+ }
+
+ ScaleMode scaleMode_;
+ std::mt19937 gen_;
+ std::normal_distribution<real_t> elongationDist_;
+ std::normal_distribution<real_t> flatnessDist_;
+ real_t normalVolume_;
+ real_t maxDiameterScalingFactor_;
+ FormParameters normalFormParameters_;
+
+ const real_t eps_ = 0.1_r;
+};
+
+
+class ShapeGenerator
+{
+public:
+ virtual void setShape(real_t diameter, real_t maximumAllowedInteractionRadius, data::ParticleStorage::baseShape_type& shape, real_t& interactionRadius) = 0;
+ virtual real_t getMaxDiameterScalingFactor() = 0; // used to estimate maximum particle extent
+ virtual real_t getNormalVolume() = 0; // volume of a typical particle with diameter 1
+ virtual FormParameters getNormalFormParameters() = 0; // form (L,I,S) of a typical particle with diameter 1
+ virtual bool generatesSingleShape() = 0; // generate a single shape or multiple shapes?
+ virtual ~ShapeGenerator() = default;
+};
+
+class SphereGenerator : public ShapeGenerator
+{
+public:
+ explicit SphereGenerator()
+ {
+ maxDiameterScalingFactor_ = 1_r;
+ normalVolume_ = math::pi / 6_r;
+ WALBERLA_LOG_INFO_ON_ROOT("Will create spheres");
+ }
+
+ void setShape(real_t diameter, real_t maximumAllowedInteractionRadius, data::ParticleStorage::baseShape_type& shape, real_t& interactionRadius) override
+ {
+ real_t radius = diameter * 0.5_r;
+ if(radius > maximumAllowedInteractionRadius) radius = maximumAllowedInteractionRadius;
+ shape = std::make_shared<data::Sphere>(radius);
+ interactionRadius = radius;
+ }
+
+ real_t getMaxDiameterScalingFactor() override { return maxDiameterScalingFactor_; }
+ real_t getNormalVolume() override { return normalVolume_; }
+ FormParameters getNormalFormParameters() override {return FormParameters(1_r);}
+ bool generatesSingleShape() override { return true; }
+
+private:
+ real_t maxDiameterScalingFactor_;
+ real_t normalVolume_;
+};
+
+
+class EllipsoidGenerator : public ShapeGenerator
+{
+public:
+ EllipsoidGenerator(shared_ptr<NormalizedFormGenerator> normalizedFormGenerator) :
+ normalizedFormGenerator_(normalizedFormGenerator){}
+
+ void setShape(real_t diameter, real_t maximumAllowedInteractionRadius, data::ParticleStorage::baseShape_type& shape, real_t& interactionRadius) override
+ {
+ Vector3<real_t> semiAxes = 0.5_r * normalizedFormGenerator_->get() * diameter;
+ sortVector(semiAxes); // we want to have particles that are longest in z-direction
+
+ real_t maxAxis = semiAxes.max();
+ if(maxAxis > maximumAllowedInteractionRadius)
+ {
+ semiAxes *= (maximumAllowedInteractionRadius / maxAxis);
+ }
+
+ shape = std::make_shared<data::Ellipsoid>(semiAxes);
+ interactionRadius = semiAxes.max();
+ }
+
+ real_t getMaxDiameterScalingFactor() override { return normalizedFormGenerator_->getMaxDiameterScalingFactor(); }
+ real_t getNormalVolume() override { return normalizedFormGenerator_->getNormalVolume(); }
+ FormParameters getNormalFormParameters() override {return normalizedFormGenerator_->getNormalFormParameters();}
+ bool generatesSingleShape() override { return normalizedFormGenerator_->generatesSingleShape(); }
+
+private:
+ shared_ptr<NormalizedFormGenerator> normalizedFormGenerator_;
+};
+
+
+class MeshesGenerator : public ShapeGenerator
+{
+public:
+ MeshesGenerator(const std::vector<std::string> & meshFiles, ScaleMode scaleMode, shared_ptr<NormalizedFormGenerator> normalizedFormGenerator) :
+ normalizedFormGenerator_(normalizedFormGenerator), gen_(static_cast<unsigned long>(walberla::mpi::MPIManager::instance()->rank()))
+ {
+ WALBERLA_CHECK(!meshFiles.empty());
+
+ maxDiameterScalingFactor_ = 1_r;
+ normalVolume_ = 0_r;
+ for(const auto& meshFileName : meshFiles)
+ {
+ mesh::TriangleMesh mesh;
+ mesh::readAndBroadcast<mesh::TriangleMesh>(meshFileName, mesh);
+ mesh::translate(mesh, mesh::toWalberla(-mesh::computeCentroid(mesh)));
+
+ if(normalizedFormGenerator_->generatesSingleShape())
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Read mesh file: " << meshFileName << " (" << mesh.n_vertices() << " vertices, "
+ << mesh.n_faces()
+ << " faces, volume = " << mesh::computeVolume(mesh) << ")");
+
+ // all shape scaling is handled as given by the meshes and thus by this class internally
+ switch( scaleMode )
+ {
+ case ScaleMode::sphereEquivalent:
+ {
+ real_t meshScalingFactor = std::cbrt(math::pi / (6_r * mesh::computeVolume(mesh)));
+ mesh::scale(mesh, Vec3(meshScalingFactor));
+ break;
+ }
+ case ScaleMode::sieveLike:
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Using sieve-like scaling! Assuming that the mesh is scaled such that a scaling with the mesh size results in the correct size!");
+ break;
+ }
+ default: WALBERLA_ABORT("Unknown shape scale mode");
+ }
+
+ maxDiameterScalingFactor_ = std::max( maxDiameterScalingFactor_, 2_r * mesh::computeBoundingSphereRadius(mesh, mesh::computeCentroid(mesh)) );
+ normalVolume_ += mesh::computeVolume(mesh);
+
+ } else
+ {
+ // meshes are scaled via form factors provided from outside
+ // -> scale mesh such that the semi axes are 0.5, i.e. removing all shape features from it
+
+ auto semiAxes = getSemiAxesFromMesh(mesh);
+ auto meshScalingFactor = Vec3(1_r / (2_r * semiAxes[0]), 1_r / (2_r * semiAxes[1]), 1_r / (2_r * semiAxes[2]));
+ mesh::scale(mesh, meshScalingFactor);
+
+ maxDiameterScalingFactor_ = std::max( maxDiameterScalingFactor_,
+ 2_r * mesh::computeBoundingSphereRadius(mesh, mesh::computeCentroid(mesh)) * normalizedFormGenerator_->getMaxDiameterScalingFactor() );
+ auto meshCopy = mesh;
+ auto normalFormParameters = normalizedFormGenerator_->getNormalFormParameters();
+ sortVector(normalFormParameters);
+ mesh::scale(meshCopy, normalFormParameters);
+ normalVolume_ += mesh::computeVolume(meshCopy);
+ }
+
+ particleMeshes_.push_back(mesh);
+
+ }
+
+ WALBERLA_CHECK(!particleMeshes_.empty());
+
+ normalVolume_ /= real_c(particleMeshes_.size()); // use average normal mesh volume here as estimation
+
+ dist_ = std::uniform_int_distribution<uint_t>(0, particleMeshes_.size()-1);
+
+ WALBERLA_LOG_INFO_ON_ROOT("Read and stored in total " << particleMeshes_.size() << " meshes, and will randomly pick one for each generated particle.")
+ }
+
+ void setShape(real_t diameter, real_t maximumAllowedInteractionRadius, data::ParticleStorage::baseShape_type& shape, real_t& interactionRadius) override
+ {
+ auto meshCopy = particleMeshes_[dist_(gen_)];
+ auto normalizedFormParameters = normalizedFormGenerator_->get();
+ sortVector(normalizedFormParameters); // we want to have particles that are longest in z-direction
+
+ mesh::scale(meshCopy, normalizedFormParameters * diameter);
+ auto convexPolyPtr = std::make_shared<data::ConvexPolyhedron>(meshCopy);
+ convexPolyPtr->updateMeshQuantities();
+ real_t boundingRadius = convexPolyPtr->getBoundingSphereRadius();
+
+ if(boundingRadius > maximumAllowedInteractionRadius)
+ {
+ // scale to match limiting radius
+ mesh::scale(meshCopy, Vec3(maximumAllowedInteractionRadius / boundingRadius) );
+ convexPolyPtr = std::make_shared<data::ConvexPolyhedron>(meshCopy);
+ convexPolyPtr->updateMeshQuantities();
+ }
+ shape = convexPolyPtr;
+ interactionRadius = convexPolyPtr->getBoundingSphereRadius();
+ }
+
+ real_t getMaxDiameterScalingFactor() override { return maxDiameterScalingFactor_; }
+ real_t getNormalVolume() override { return normalVolume_; }
+ FormParameters getNormalFormParameters() override {return normalizedFormGenerator_->getNormalFormParameters();}
+ bool generatesSingleShape() override { return particleMeshes_.size() == 1 && normalizedFormGenerator_->generatesSingleShape(); }
+
+private:
+ shared_ptr<NormalizedFormGenerator> normalizedFormGenerator_;
+
+ std::vector<mesh::TriangleMesh> particleMeshes_;
+ std::mt19937 gen_;
+ std::uniform_int_distribution<uint_t> dist_;
+
+ real_t maxDiameterScalingFactor_;
+ real_t normalVolume_;
+};
+
+class UnscaledMeshesPerFractionGenerator : public ShapeGenerator
+{
+public:
+ UnscaledMeshesPerFractionGenerator(const Config::BlockHandle & shapeConfig, const std::vector<real_t> & massFractions) :
+ gen_(static_cast<unsigned long>(walberla::mpi::MPIManager::instance()->rank()))
+ {
+
+ auto meshesConfig = shapeConfig.getBlock("UnscaledMeshesPerFraction");
+ std::string meshesTopFolder = meshesConfig.getParameter<std::string>("folder");
+
+ std::vector<real_t> avgVolumesPerFraction(massFractions.size(), 0_r);
+
+ maxDiameterScalingFactor_ = 1_r;
+ generatesSingleShape_ = true;
+
+ auto meshesTopFolderPath = filesystem::path(meshesTopFolder);
+ for(uint_t fractionIdx = 0; fractionIdx < massFractions.size(); ++fractionIdx)
+ {
+ auto meshesFolder = meshesTopFolderPath / std::to_string(fractionIdx);
+
+ if(!filesystem::exists(meshesFolder))
+ {
+ WALBERLA_ABORT("Path " << meshesFolder.string() << " expected but does not exist.");
+ }
+
+ std::vector<mesh::TriangleMesh> meshesVector;
+
+ for (const auto &entry : filesystem::directory_iterator(meshesFolder)) {
+ std::string meshFileName = entry.path();
+ if (meshFileName.find(".obj") == std::string::npos) {
+ // open mesh seemingly can only read .obj files reliably, so skip all others
+ continue;
+ }
+ mesh::TriangleMesh mesh;
+ mesh::readAndBroadcast<mesh::TriangleMesh>(meshFileName, mesh);
+
+ WALBERLA_LOG_INFO_ON_ROOT("Read mesh file: " << meshFileName << " (" << mesh.n_vertices() << " vertices, "
+ << mesh.n_faces()
+ << " faces, volume = " << mesh::computeVolume(mesh) << ")");
+ mesh::translate(mesh, mesh::toWalberla(-mesh::computeCentroid(mesh)));
+ meshesVector.push_back(mesh);
+
+ avgVolumesPerFraction[fractionIdx] += mesh::computeVolume(mesh);
+
+ /*
+ auto meshAABB = mesh::computeAABB(mesh);
+ auto aabbHeight = meshAABB.zSize();
+ auto aabbHorizontalDiagonal = std::sqrt(meshAABB.xSize()*meshAABB.xSize() + meshAABB.ySize()*meshAABB.ySize());
+ maxDiameterScalingFactor_ = std::max(maxDiameterScalingFactor_, aabbHeight / aabbHorizontalDiagonal);
+ */
+ maxDiameterScalingFactor_ = std::max( maxDiameterScalingFactor_, 2_r * mesh::computeBoundingSphereRadius(mesh, mesh::computeCentroid(mesh)) );
+ }
+
+ WALBERLA_CHECK(!particleMeshPerFractionVector_.empty());
+
+ particleMeshPerFractionVector_.push_back(meshesVector);
+
+ if(particleMeshPerFractionVector_.size() != 1) generatesSingleShape_ = false;
+
+ avgVolumesPerFraction[fractionIdx] /= real_c(meshesVector.size()); // use average normal mesh volume here as estimation
+
+ distsPerFraction_.push_back(std::uniform_int_distribution<uint_t>(0, meshesVector.size()-1));
+
+ }
+
+ auto particleNumbers = transferMassFractionsToParticleNumbersFromAvgVolumes(massFractions, avgVolumesPerFraction);
+ auto totalParticles = std::accumulate(particleNumbers.begin(), particleNumbers.end(), real_t(0));
+ std::string outString = "Particle probabilities per size fraction: ";
+ for(const auto & p : particleNumbers) outString += std::to_string(p/totalParticles) + " | ";
+ WALBERLA_LOG_INFO_ON_ROOT(outString);
+ distForSizeFraction_ = std::discrete_distribution<uint_t>(particleNumbers.begin(), particleNumbers.end());
+
+ WALBERLA_LOG_INFO_ON_ROOT("Read and stored in total " << particleMeshPerFractionVector_.size() << " size fractions and their meshes, and will randomly pick one for each to-be generated size fraction.")
+ }
+
+ void setShape(real_t /*diameter*/, real_t maximumAllowedInteractionRadius, data::ParticleStorage::baseShape_type& shape, real_t& interactionRadius) override
+ {
+ auto sizeFractionIdx = distForSizeFraction_(gen_);
+ auto meshCopy = particleMeshPerFractionVector_[sizeFractionIdx][distsPerFraction_[sizeFractionIdx](gen_)];
+ auto convexPolyPtr = std::make_shared<data::ConvexPolyhedron>(meshCopy);
+ convexPolyPtr->updateMeshQuantities();
+ shape = convexPolyPtr;
+ interactionRadius = convexPolyPtr->getBoundingSphereRadius();
+ WALBERLA_CHECK_GREATER_EQUAL(maximumAllowedInteractionRadius, interactionRadius, "Particle shape larger than allowed radius")
+ }
+
+ real_t getMaxDiameterScalingFactor() override { return maxDiameterScalingFactor_; }
+ real_t getNormalVolume() override { return 1_r; } //dummy value
+ FormParameters getNormalFormParameters() override {return FormParameters(1_r);} // dummy value
+ bool generatesSingleShape() override { return generatesSingleShape_; }
+
+private:
+ std::vector<std::vector<mesh::TriangleMesh>> particleMeshPerFractionVector_;
+ std::mt19937 gen_;
+ std::discrete_distribution<uint_t> distForSizeFraction_;
+ std::vector<std::uniform_int_distribution<uint_t>> distsPerFraction_;
+
+ real_t maxDiameterScalingFactor_;
+ bool generatesSingleShape_;
+};
+
+
+} // namespace msa_pd
+} // namespace walberla
diff --git a/apps/showcases/ParticlePacking/Utility.h b/apps/showcases/ParticlePacking/Utility.h
new file mode 100644
index 000000000..5e0f0b4c6
--- /dev/null
+++ b/apps/showcases/ParticlePacking/Utility.h
@@ -0,0 +1,221 @@
+//======================================================================================================================
+//
+// 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 Utility.h
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "mesa_pd/data/ParticleStorage.h"
+
+#include <iterator>
+#include <algorithm>
+#include <functional>
+
+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)
+{
+ std::istringstream iss(inputStr);
+ return std::vector<T>{std::istream_iterator<T>(iss),std::istream_iterator<T>()};
+}
+
+real_t radiusFromSphereVolume(real_t volume)
+{
+ return std::cbrt(real_t(3) / ( real_t(4) * math::pi) * volume);
+}
+
+real_t diameterFromSphereVolume(real_t volume)
+{
+ return real_t(2) * radiusFromSphereVolume(volume);
+}
+
+uint_t getIndexOfSecondSemiAxis(Vector3<real_t> semiAxes)
+{
+ std::set<uint_t> indices = {0,1,2};
+ indices.erase(semiAxes.indexOfMin());
+ indices.erase(semiAxes.indexOfMax());
+ return *indices.begin();
+}
+
+void sortVector(Vector3<real_t> & v)
+{
+ std::sort(v.data(),v.data()+3);
+}
+
+real_t sizeFromSemiAxes(Vec3 semiAxes)
+{
+ sortVector(semiAxes);
+ return 2_r * std::sqrt((semiAxes[0] * semiAxes[0] + semiAxes[1] * semiAxes[1]) / 2_r); // factor of 2 because those are SEMI axis
+}
+
+// assuming a mass-equivalent ellipsoid
+Vec3 semiAxesFromInertiaTensor(const Matrix3<real_t> & inertiaTensor, real_t mass)
+{
+ Vec3 semiAxes;
+ semiAxes[0] = std::sqrt( (- inertiaTensor(0,0) + inertiaTensor(1,1) + inertiaTensor(2,2)) / (2_r * mass / 5_r) );
+ semiAxes[1] = std::sqrt( (- inertiaTensor(1,1) + inertiaTensor(2,2) + inertiaTensor(0,0)) / (2_r * mass / 5_r) );
+ semiAxes[2] = std::sqrt( (- inertiaTensor(2,2) + inertiaTensor(0,0) + inertiaTensor(1,1)) / (2_r * mass / 5_r) );
+ return semiAxes;
+}
+
+
+std::vector<real_t> getMeanDiametersFromSieveSizes(std::vector<real_t> sieveSizes)
+{
+ //since grain sizes are logarithmically distributed, it is practice to use the geometric mean
+ std::vector<real_t> meanDiameters(sieveSizes.size()-1,0_r);
+ for(uint_t i = 0; i < meanDiameters.size(); ++i)
+ {
+ meanDiameters[i] = std::sqrt(sieveSizes[i] * sieveSizes[i+1]);
+ }
+ return meanDiameters;
+}
+
+// if totalMass and density are given actual numbers, the resulting particle numbers are a good estimate for the actual numbers
+// else, the resulting numbers are directly proportional to the actual ones, which is sufficient to define the distributions
+std::vector<real_t> transferMassFractionsToParticleNumbersFromAvgVolumes(std::vector<real_t> massFractions, std::vector<real_t> avgVolumePerSizeFraction, real_t totalMass = real_t(1), real_t density = real_t(1) )
+{
+ WALBERLA_CHECK_EQUAL(avgVolumePerSizeFraction.size(), massFractions.size(), "Number of entries in volume and mass-fraction array has to be the same!");
+ std::vector<real_t> particleNumbers(massFractions.size(), real_t(0));
+ for(uint_t n = 0; n < massFractions.size(); ++n )
+ {
+ if(avgVolumePerSizeFraction[n] > 0_r) particleNumbers[n] = totalMass * massFractions[n] / (density * avgVolumePerSizeFraction[n]);
+ }
+ return particleNumbers;
+}
+
+// note: normalVolume is the volume of a typical particle with diameter = 1. For sphere: PI / 6
+std::vector<real_t> transferMassFractionsToParticleNumbers(std::vector<real_t> massFractions, std::vector<real_t> diameters, real_t normalVolume = math::pi / real_t(6), real_t totalMass = real_t(1), real_t density = real_t(1) )
+{
+ WALBERLA_CHECK_EQUAL(diameters.size(), massFractions.size(), "Number of entries in diameter and mass-fraction array has to be the same!");
+ std::vector<real_t> avgVolumePerSizeFraction(massFractions.size(), real_t(0));
+ for(uint_t n = 0; n < massFractions.size(); ++n )
+ {
+ avgVolumePerSizeFraction[n] = normalVolume * diameters[n] * diameters[n] * diameters[n];
+ }
+ return transferMassFractionsToParticleNumbersFromAvgVolumes(massFractions, avgVolumePerSizeFraction, totalMass, density);
+}
+
+real_t computePercentileFromSieveDistribution(std::vector<real_t> diameters, std::vector<real_t> massFractions, real_t percentile)
+{
+ WALBERLA_CHECK(percentile > 0_r && percentile < 100_r, "Percentile is a value between 0 and 100");
+ if(diameters[0] > diameters[1])
+ {
+ // reverse order to have it ascending
+ std::reverse(diameters.begin(), diameters.end());
+ std::reverse(massFractions.begin(), massFractions.end());
+ }
+
+ std::vector<real_t> cdf(massFractions.size(),0_r);
+ std::partial_sum(massFractions.begin(), massFractions.end(), cdf.begin());
+
+ for(uint_t i = 0; i < cdf.size()-1; ++i)
+ {
+ if(cdf[i] <= percentile/100_r && percentile/100_r <= cdf[i+1] )
+ {
+ real_t f_small = cdf[i];
+ real_t f_large = cdf[i+1];
+ if(f_small <= 0.000001_r && f_large>= 0.999999_r)
+ {
+ // special case of uniform distribution -> percentile is this value
+ return diameters[i+1];
+ }
+ real_t phi_small = - std::log2(diameters[i]);
+ real_t phi_large = - std::log2(diameters[i+1]);
+ // logarithmic interpolation of diameter value
+ real_t phi_percentile = phi_small + (phi_large - phi_small) / (f_large - f_small) * (percentile / 100_r - f_small);
+ return std::pow(2_r,-phi_percentile);
+ }
+ }
+ return diameters[0];
+}
+
+
+auto createPlane( std::shared_ptr<data::ParticleStorage> particleStorage,
+ const Vec3& pos,
+ const Vec3& normal)
+{
+ auto p = particleStorage->create(true);
+ p->setPosition( pos );
+ p->setBaseShape( std::make_shared<data::HalfSpace>( normal ) );
+ p->getBaseShapeRef()->updateMassAndInertia(real_t(1));
+ p->setOwner( walberla::mpi::MPIManager::instance()->rank() );
+ p->setType( 0 );
+ p->setInteractionRadius(std::numeric_limits<real_t>::infinity());
+ data::particle_flags::set(p->getFlagsRef(), data::particle_flags::GLOBAL);
+ data::particle_flags::set(p->getFlagsRef(), data::particle_flags::INFINITE);
+ data::particle_flags::set(p->getFlagsRef(), data::particle_flags::FIXED);
+ data::particle_flags::set(p->getFlagsRef(), data::particle_flags::NON_COMMUNICATING);
+ return p;
+}
+
+auto createCylindricalBoundary( std::shared_ptr<data::ParticleStorage> particleStorage,
+ const Vec3& pos,
+ const Vec3& axis, real_t radius)
+{
+ auto p = particleStorage->create(true);
+ p->setPosition( pos );
+ p->setBaseShape( std::make_shared<data::CylindricalBoundary>( radius, axis ) );
+ p->getBaseShapeRef()->updateMassAndInertia(real_t(1));
+ p->setOwner( walberla::mpi::MPIManager::instance()->rank() );
+ p->setType( 0 );
+ p->setInteractionRadius(std::numeric_limits<real_t>::infinity());
+ data::particle_flags::set(p->getFlagsRef(), data::particle_flags::GLOBAL);
+ data::particle_flags::set(p->getFlagsRef(), data::particle_flags::INFINITE);
+ data::particle_flags::set(p->getFlagsRef(), data::particle_flags::FIXED);
+ data::particle_flags::set(p->getFlagsRef(), data::particle_flags::NON_COMMUNICATING);
+ return p;
+}
+
+
+void writeParticleInformationToFile(const std::string& filename, const std::string& particleInfoStr, bool logToProcessLocalFiles)
+{
+
+ if(logToProcessLocalFiles)
+ {
+ std::ofstream file;
+ file.open( filename.c_str());
+ file << particleInfoStr;
+ file.close();
+ }else{
+ walberla::mpi::writeMPITextFile( filename, particleInfoStr );
+ }
+
+}
+
+
+
+} // namespace mesa_pd
+} // namespace walberla
diff --git a/apps/showcases/ParticlePacking/example_meshes/sediment_scan_0.obj b/apps/showcases/ParticlePacking/example_meshes/sediment_scan_0.obj
new file mode 100644
index 000000000..3e112080b
--- /dev/null
+++ b/apps/showcases/ParticlePacking/example_meshes/sediment_scan_0.obj
@@ -0,0 +1,605 @@
+# https://github.com/mikedh/trimesh
+v -0.02417045 0.06592044 -1.07176822
+v 0.13701385 0.54986295 -0.32406061
+v -0.21109703 0.10191524 -1.04785955
+v -0.23173924 -0.11902797 -0.31848047
+v -0.20100322 0.46167435 -0.52743332
+v -0.16944822 -0.52708383 0.72240008
+v -0.18419047 0.50212121 -0.03590216
+v -0.14945623 0.65596573 0.15793829
+v 0.08926398 0.20941177 0.88442000
+v 0.12003405 0.40262039 -0.64842151
+v -0.10598839 0.54413552 -0.40983931
+v -0.06298423 -0.61676329 0.63621720
+v 0.20271825 -0.30186448 -0.83229467
+v -0.24879295 -0.27231811 -0.86163188
+v 0.04391830 -0.51534681 0.74546334
+v -0.18286594 -0.17005737 1.03575934
+v -0.15656668 0.64640682 0.22858417
+v -0.18751951 -0.10724870 0.89847957
+v -0.20935047 -0.26993524 0.29303015
+v -0.14595073 -0.60504029 0.54468882
+v -0.20520294 -0.49832057 -0.18549337
+v -0.18073674 -0.27971435 0.96059046
+v 0.25499010 -0.31437120 -0.30803213
+v 0.22693859 -0.04129109 -0.85514519
+v 0.18224998 -0.38697172 -0.49015189
+v -0.13577139 -0.31664914 1.13875314
+v -0.17297744 -0.23922282 -1.11252084
+v 0.30411918 -0.01793443 0.10738888
+v 0.00522421 -0.05353117 1.09420037
+v -0.15517983 0.47787656 -0.62953709
+v -0.07065463 0.72914195 0.46672341
+v -0.15264593 0.13743894 1.11535444
+v 0.19839563 -0.43213791 -0.05483209
+v 0.21830344 0.29043442 -0.57605616
+v -0.02551993 0.72007305 0.25401898
+v -0.11729200 0.47431677 -0.64609953
+v -0.23591093 0.01443627 -0.45699604
+v -0.17296194 0.10061511 -1.06053375
+v -0.05332434 0.63163879 0.65227659
+v 0.18244612 0.47495787 -0.42886181
+v 0.01932635 0.38064080 0.86412652
+v -0.10669568 -0.54256761 -0.09785274
+v -0.14029039 0.70525294 0.47714366
+v -0.11048504 0.07172404 1.15524953
+v -0.13254648 0.50390584 0.85700419
+v -0.25643367 -0.07026606 -1.03665507
+v -0.00709990 -0.19368527 -1.16722244
+v -0.08497270 0.23838325 1.06455407
+v -0.14419796 -0.09753769 -1.14501427
+v -0.20162666 -0.16344693 -1.10403124
+v -0.18373502 -0.51868949 0.44525045
+v -0.04093191 0.43992541 -0.71114451
+v -0.13015878 -0.60997436 0.76568277
+v 0.30602227 0.02872141 0.33874770
+v 0.15045553 -0.51994738 0.39987427
+v -0.05533047 0.49706873 0.84861910
+v -0.11896789 0.65548723 0.02219109
+v -0.23234262 -0.41114605 -0.59733837
+v -0.00665313 -0.59700535 0.66833287
+v -0.00270839 -0.49132225 0.87739673
+v 0.27954805 -0.28229195 0.14280300
+v 0.06894202 0.55778562 0.44074600
+v -0.24950699 -0.23864670 -1.02268686
+v 0.02060672 -0.44591946 -0.49414559
+v 0.05029771 0.01112007 -1.07017569
+v -0.12650600 -0.20506517 1.17733673
+v -0.11433789 -0.27124076 -1.13051493
+v -0.16746322 -0.30497314 1.09117918
+v -0.21913297 0.33229493 -0.69143837
+v -0.10547765 0.69116228 0.63098075
+v 0.09790875 0.60154710 0.14037367
+v 0.11248413 -0.41007869 -0.51836913
+v 0.26192693 -0.21905407 -0.57860031
+v -0.06747283 -0.00858800 -1.14834922
+v 0.17666887 0.53758446 -0.19071287
+v 0.29127225 0.23623759 -0.11927450
+v -0.07475661 0.39767138 0.96159716
+v -0.17970170 0.11463361 0.92730013
+v -0.19322328 0.47308074 -0.59480505
+v 0.20881343 -0.23226355 0.70832137
+v 0.23079811 -0.22324964 -0.87991447
+v 0.16303689 -0.15794833 0.84703971
+v 0.15223775 0.42899705 -0.57884201
+v 0.27494377 0.26224738 -0.38610899
+v 0.12091936 -0.32226195 -0.88274361
+v 0.17083455 -0.44229374 0.60920733
+v 0.30363055 -0.14166875 0.39045921
+v -0.16598957 -0.22112602 1.14655259
+v 0.16150960 0.06499542 0.83793685
+v 0.21322138 -0.35277946 0.59921303
+v 0.29166533 -0.21245815 0.04678303
+v -0.13326804 -0.59635217 0.36074351
+v 0.15811502 0.23615215 -0.75820653
+v -0.08677741 0.25878461 -0.90656025
+v 0.17578636 -0.34570283 0.74041289
+v -0.20488176 0.02342118 -1.08926242
+v -0.19583805 -0.49081559 -0.61378289
+v 0.26619287 -0.27944981 0.39809376
+v 0.28997738 0.13449034 -0.18954945
+v 0.23485040 -0.26667359 -0.77750266
+v 0.14171863 -0.18881531 -1.07053739
+v 0.26020932 0.29717898 -0.42937116
+v 0.22072077 -0.05363683 0.72720844
+v -0.23238384 -0.28456313 -1.04791497
+v 0.01072801 0.68660295 0.09772671
+v 0.08650958 -0.53816191 0.34234171
+v -0.22382958 -0.43557221 -0.29750290
+v 0.22676622 0.10212319 0.64750130
+v -0.22356174 -0.38535323 -0.84416995
+v 0.20681067 -0.13184801 -0.95386008
+v 0.13829534 -0.08050160 -1.06378358
+v 0.29205244 -0.10588419 0.48693610
+v -0.15390263 -0.56164323 0.79692871
+v 0.01485938 -0.28267833 -1.09128489
+v 0.21041291 -0.42180924 0.44285785
+v -0.16697704 -0.48984718 -0.59114652
+v 0.17105320 0.19876175 0.74776826
+v 0.14881639 -0.25359354 -1.03958479
+v 0.25881792 0.35661871 -0.23267050
+v 0.04581930 -0.57807580 0.58409278
+v -0.01498864 0.31560793 0.95852587
+v 0.10764352 0.37538363 0.67767586
+v -0.02773468 -0.34049219 1.03697828
+v -0.14434263 0.65781187 0.64086780
+v 0.17675836 -0.34998568 -0.72226417
+v 0.21883715 -0.35063740 -0.48583098
+v 0.00865317 -0.04599846 -1.14435263
+v -0.18982153 -0.46227809 -0.70143077
+v -0.14184498 0.00675651 1.17630888
+v -0.20488419 0.37660339 -0.73097947
+v -0.03346433 -0.15354642 1.13243989
+v -0.15525674 0.47947441 0.79749020
+v 0.06221444 0.46125796 -0.62863527
+v -0.17124050 0.07423043 1.10298066
+v 0.06083507 0.57686678 -0.29321431
+v -0.21256637 -0.36744179 0.12204940
+v -0.16884064 -0.57287792 0.49333432
+v 0.24923133 -0.09917100 -0.81962288
+v 0.21916392 -0.41684200 0.09972426
+v -0.04521411 -0.06914949 1.14305424
+v -0.21186959 -0.48442159 -0.53062260
+v -0.25292086 0.01098005 -1.02206174
+v -0.11590075 0.72958551 0.35863016
+v -0.11873659 -0.52099838 0.89883145
+v 0.00695034 0.15864935 1.02882417
+v 0.03682127 0.35646576 -0.75632509
+v -0.11271422 -0.62501167 0.67733490
+v -0.20909119 -0.35467410 -0.98251866
+v 0.25695214 -0.14944222 -0.75090442
+v 0.24004279 -0.37888717 0.20388126
+v -0.08009215 -0.42143618 1.00480173
+v 0.00294116 -0.59436910 0.54371550
+vn 0.30662637 0.63673860 -0.70749150
+vn 0.48594324 0.84136855 -0.23655471
+vn -0.59545613 0.51373335 -0.61766499
+vn -0.99938379 -0.00168868 0.03505984
+vn -0.96855069 0.24824074 -0.01691430
+vn -0.97817234 -0.19755919 0.06441457
+vn -0.99321896 0.11271673 0.02847856
+vn -0.84494664 0.52912626 -0.07804217
+vn 0.80033438 0.31554035 0.50980307
+vn 0.57639549 0.59566992 -0.55941539
+vn -0.07682513 0.96072533 -0.26665473
+vn 0.20973173 -0.97774500 -0.00522644
+vn 0.67227299 -0.69268296 -0.26122662
+vn -0.99560347 -0.09366237 -0.00104855
+vn 0.67855005 -0.65399043 0.33446427
+vn -0.99572524 -0.01777373 0.09063849
+vn -0.96059038 0.27789867 -0.00620095
+vn -0.99902227 0.01873866 0.04004216
+vn -0.99906927 -0.02064105 0.03787532
+vn -0.61765474 -0.78634387 -0.01288197
+vn -0.93887475 -0.34407112 0.01136965
+vn -0.99452819 -0.08314658 0.06324807
+vn 0.92603056 -0.37354840 -0.05412016
+vn 0.90032444 0.26069624 -0.34850162
+vn 0.54581582 -0.82804713 -0.12815242
+vn -0.16808277 -0.54279336 0.82287517
+vn -0.50800733 -0.25080002 -0.82403150
+vn 0.99956824 -0.01139242 -0.02708424
+vn 0.75864421 0.07384395 0.64730676
+vn -0.27055095 0.86163075 -0.42941195
+vn 0.29515647 0.93699045 0.18689988
+vn -0.60759517 0.40507381 0.68318615
+vn 0.70298663 -0.70621818 -0.08405764
+vn 0.88388634 0.33972801 -0.32144954
+vn 0.34973198 0.93554882 -0.04935537
+vn -0.06978817 0.87938846 -0.47096236
+vn -0.99858172 0.04113282 0.03380285
+vn -0.13190192 0.61138430 -0.78026349
+vn 0.69800042 0.61940028 0.35935874
+vn 0.75301498 0.60616773 -0.25598656
+vn 0.75137678 0.46536248 0.46783619
+vn 0.19396447 -0.97341059 -0.12185899
+vn -0.75112213 0.65376548 0.09168557
+vn 0.22904987 0.33481720 0.91402057
+vn -0.50943310 0.62496899 0.59151642
+vn -0.94139644 0.02383883 -0.33645868
+vn 0.20787949 -0.21454446 -0.95433579
+vn 0.23423045 0.49055294 0.83934135
+vn -0.37690130 0.04770761 -0.92502400
+vn -0.67379006 -0.05954556 -0.73651971
+vn -0.98508399 -0.16778811 0.03816641
+vn 0.08891564 0.80043315 -0.59279068
+vn -0.32176027 -0.88846570 0.32725987
+vn 0.97974740 0.17449606 0.09821488
+vn 0.64314771 -0.76573868 -0.00230190
+vn 0.47978869 0.66906263 0.56758966
+vn -0.32778768 0.91972692 -0.21600375
+vn -0.97869939 -0.20444580 -0.01869291
+vn 0.47121638 -0.85492713 0.21692101
+vn 0.55488371 -0.65854308 0.50835526
+vn 0.95849094 -0.28511640 0.00193642
+vn 0.84348058 0.50288737 0.18879831
+vn -0.96673886 -0.12838960 -0.22120597
+vn 0.31267877 -0.93792048 -0.15012384
+vn 0.52280523 0.54904458 -0.65209258
+vn -0.02287402 -0.18733950 0.98202886
+vn -0.15659212 -0.57346133 -0.80412749
+vn -0.91801582 -0.27670426 0.28404525
+vn -0.98725038 0.15529714 -0.03492105
+vn 0.02124534 0.91044119 0.41309257
+vn 0.82283350 0.55893127 0.10266867
+vn 0.34069545 -0.92797505 -0.15095998
+vn 0.98593575 -0.15571115 -0.06070204
+vn 0.01006380 0.39894944 -0.91691770
+vn 0.79270597 0.60944703 -0.01384082
+vn 0.97771674 0.20623826 0.03918873
+vn 0.19563185 0.59992861 0.77576662
+vn -0.99679529 0.06110720 0.05162413
+vn -0.73146682 0.64445182 -0.22279620
+vn 0.93668366 -0.15101473 0.31594031
+vn 0.92633858 -0.24375656 -0.28719258
+vn 0.86539519 -0.11511735 0.48768756
+vn 0.68987646 0.60781237 -0.39323606
+vn 0.98262437 0.15803819 -0.09733080
+vn 0.32379119 -0.89868461 -0.29584666
+vn 0.80146488 -0.54682201 0.24215643
+vn 0.98816009 -0.12587588 0.08772062
+vn -0.84020002 -0.12262913 0.52822914
+vn 0.87146083 0.14220982 0.46939577
+vn 0.92836465 -0.30620057 0.21066630
+vn 0.99043035 -0.13397675 -0.03313533
+vn -0.14572037 -0.98496948 -0.09273991
+vn 0.74761288 0.44581314 -0.49226581
+vn 0.02763226 0.72509548 -0.68809374
+vn 0.84949950 -0.34929949 0.39539911
+vn -0.57824862 0.26085485 -0.77303511
+vn -0.36990713 -0.90252956 -0.22047476
+vn 0.95093278 -0.29340045 0.09819892
+vn 0.99851673 0.00912691 -0.05367526
+vn 0.91993444 -0.37091489 -0.12705424
+vn 0.71311266 -0.10112813 -0.69371712
+vn 0.93460204 0.30105752 -0.18942911
+vn 0.93727180 0.01940533 0.34805893
+vn -0.74636140 -0.42855057 -0.50920435
+vn 0.30900704 0.94044980 -0.14166447
+vn 0.34783142 -0.92945708 -0.12297495
+vn -0.98742855 -0.15656319 0.02174452
+vn 0.94215116 0.23785132 0.23617353
+vn -0.86769026 -0.47449536 -0.14821527
+vn 0.90757137 0.07094594 -0.41386095
+vn 0.71266689 0.26179335 -0.65082266
+vn 0.97628762 -0.04015657 0.21272032
+vn -0.85423580 -0.44693413 0.26557687
+vn 0.26380333 -0.77475362 -0.57459955
+vn 0.88854963 -0.44931689 0.09270320
+vn 0.22287323 -0.95727949 -0.18423760
+vn 0.88753949 0.32528995 0.32628225
+vn 0.64842674 -0.54637243 -0.53011313
+vn 0.93712581 0.34897985 0.00287776
+vn 0.54842389 -0.82876883 0.11123607
+vn 0.68124287 0.40462563 0.61007070
+vn 0.85131840 0.45279272 0.26502024
+vn 0.65647182 -0.42398144 0.62392667
+vn -0.82229911 0.49862678 0.27421800
+vn 0.51077804 -0.83798980 -0.19203879
+vn 0.80069964 -0.59082032 -0.09905273
+vn 0.38173479 0.31539478 -0.86879496
+vn -0.05886705 -0.94585153 -0.31921710
+vn -0.36037258 0.12736400 0.92407252
+vn -0.69977815 0.57113206 -0.42909056
+vn 0.65599839 -0.15869877 0.73788943
+vn -0.97341343 0.17384688 0.14914276
+vn 0.32632658 0.82482839 -0.46170237
+vn -0.95516748 0.10757692 0.27583017
+vn 0.12572482 0.95370719 -0.27319564
+vn -0.99726510 -0.06266898 0.03917802
+vn -0.90568191 -0.42360360 0.01732836
+vn 0.97323408 0.07173797 -0.21833254
+vn 0.84544938 -0.53321332 -0.02998159
+vn 0.54500673 0.07563401 0.83501327
+vn -0.82698559 -0.56088103 -0.03882392
+vn -0.94921847 0.15349068 -0.27463595
+vn -0.43859486 0.89690019 -0.05660925
+vn -0.18584224 -0.78958380 0.58482483
+vn 0.70439878 0.28253287 0.65115093
+vn 0.40885628 0.64035167 -0.65022017
+vn -0.15328284 -0.98534253 0.07486300
+vn -0.43805131 -0.81142069 -0.38692056
+vn 0.99296602 -0.06522323 -0.09881504
+vn 0.90316306 -0.42883584 0.01990758
+vn 0.34028586 -0.68156507 0.64782296
+vn 0.29522591 -0.95267460 -0.07247596
+f 114//114 47//47 118//118
+f 129//129 32//32 134//134
+f 129//129 44//44 32//32
+f 115//115 55//55 139//139
+f 123//123 131//131 26//26
+f 123//123 82//82 131//131
+f 123//123 95//95 82//82
+f 147//147 59//59 53//53
+f 49//49 74//74 47//47
+f 133//133 2//2 83//83
+f 75//75 2//2 105//105
+f 57//57 79//79 143//143
+f 43//43 31//31 143//143
+f 36//36 133//133 52//52
+f 67//67 47//47 114//114
+f 67//67 49//49 47//47
+f 85//85 114//114 118//118
+f 48//48 32//32 44//44
+f 48//48 77//77 32//32
+f 48//48 44//44 145//145
+f 132//132 134//134 32//32
+f 132//132 78//78 134//134
+f 29//29 89//89 145//145
+f 29//29 82//82 89//89
+f 29//29 131//131 82//82
+f 88//88 129//129 134//134
+f 66//66 26//26 131//131
+f 66//66 88//88 26//26
+f 66//66 129//129 88//88
+f 140//140 44//44 129//129
+f 140//140 66//66 131//131
+f 140//140 129//129 66//66
+f 140//140 131//131 29//29
+f 140//140 29//29 145//145
+f 140//140 145//145 44//44
+f 33//33 139//139 55//55
+f 33//33 55//55 25//25
+f 98//98 61//61 87//87
+f 86//86 55//55 115//115
+f 150//150 115//115 139//139
+f 150//150 98//98 115//115
+f 150//150 61//61 98//98
+f 91//91 87//87 61//61
+f 100//100 149//149 73//73
+f 112//112 87//87 54//54
+f 60//60 53//53 59//59
+f 60//60 95//95 123//123
+f 60//60 86//86 95//95
+f 152//152 42//42 64//64
+f 116//116 64//64 42//42
+f 116//116 42//42 97//97
+f 116//116 125//125 64//64
+f 116//116 85//85 125//125
+f 116//116 97//97 85//85
+f 72//72 64//64 125//125
+f 72//72 125//125 25//25
+f 92//92 97//97 42//42
+f 16//16 88//88 134//134
+f 16//16 134//134 78//78
+f 101//101 118//118 47//47
+f 101//101 47//47 111//111
+f 127//127 111//111 47//47
+f 127//127 47//47 74//74
+f 24//24 102//102 138//138
+f 10//10 133//133 83//83
+f 84//84 138//138 102//102
+f 84//84 149//149 138//138
+f 117//117 54//54 76//76
+f 40//40 83//83 2//2
+f 40//40 2//2 75//75
+f 35//35 143//143 31//31
+f 35//35 75//75 105//105
+f 35//35 57//57 143//143
+f 35//35 105//105 57//57
+f 70//70 43//43 124//124
+f 70//70 31//31 43//43
+f 130//130 3//3 142//142
+f 130//130 69//69 79//79
+f 130//130 142//142 69//69
+f 94//94 130//130 52//52
+f 94//94 3//3 130//130
+f 96//96 74//74 49//49
+f 96//96 49//49 50//50
+f 96//96 142//142 3//3
+f 8//8 43//43 143//143
+f 8//8 143//143 79//79
+f 37//37 69//69 142//142
+f 5//5 79//79 69//69
+f 5//5 8//8 79//79
+f 11//11 57//57 105//105
+f 135//135 105//105 2//2
+f 135//135 11//11 105//105
+f 135//135 36//36 11//11
+f 135//135 2//2 133//133
+f 135//135 133//133 36//36
+f 14//14 58//58 107//107
+f 141//141 107//107 58//58
+f 141//141 109//109 97//97
+f 141//141 58//58 109//109
+f 20//20 147//147 53//53
+f 20//20 53//53 137//137
+f 20//20 92//92 147//147
+f 20//20 137//137 141//141
+f 20//20 141//141 97//97
+f 20//20 97//97 92//92
+f 148//148 67//67 114//114
+f 148//148 97//97 109//109
+f 148//148 109//109 104//104
+f 148//148 104//104 67//67
+f 128//128 85//85 97//97
+f 128//128 114//114 85//85
+f 128//128 148//148 114//114
+f 128//128 97//97 148//148
+f 13//13 85//85 118//118
+f 13//13 125//125 85//85
+f 121//121 48//48 145//145
+f 121//121 77//77 48//48
+f 121//121 41//41 77//77
+f 45//45 132//132 32//32
+f 45//45 32//32 77//77
+f 45//45 124//124 132//132
+f 45//45 70//70 124//124
+f 7//7 78//78 132//132
+f 7//7 132//132 124//124
+f 7//7 37//37 78//78
+f 7//7 69//69 37//37
+f 7//7 5//5 69//69
+f 113//113 137//137 53//53
+f 90//90 98//98 87//87
+f 90//90 87//87 112//112
+f 90//90 95//95 86//86
+f 90//90 86//86 115//115
+f 90//90 115//115 98//98
+f 120//120 55//55 86//86
+f 120//120 152//152 55//55
+f 120//120 59//59 152//152
+f 23//23 150//150 139//139
+f 23//23 61//61 150//150
+f 23//23 100//100 73//73
+f 23//23 139//139 33//33
+f 23//23 91//91 61//61
+f 23//23 73//73 91//91
+f 28//28 54//54 87//87
+f 28//28 87//87 91//91
+f 28//28 76//76 54//54
+f 28//28 91//91 73//73
+f 28//28 73//73 149//149
+f 81//81 13//13 118//118
+f 81//81 100//100 13//13
+f 81//81 138//138 149//149
+f 81//81 149//149 100//100
+f 80//80 82//82 95//95
+f 80//80 90//90 112//112
+f 80//80 95//95 90//90
+f 103//103 89//89 82//82
+f 103//103 80//80 112//112
+f 103//103 82//82 80//80
+f 151//151 60//60 123//123
+f 151//151 123//123 26//26
+f 15//15 60//60 59//59
+f 15//15 86//86 60//60
+f 15//15 120//120 86//86
+f 15//15 59//59 120//120
+f 12//12 152//152 59//59
+f 12//12 59//59 147//147
+f 12//12 147//147 92//92
+f 12//12 92//92 42//42
+f 12//12 42//42 152//152
+f 106//106 152//152 64//64
+f 106//106 64//64 72//72
+f 106//106 55//55 152//152
+f 106//106 72//72 25//25
+f 106//106 25//25 55//55
+f 6//6 137//137 113//113
+f 46//46 50//50 104//104
+f 46//46 96//96 50//50
+f 46//46 142//142 96//96
+f 46//46 37//37 142//142
+f 27//27 67//67 104//104
+f 27//27 104//104 50//50
+f 27//27 50//50 49//49
+f 27//27 49//49 67//67
+f 1//1 94//94 52//52
+f 1//1 127//127 74//74
+f 1//1 74//74 94//94
+f 93//93 24//24 111//111
+f 93//93 10//10 83//83
+f 93//93 111//111 10//10
+f 110//110 24//24 138//138
+f 110//110 111//111 24//24
+f 110//110 138//138 81//81
+f 110//110 101//101 111//111
+f 110//110 81//81 118//118
+f 110//110 118//118 101//101
+f 65//65 10//10 111//111
+f 65//65 111//111 127//127
+f 65//65 127//127 1//1
+f 119//119 76//76 84//84
+f 119//119 84//84 102//102
+f 119//119 117//117 76//76
+f 119//119 75//75 117//117
+f 119//119 40//40 75//75
+f 119//119 102//102 40//40
+f 99//99 84//84 76//76
+f 99//99 149//149 84//84
+f 99//99 28//28 149//149
+f 99//99 76//76 28//28
+f 9//9 117//117 41//41
+f 9//9 121//121 145//145
+f 9//9 41//41 121//121
+f 9//9 145//145 89//89
+f 9//9 89//89 117//117
+f 108//108 112//112 54//54
+f 108//108 54//54 117//117
+f 108//108 103//103 112//112
+f 108//108 117//117 89//89
+f 108//108 89//89 103//103
+f 34//34 40//40 102//102
+f 34//34 83//83 40//40
+f 34//34 93//93 83//83
+f 34//34 102//102 24//24
+f 34//34 24//24 93//93
+f 56//56 77//77 41//41
+f 56//56 45//45 77//77
+f 56//56 70//70 45//45
+f 30//30 130//130 79//79
+f 30//30 11//11 36//36
+f 30//30 36//36 52//52
+f 30//30 52//52 130//130
+f 30//30 79//79 57//57
+f 30//30 57//57 11//11
+f 38//38 94//94 74//74
+f 38//38 3//3 94//94
+f 38//38 96//96 3//3
+f 38//38 74//74 96//96
+f 18//18 16//16 78//78
+f 18//18 78//78 37//37
+f 17//17 43//43 8//8
+f 17//17 8//8 5//5
+f 17//17 5//5 7//7
+f 17//17 7//7 124//124
+f 17//17 124//124 43//43
+f 63//63 109//109 58//58
+f 63//63 58//58 14//14
+f 63//63 104//104 109//109
+f 63//63 46//46 104//104
+f 63//63 14//14 46//46
+f 21//21 141//141 137//137
+f 21//21 107//107 141//141
+f 126//126 13//13 100//100
+f 126//126 100//100 23//23
+f 126//126 25//25 125//125
+f 126//126 125//125 13//13
+f 126//126 23//23 33//33
+f 126//126 33//33 25//25
+f 144//144 113//113 53//53
+f 144//144 26//26 113//113
+f 144//144 151//151 26//26
+f 144//144 53//53 60//60
+f 144//144 60//60 151//151
+f 68//68 113//113 26//26
+f 68//68 26//26 88//88
+f 68//68 6//6 113//113
+f 68//68 88//88 16//16
+f 22//22 68//68 16//16
+f 22//22 6//6 68//68
+f 51//51 137//137 6//6
+f 51//51 6//6 22//22
+f 51//51 21//21 137//137
+f 51//51 107//107 21//21
+f 146//146 1//1 52//52
+f 146//146 52//52 133//133
+f 146//146 133//133 10//10
+f 146//146 65//65 1//1
+f 146//146 10//10 65//65
+f 71//71 62//62 117//117
+f 71//71 117//117 75//75
+f 71//71 75//75 35//35
+f 71//71 35//35 31//31
+f 71//71 31//31 62//62
+f 122//122 41//41 117//117
+f 122//122 117//117 62//62
+f 39//39 62//62 31//31
+f 39//39 31//31 70//70
+f 39//39 122//122 62//62
+f 39//39 70//70 56//56
+f 39//39 56//56 41//41
+f 39//39 41//41 122//122
+f 4//4 18//18 37//37
+f 4//4 16//16 18//18
+f 4//4 19//19 16//16
+f 4//4 37//37 46//46
+f 4//4 46//46 14//14
+f 4//4 14//14 19//19
+f 136//136 19//19 14//14
+f 136//136 14//14 107//107
+f 136//136 22//22 16//16
+f 136//136 16//16 19//19
+f 136//136 51//51 22//22
+f 136//136 107//107 51//51
\ No newline at end of file
diff --git a/apps/showcases/ParticlePacking/example_meshes/sediment_scan_1.obj b/apps/showcases/ParticlePacking/example_meshes/sediment_scan_1.obj
new file mode 100644
index 000000000..353ae8a0f
--- /dev/null
+++ b/apps/showcases/ParticlePacking/example_meshes/sediment_scan_1.obj
@@ -0,0 +1,545 @@
+# https://github.com/mikedh/trimesh
+v -0.01682526 -0.47405221 0.48244358
+v 0.00436325 -0.47295352 0.54177185
+v -0.01259683 -0.47078728 0.68162064
+v -0.13549947 -0.47411451 0.52482807
+v -0.15245956 -0.47194827 0.66467686
+v -0.13127104 -0.47084958 0.72400513
+v 0.14958462 0.00827245 -0.92423798
+v -0.38444932 0.00799211 -0.73350780
+v -0.42259791 0.00905966 -0.65298728
+v 0.41230362 0.02798664 0.18605542
+v 0.28512619 0.03662044 0.76664282
+v 0.18803936 -0.02534154 0.86937979
+v -0.43493070 -0.04850755 -0.31293722
+v -0.12983468 -0.02659597 0.91601273
+v 0.44239422 -0.04260916 -0.29179501
+v 0.23050911 -0.05359609 -0.88507774
+v 0.18062000 0.47538058 0.66200906
+v 0.06194579 0.47531828 0.70439354
+v 0.09168388 0.45249490 -0.71103814
+v 0.21881495 0.45908706 -0.35506850
+v -0.12452239 0.46869497 0.36961614
+v 0.26542041 0.46454938 -0.03723489
+v 0.01115827 0.45136506 -0.74917417
+v -0.12029396 0.47195990 0.56879321
+v -0.09910544 0.47305860 0.62812148
+v 0.29083735 0.46891301 0.22127045
+v -0.04817884 0.45133391 -0.72798193
+v -0.08632743 0.45240146 -0.64746141
+v 0.29506578 0.47217794 0.42044751
+v 0.27810570 0.47434418 0.56029630
+v 0.12910112 -0.08837107 0.89159624
+v -0.35828101 -0.09841505 0.46360299
+v 0.20220738 0.41348089 0.72236154
+v 0.25736241 0.39502114 -0.43456481
+v 0.13023134 0.38842897 -0.79053444
+v 0.08353317 0.41341859 0.76474602
+v 0.30396787 0.40048346 -0.11673120
+v -0.03081989 0.38616929 -0.86680651
+v 0.32938481 0.40484708 0.14177414
+v -0.24702616 0.40236936 0.21384778
+v 0.33361324 0.40811202 0.34095120
+v -0.09015699 0.38613814 -0.84561426
+v -0.23856930 0.40889923 0.61220190
+v -0.19619228 0.41109661 0.73085845
+v 0.29969308 0.41244449 0.62064878
+v -0.16645418 0.38827323 -0.68457323
+v 0.18883709 -0.15133830 0.87142821
+v -0.40021280 -0.17883679 -0.59058638
+v 0.30328287 -0.15454094 0.62986666
+v -0.18837405 -0.15262389 0.93925339
+v 0.46438047 -0.16750723 -0.23041832
+v -0.30168961 0.32737879 -0.50131575
+v -0.31864969 0.32954503 -0.36146696
+v 0.34674376 0.33968246 0.00294956
+v -0.33138134 0.33497620 -0.02244111
+v 0.35520062 0.34621233 0.40130368
+v 0.33824053 0.34837857 0.54115247
+v -0.31869605 0.34477100 0.57509008
+v -0.25513052 0.34806708 0.75307490
+v -0.32351675 -0.24397027 -0.75060320
+v 0.22738455 -0.21540423 0.79193190
+v -0.36166534 -0.24290272 -0.67008269
+v 0.03250588 -0.24378338 -0.87775666
+v -0.36171171 -0.22767675 0.26647435
+v -0.33629477 -0.22331312 0.52497969
+v -0.18797518 -0.21562227 0.94027760
+v 0.26553314 -0.21647177 0.71141138
+v 0.46477933 -0.23050561 -0.22939411
+v 0.25289422 -0.24149254 -0.82267684
+v 0.22419362 0.28858282 0.78373823
+v 0.16485652 0.28855167 0.80493047
+v -0.34366776 0.26218302 -0.61894809
+v 0.08865205 0.26023482 -0.90714257
+v -0.39877644 0.26541681 -0.39857878
+v -0.39454801 0.26868174 -0.19940172
+v -0.36067421 0.27957523 0.45745774
+v -0.26737058 0.26004793 -0.77998912
+v 0.31326882 0.26579059 -0.65288569
+v 0.24970329 0.26249451 -0.83087051
+v -0.35703805 -0.30263617 -0.46988142
+v 0.41429795 -0.28700527 0.19117647
+v 0.32104069 -0.28270394 0.49206629
+v -0.31470739 -0.28521281 0.58533217
+v 0.44398969 -0.29460268 -0.28769817
+v -0.33485840 0.22094048 0.71698729
+v 0.38991851 0.21588309 0.12365452
+v -0.44075460 0.20022104 -0.51621112
+v -0.43652617 0.20348597 -0.31703406
+v -0.40688079 0.21111453 0.14064834
+v 0.32207819 0.22454804 0.68304968
+v -0.23314427 0.22316902 0.81445159
+v -0.38146385 0.21547816 0.39915368
+v 0.28392959 0.22561559 0.76357019
+v 0.20699377 -0.34249969 0.73465204
+v -0.28034200 -0.36776964 -0.62989824
+v 0.38505145 -0.35763221 -0.26548172
+v -0.31849059 -0.36670210 -0.54937772
+v 0.16884518 -0.34143214 0.81517256
+v 0.13501774 -0.36755160 -0.77824394
+v 0.37231980 -0.35220104 0.07354413
+v 0.35535972 -0.35003480 0.21339292
+v -0.32699381 -0.35800599 -0.01117481
+v -0.06850324 -0.34155673 0.89994153
+v 0.30025104 -0.34680102 0.43376222
+v -0.12351919 -0.36877489 -0.75280324
+v -0.31853695 -0.35147612 0.38717931
+v -0.18285630 -0.36880604 -0.73161100
+v -0.27193150 -0.34601380 0.70501292
+v -0.22955447 -0.34381642 0.82366947
+v 0.29606897 -0.36529192 -0.70197188
+v 0.38082302 -0.36089714 -0.46465878
+v -0.31327102 0.15904079 0.77733977
+v 0.14878689 0.13426921 -0.92628640
+v -0.38524705 0.13398887 -0.73555621
+v -0.46154424 0.13612397 -0.57451518
+v -0.45731581 0.13938890 -0.37533812
+v -0.15221979 0.16130048 0.85361183
+v 0.37763209 0.14308991 -0.47285246
+v 0.22931251 0.13539905 -0.88815037
+v -0.23756611 -0.42857064 -0.51021748
+v 0.23713073 -0.42832145 -0.67975543
+v 0.30069627 -0.42502537 -0.50177061
+v 0.14805553 -0.40552922 0.75686850
+v 0.15660512 -0.42945129 -0.71789146
+v -0.27148627 -0.42423816 -0.23051991
+v 0.28796462 -0.41959420 -0.16274475
+v -0.28421792 -0.41880698 0.10850594
+v -0.27576106 -0.41227712 0.50686007
+v -0.14862999 -0.40568496 0.86282971
+v -0.19941751 -0.42963818 -0.59073800
+v 0.20316421 -0.40876300 0.53649919
+v -0.33406067 0.09494372 0.71903570
+v -0.31287216 0.09604241 0.77836398
+v -0.46114537 0.07312558 -0.57349097
+v -0.45691694 0.07639052 -0.37431391
+v 0.32287592 0.09855128 0.68509810
+v -0.38066612 0.08948140 0.40120209
+vn 0.03867952 -0.99911976 -0.01623593
+vn 0.15930405 -0.98678095 0.02975867
+vn 0.16163142 -0.96738383 0.19504823
+vn -0.12372163 -0.99225804 -0.01081357
+vn -0.36125043 -0.92405788 0.12496062
+vn -0.24494187 -0.92934298 0.27626997
+vn 0.05986440 -0.07675815 -0.99525094
+vn -0.68893252 -0.06911985 -0.72152230
+vn -0.91832157 -0.12488814 -0.37561743
+vn 0.99111731 0.06827206 0.11412890
+vn 0.85375305 -0.05759622 0.51748276
+vn 0.45418873 0.10702261 0.88445393
+vn -0.97526825 -0.21704931 0.04173060
+vn -0.03692835 0.23327189 0.97171010
+vn 0.97740839 0.12998475 -0.16666374
+vn 0.66331920 -0.12460515 -0.73788969
+vn 0.24217923 0.90666598 0.34540704
+vn 0.01673177 0.90746071 0.41980366
+vn 0.31981763 0.91352107 -0.25138801
+vn 0.42774618 0.89813199 -0.10194176
+vn -0.24697599 0.96892276 -0.01383997
+vn 0.44644355 0.89248420 -0.06449886
+vn 0.04839782 0.94388992 -0.32669477
+vn -0.25779933 0.96584164 0.02625719
+vn -0.20169271 0.96740591 0.15312038
+vn 0.47817023 0.87762758 -0.03351220
+vn -0.24528562 0.94373575 -0.22180575
+vn -0.31520240 0.94278521 -0.10864384
+vn 0.53887915 0.84165165 0.03509642
+vn 0.52488104 0.83589770 0.16054572
+vn 0.20978202 0.04545470 0.97669103
+vn -0.98732849 -0.09706312 0.12554363
+vn 0.44630615 0.57834877 0.68287885
+vn 0.74142022 0.64477157 -0.18591845
+vn 0.46946808 0.71835340 -0.51338886
+vn 0.08695311 0.55455107 0.82759427
+vn 0.83826334 0.53463002 -0.10716961
+vn -0.09289972 0.53685165 -0.83854633
+vn 0.88775783 0.45700358 -0.05507958
+vn -0.53777720 0.84294757 -0.01533218
+vn 0.90456280 0.42629182 0.00643565
+vn -0.43015695 0.54629327 -0.71869929
+vn -0.56738682 0.81547809 0.11431396
+vn -0.36656809 0.73725853 0.56751889
+vn 0.78983390 0.47008338 0.39393405
+vn -0.58575761 0.77425753 -0.23961072
+vn 0.56226835 -0.14620880 0.81392708
+vn -0.94847711 -0.28799230 -0.13210450
+vn 0.94137165 -0.15558803 0.29935227
+vn -0.36161674 0.06926646 0.92975023
+vn 0.99364539 0.04658938 -0.10246109
+vn -0.59741945 0.79053850 -0.13468064
+vn -0.56565976 0.82288125 -0.05380965
+vn 0.92966571 0.36202705 -0.06825017
+vn -0.65032414 0.75952163 -0.01433244
+vn 0.96191797 0.26153542 0.07945461
+vn 0.95747383 0.23813284 0.16290061
+vn -0.80487614 0.57856395 0.13205362
+vn -0.55010983 0.53738568 0.63921499
+vn -0.61943362 -0.34427700 -0.70553195
+vn 0.84212724 -0.23425598 0.48574257
+vn -0.86838387 -0.39253097 -0.30303280
+vn -0.08877382 -0.34715060 -0.93359824
+vn -0.97456685 -0.21421881 0.06580090
+vn -0.97601095 -0.16325619 0.14404874
+vn -0.41983881 -0.15970579 0.89343686
+vn 0.90039706 -0.24161190 0.36181324
+vn 0.98923166 -0.12525685 -0.07570631
+vn 0.64722287 -0.27865825 -0.70954361
+vn 0.47719215 0.35685252 0.80308401
+vn 0.20498441 0.32866826 0.92193197
+vn -0.68408491 0.67380551 -0.27930981
+vn 0.05033073 0.27936495 -0.95886497
+vn -0.73290310 0.67804032 -0.05580662
+vn -0.80537039 0.59247607 0.01872543
+vn -0.93368873 0.34587194 0.09272516
+vn -0.58353743 0.39458081 -0.70978169
+vn 0.88165995 0.38451225 -0.27354352
+vn 0.72675249 0.40129120 -0.55749098
+vn -0.91197354 -0.40985374 -0.01800487
+vn 0.92394123 -0.34724769 0.16047321
+vn 0.89194597 -0.35034784 0.28581248
+vn -0.95269935 -0.24255164 0.18311924
+vn 0.92067529 -0.36458033 -0.13942094
+vn -0.91318393 0.21769648 0.34453352
+vn 0.97492242 0.22243204 0.00708933
+vn -0.83664485 0.52981123 -0.13901603
+vn -0.91690610 0.39789377 0.03104432
+vn -0.97883137 0.18936138 0.07766225
+vn 0.94208818 0.10397810 0.31883918
+vn -0.35654313 0.33039312 0.87390926
+vn -0.98761120 0.10848283 0.11338250
+vn 0.73864305 0.18817626 0.64729911
+vn 0.81493813 -0.47566640 0.33108506
+vn -0.58424222 -0.72476387 -0.36521002
+vn 0.65763337 -0.75291605 0.02521454
+vn -0.73974743 -0.65367287 -0.15964185
+vn 0.56902022 -0.48579253 0.66349198
+vn 0.13730444 -0.62830641 -0.76575358
+vn 0.65187302 -0.75531511 0.06753257
+vn 0.63945773 -0.75882094 0.12363084
+vn -0.89328716 -0.44929829 0.01300351
+vn 0.15667186 -0.43614217 0.88613426
+vn 0.74681672 -0.62750158 0.22024204
+vn -0.20225875 -0.71781045 -0.66621284
+vn -0.91668696 -0.39417646 0.06565006
+vn -0.32206350 -0.72120062 -0.61330642
+vn -0.85628468 -0.43866109 0.27267746
+vn -0.79343616 -0.43068570 0.43008009
+vn 0.74871039 -0.43640460 -0.49898274
+vn 0.86875939 -0.44871897 -0.20954331
+vn -0.78485077 0.13712117 0.60432363
+vn 0.09471409 0.09685195 -0.99078199
+vn -0.70264599 0.21020701 -0.67978057
+vn -0.93834752 0.24686325 -0.24199684
+vn -0.98962864 0.13660045 0.04444626
+vn -0.11210481 0.33224636 0.93650674
+vn 0.96267282 0.17942446 -0.20265218
+vn 0.70419611 0.10540349 -0.70213812
+vn -0.36644782 -0.92566082 -0.09417023
+vn 0.35053042 -0.86486892 -0.35934688
+vn 0.43002038 -0.89852081 -0.08799332
+vn 0.46243757 -0.79548701 0.39160172
+vn 0.11225328 -0.87569947 -0.46962712
+vn -0.51399412 -0.85758513 -0.01891540
+vn 0.40990959 -0.91147222 0.03453273
+vn -0.63941968 -0.76876111 0.01219916
+vn -0.73303789 -0.67192306 0.10571118
+vn -0.27685542 -0.72358484 0.63227847
+vn -0.25555454 -0.93924613 -0.22914750
+vn 0.49109354 -0.86315003 0.11746984
+vn -0.96568756 -0.06131610 0.25236455
+vn -0.88400147 -0.04516412 0.46529732
+vn -0.97788060 -0.09196171 -0.18786316
+vn -0.99565776 -0.07621673 0.05344753
+vn 0.95158393 -0.03430538 0.30546877
+vn -0.99193362 -0.04157605 0.11974610
+f 40//40 28//28 53//53
+f 109//109 133//133 108//108
+f 52//52 53//53 28//28
+f 74//74 53//53 52//52
+f 48//48 13//13 134//134
+f 135//135 134//134 13//13
+f 132//132 108//108 133//133
+f 101//101 100//100 81//81
+f 98//98 47//47 103//103
+f 98//98 103//103 129//129
+f 98//98 129//129 123//123
+f 23//23 35//35 38//38
+f 23//23 19//19 35//35
+f 84//84 68//68 81//81
+f 84//84 81//81 100//100
+f 21//21 28//28 40//40
+f 43//43 21//21 40//40
+f 43//43 24//24 21//21
+f 66//66 109//109 129//129
+f 66//66 129//129 103//103
+f 66//66 133//133 109//109
+f 66//66 31//31 50//50
+f 66//66 103//103 47//47
+f 66//66 47//47 31//31
+f 107//107 95//95 60//60
+f 107//107 130//130 95//95
+f 121//121 110//110 111//111
+f 121//121 111//111 122//122
+f 96//96 122//122 111//111
+f 96//96 100//100 122//122
+f 96//96 84//84 100//100
+f 96//96 111//111 84//84
+f 72//72 74//74 52//52
+f 72//72 87//87 74//74
+f 55//55 40//40 53//53
+f 55//55 53//53 74//74
+f 55//55 43//43 40//40
+f 55//55 58//58 43//43
+f 92//92 132//132 85//85
+f 92//92 137//137 132//132
+f 59//59 58//58 85//85
+f 59//59 43//43 58//58
+f 59//59 44//44 43//43
+f 76//76 88//88 89//89
+f 76//76 89//89 92//92
+f 76//76 92//92 85//85
+f 76//76 85//85 58//58
+f 75//75 76//76 58//58
+f 75//75 88//88 76//76
+f 75//75 55//55 74//74
+f 75//75 58//58 55//55
+f 75//75 74//74 87//87
+f 75//75 87//87 88//88
+f 112//112 50//50 91//91
+f 112//112 91//91 59//59
+f 112//112 66//66 50//50
+f 112//112 133//133 66//66
+f 112//112 59//59 85//85
+f 112//112 132//132 133//133
+f 112//112 85//85 132//132
+f 36//36 18//18 44//44
+f 62//62 95//95 97//97
+f 62//62 60//60 95//95
+f 32//32 137//137 135//135
+f 32//32 132//132 137//137
+f 32//32 65//65 132//132
+f 64//64 135//135 13//13
+f 64//64 32//32 135//135
+f 64//64 65//65 32//32
+f 83//83 106//106 108//108
+f 83//83 132//132 65//65
+f 83//83 108//108 132//132
+f 83//83 64//64 106//106
+f 83//83 65//65 64//64
+f 80//80 62//62 97//97
+f 80//80 48//48 62//62
+f 128//128 108//108 106//106
+f 128//128 5//5 6//6
+f 128//128 109//109 108//108
+f 128//128 6//6 129//129
+f 128//128 129//129 109//109
+f 104//104 101//101 81//81
+f 131//131 104//104 123//123
+f 131//131 101//101 104//104
+f 61//61 47//47 98//98
+f 61//61 11//11 47//47
+f 61//61 67//67 11//11
+f 136//136 67//67 49//49
+f 136//136 49//49 81//81
+f 136//136 10//10 90//90
+f 136//136 81//81 10//10
+f 136//136 93//93 11//11
+f 136//136 90//90 93//93
+f 136//136 11//11 67//67
+f 7//7 69//69 63//63
+f 7//7 16//16 69//69
+f 119//119 16//16 7//7
+f 119//119 7//7 113//113
+f 73//73 119//119 113//113
+f 73//73 35//35 79//79
+f 73//73 38//38 35//35
+f 73//73 79//79 119//119
+f 14//14 31//31 71//71
+f 14//14 50//50 31//31
+f 45//45 93//93 90//90
+f 12//12 31//31 47//47
+f 12//12 71//71 31//31
+f 12//12 47//47 11//11
+f 12//12 11//11 93//93
+f 25//25 44//44 18//18
+f 25//25 43//43 44//44
+f 25//25 24//24 43//43
+f 105//105 107//107 60//60
+f 105//105 60//60 63//63
+f 105//105 124//124 130//130
+f 105//105 130//130 107//107
+f 120//120 125//125 97//97
+f 120//120 95//95 130//130
+f 120//120 97//97 95//95
+f 120//120 130//130 4//4
+f 120//120 4//4 5//5
+f 120//120 5//5 125//125
+f 99//99 63//63 69//69
+f 99//99 121//121 124//124
+f 99//99 105//105 63//63
+f 99//99 124//124 105//105
+f 99//99 69//69 110//110
+f 99//99 110//110 121//121
+f 46//46 72//72 52//52
+f 46//46 52//52 28//28
+f 116//116 115//115 134//134
+f 116//116 134//134 135//135
+f 116//116 92//92 89//89
+f 116//116 89//89 88//88
+f 116//116 135//135 137//137
+f 116//116 137//137 92//92
+f 116//116 88//88 87//87
+f 116//116 87//87 115//115
+f 9//9 62//62 48//48
+f 9//9 48//48 134//134
+f 9//9 134//134 115//115
+f 9//9 8//8 60//60
+f 9//9 60//60 62//62
+f 9//9 115//115 114//114
+f 9//9 114//114 8//8
+f 102//102 80//80 97//97
+f 102//102 106//106 64//64
+f 102//102 64//64 13//13
+f 102//102 13//13 48//48
+f 102//102 48//48 80//80
+f 127//127 125//125 5//5
+f 127//127 5//5 128//128
+f 127//127 102//102 97//97
+f 127//127 97//97 125//125
+f 127//127 128//128 106//106
+f 127//127 106//106 102//102
+f 94//94 98//98 123//123
+f 94//94 123//123 104//104
+f 94//94 61//61 98//98
+f 94//94 67//67 61//61
+f 126//126 122//122 100//100
+f 126//126 2//2 122//122
+f 126//126 100//100 101//101
+f 126//126 101//101 131//131
+f 117//117 14//14 71//71
+f 117//117 71//71 36//36
+f 117//117 59//59 91//91
+f 117//117 91//91 50//50
+f 117//117 50//50 14//14
+f 117//117 36//36 44//44
+f 117//117 44//44 59//59
+f 57//57 45//45 90//90
+f 57//57 90//90 10//10
+f 70//70 12//12 93//93
+f 70//70 71//71 12//12
+f 70//70 33//33 36//36
+f 70//70 36//36 71//71
+f 70//70 93//93 45//45
+f 70//70 45//45 33//33
+f 1//1 4//4 130//130
+f 1//1 130//130 124//124
+f 1//1 124//124 121//121
+f 1//1 121//121 122//122
+f 1//1 122//122 2//2
+f 3//3 131//131 123//123
+f 3//3 126//126 131//131
+f 3//3 2//2 126//126
+f 3//3 123//123 129//129
+f 3//3 129//129 6//6
+f 3//3 1//1 2//2
+f 3//3 6//6 5//5
+f 3//3 5//5 4//4
+f 3//3 4//4 1//1
+f 42//42 46//46 28//28
+f 42//42 28//28 27//27
+f 42//42 23//23 38//38
+f 42//42 27//27 23//23
+f 82//82 67//67 94//94
+f 82//82 94//94 104//104
+f 82//82 104//104 81//81
+f 82//82 49//49 67//67
+f 82//82 81//81 49//49
+f 37//37 20//20 22//22
+f 51//51 10//10 81//81
+f 51//51 68//68 84//84
+f 51//51 81//81 68//68
+f 51//51 118//118 15//15
+f 51//51 84//84 111//111
+f 51//51 111//111 110//110
+f 51//51 110//110 69//69
+f 51//51 69//69 16//16
+f 51//51 16//16 119//119
+f 51//51 119//119 79//79
+f 51//51 79//79 78//78
+f 51//51 78//78 118//118
+f 86//86 51//51 15//15
+f 17//17 33//33 45//45
+f 17//17 45//45 30//30
+f 17//17 18//18 36//36
+f 17//17 36//36 33//33
+f 77//77 8//8 114//114
+f 77//77 114//114 115//115
+f 77//77 72//72 46//46
+f 77//77 46//46 42//42
+f 77//77 60//60 8//8
+f 77//77 115//115 87//87
+f 77//77 87//87 72//72
+f 77//77 63//63 60//60
+f 77//77 7//7 63//63
+f 77//77 113//113 7//7
+f 77//77 73//73 113//113
+f 77//77 42//42 38//38
+f 77//77 38//38 73//73
+f 34//34 37//37 78//78
+f 34//34 20//20 37//37
+f 34//34 19//19 20//20
+f 34//34 35//35 19//19
+f 34//34 78//78 79//79
+f 34//34 79//79 35//35
+f 56//56 41//41 29//29
+f 56//56 57//57 10//10
+f 56//56 29//29 30//30
+f 56//56 10//10 51//51
+f 56//56 51//51 86//86
+f 56//56 30//30 45//45
+f 56//56 45//45 57//57
+f 26//26 39//39 37//37
+f 26//26 22//22 20//20
+f 26//26 37//37 22//22
+f 26//26 29//29 41//41
+f 26//26 41//41 39//39
+f 26//26 17//17 30//30
+f 26//26 30//30 29//29
+f 26//26 20//20 19//19
+f 26//26 19//19 23//23
+f 26//26 23//23 27//27
+f 26//26 27//27 28//28
+f 26//26 28//28 21//21
+f 26//26 21//21 24//24
+f 26//26 24//24 25//25
+f 26//26 25//25 18//18
+f 26//26 18//18 17//17
+f 54//54 56//56 86//86
+f 54//54 86//86 15//15
+f 54//54 15//15 118//118
+f 54//54 118//118 78//78
+f 54//54 39//39 41//41
+f 54//54 41//41 56//56
+f 54//54 78//78 37//37
+f 54//54 37//37 39//39
\ No newline at end of file
diff --git a/apps/showcases/PegIntoSphereBed/PegIntoSphereBed.cpp b/apps/showcases/PegIntoSphereBed/PegIntoSphereBed.cpp
index 0a91f03c9..79df3f3eb 100644
--- a/apps/showcases/PegIntoSphereBed/PegIntoSphereBed.cpp
+++ b/apps/showcases/PegIntoSphereBed/PegIntoSphereBed.cpp
@@ -380,7 +380,7 @@ int main( int argc, char ** argv ) {
uint_t(1), vtk_out, "simulation_step");
vtkDomainOutput->write();
// mesapd mesh output
- mesa_pd::MeshParticleVTKOutput<mesh::PolyMesh> meshParticleVTK(ps, ss, "mesh", visSpacing);
+ mesa_pd::MeshParticleVTKOutput<mesh::PolyMesh> meshParticleVTK(ps, "mesh", visSpacing);
meshParticleVTK.addFaceOutput<mesa_pd::data::SelectParticleUid>("Uid");
auto surfaceVelocityDataSource = make_shared<mesa_pd::SurfaceVelocityVertexDataSource<mesh::PolyMesh,
mesa_pd::data::ParticleAccessorWithShape>>("SurfaceVelocity", ac);
@@ -415,7 +415,7 @@ int main( int argc, char ** argv ) {
}
// VTK
- meshParticleVTK();
+ meshParticleVTK(ac);
particleVtkWriter->write();
// Prepare Data Structures
diff --git a/python/mesa_pd.py b/python/mesa_pd.py
index 3161a1a27..adbd847e8 100755
--- a/python/mesa_pd.py
+++ b/python/mesa_pd.py
@@ -24,7 +24,13 @@ if __name__ == '__main__':
ps.add_property("force", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
ps.add_property("oldForce", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ON_OWNERSHIP_CHANGE")
+ # shape definition for cases with small number of different shapes
ps.add_property("shapeID", "size_t", defValue="", syncMode="ON_GHOST_CREATION")
+
+ # shape definition for cases with distinct shape per particle
+ ps.add_include("mesa_pd/data/shape/BaseShape.h")
+ ps.add_property("baseShape", "std::shared_ptr<walberla::mesa_pd::data::BaseShape>", defValue="make_shared<walberla::mesa_pd::data::BaseShape>()", syncMode="ON_GHOST_CREATION")
+
ps.add_property("rotation", "walberla::mesa_pd::Rot3", defValue="", syncMode="ALWAYS")
ps.add_property("angularVelocity", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
ps.add_property("torque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
@@ -49,6 +55,8 @@ if __name__ == '__main__':
ps.add_property("temperature", "walberla::real_t", defValue="real_t(0)", syncMode="ALWAYS")
ps.add_property("heatFlux", "walberla::real_t", defValue="real_t(0)", syncMode="NEVER")
+ ps.add_property("numContacts", "uint_t", defValue="0", syncMode="NEVER")
+
# Properties for HCSITS
ps.add_property("dv", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
ps.add_property("dw", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
@@ -144,6 +152,8 @@ if __name__ == '__main__':
hftn.add_property('hydrodynamicTorque', 'mesa_pd::Vec3', 'Vec3(real_t(0))')
hfn = mpd.add(mpi.PropertyNotification('HeatFluxNotification'))
hfn.add_property('heatFlux', 'real_t', 'real_t(0)')
+ ncn = mpd.add(mpi.PropertyNotification('NumContactNotification'))
+ ncn.add_property('numContacts', 'uint_t', '0')
mpd.add(mpi.ReduceContactHistory())
mpd.add(mpi.ReduceProperty())
mpd.add(mpi.ShapePackUnpack(ps))
diff --git a/python/mesa_pd/templates/mpi/ShapePackUnpack.templ.h b/python/mesa_pd/templates/mpi/ShapePackUnpack.templ.h
index 6ea9d158e..a1c9d79e4 100644
--- a/python/mesa_pd/templates/mpi/ShapePackUnpack.templ.h
+++ b/python/mesa_pd/templates/mpi/ShapePackUnpack.templ.h
@@ -61,6 +61,10 @@ namespace mpi {
buf >> shapeType;
switch (shapeType)
{
+ case BaseShape::INVALID_SHAPE :
+ bs = std::make_unique<mesa_pd::data::BaseShape>();
+ bs->unpack(buf);
+ break;
{%- for shape in particle.shapes %}
case {{shape}}::SHAPE_TYPE :
bs = std::make_unique<mesa_pd::data::{{shape}}>();
diff --git a/src/mesa_pd/data/ParticleAccessor.h b/src/mesa_pd/data/ParticleAccessor.h
index 14ed03dc2..f18c0105a 100644
--- a/src/mesa_pd/data/ParticleAccessor.h
+++ b/src/mesa_pd/data/ParticleAccessor.h
@@ -92,6 +92,10 @@ public:
size_t& getShapeIDRef(const size_t p_idx) {return ps_->getShapeIDRef(p_idx);}
void setShapeID(const size_t p_idx, size_t const & v) { ps_->setShapeID(p_idx, v);}
+ std::shared_ptr<walberla::mesa_pd::data::BaseShape> const & getBaseShape(const size_t p_idx) const {return ps_->getBaseShape(p_idx);}
+ std::shared_ptr<walberla::mesa_pd::data::BaseShape>& getBaseShapeRef(const size_t p_idx) {return ps_->getBaseShapeRef(p_idx);}
+ void setBaseShape(const size_t p_idx, std::shared_ptr<walberla::mesa_pd::data::BaseShape> const & v) { ps_->setBaseShape(p_idx, v);}
+
walberla::mesa_pd::Rot3 const & getRotation(const size_t p_idx) const {return ps_->getRotation(p_idx);}
walberla::mesa_pd::Rot3& getRotationRef(const size_t p_idx) {return ps_->getRotationRef(p_idx);}
void setRotation(const size_t p_idx, walberla::mesa_pd::Rot3 const & v) { ps_->setRotation(p_idx, v);}
@@ -140,6 +144,10 @@ public:
walberla::real_t& getHeatFluxRef(const size_t p_idx) {return ps_->getHeatFluxRef(p_idx);}
void setHeatFlux(const size_t p_idx, walberla::real_t const & v) { ps_->setHeatFlux(p_idx, v);}
+ uint_t const & getNumContacts(const size_t p_idx) const {return ps_->getNumContacts(p_idx);}
+ uint_t& getNumContactsRef(const size_t p_idx) {return ps_->getNumContactsRef(p_idx);}
+ void setNumContacts(const size_t p_idx, uint_t const & v) { ps_->setNumContacts(p_idx, v);}
+
walberla::mesa_pd::Vec3 const & getDv(const size_t p_idx) const {return ps_->getDv(p_idx);}
walberla::mesa_pd::Vec3& getDvRef(const size_t p_idx) {return ps_->getDvRef(p_idx);}
void setDv(const size_t p_idx, walberla::mesa_pd::Vec3 const & v) { ps_->setDv(p_idx, v);}
@@ -293,6 +301,10 @@ public:
void setShapeID(const size_t /*p_idx*/, size_t const & v) { shapeID_ = v;}
size_t& getShapeIDRef(const size_t /*p_idx*/) {return shapeID_;}
+ std::shared_ptr<walberla::mesa_pd::data::BaseShape> const & getBaseShape(const size_t /*p_idx*/) const {return baseShape_;}
+ void setBaseShape(const size_t /*p_idx*/, std::shared_ptr<walberla::mesa_pd::data::BaseShape> const & v) { baseShape_ = v;}
+ std::shared_ptr<walberla::mesa_pd::data::BaseShape>& getBaseShapeRef(const size_t /*p_idx*/) {return baseShape_;}
+
walberla::mesa_pd::Rot3 const & getRotation(const size_t /*p_idx*/) const {return rotation_;}
void setRotation(const size_t /*p_idx*/, walberla::mesa_pd::Rot3 const & v) { rotation_ = v;}
walberla::mesa_pd::Rot3& getRotationRef(const size_t /*p_idx*/) {return rotation_;}
@@ -341,6 +353,10 @@ public:
void setHeatFlux(const size_t /*p_idx*/, walberla::real_t const & v) { heatFlux_ = v;}
walberla::real_t& getHeatFluxRef(const size_t /*p_idx*/) {return heatFlux_;}
+ uint_t const & getNumContacts(const size_t /*p_idx*/) const {return numContacts_;}
+ void setNumContacts(const size_t /*p_idx*/, uint_t const & v) { numContacts_ = v;}
+ uint_t& getNumContactsRef(const size_t /*p_idx*/) {return numContacts_;}
+
walberla::mesa_pd::Vec3 const & getDv(const size_t /*p_idx*/) const {return dv_;}
void setDv(const size_t /*p_idx*/, walberla::mesa_pd::Vec3 const & v) { dv_ = v;}
walberla::mesa_pd::Vec3& getDvRef(const size_t /*p_idx*/) {return dv_;}
@@ -427,6 +443,7 @@ private:
walberla::mesa_pd::Vec3 force_;
walberla::mesa_pd::Vec3 oldForce_;
size_t shapeID_;
+ std::shared_ptr<walberla::mesa_pd::data::BaseShape> baseShape_;
walberla::mesa_pd::Rot3 rotation_;
walberla::mesa_pd::Vec3 angularVelocity_;
walberla::mesa_pd::Vec3 torque_;
@@ -439,6 +456,7 @@ private:
std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory> newContactHistory_;
walberla::real_t temperature_;
walberla::real_t heatFlux_;
+ uint_t numContacts_;
walberla::mesa_pd::Vec3 dv_;
walberla::mesa_pd::Vec3 dw_;
walberla::mesa_pd::Vec3 hydrodynamicForce_;
diff --git a/src/mesa_pd/data/ParticleStorage.h b/src/mesa_pd/data/ParticleStorage.h
index 4d64d3405..a25c9d29e 100644
--- a/src/mesa_pd/data/ParticleStorage.h
+++ b/src/mesa_pd/data/ParticleStorage.h
@@ -38,6 +38,7 @@
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/IAccessor.h>
#include <mesa_pd/data/Flags.h>
+#include <mesa_pd/data/shape/BaseShape.h>
#include <blockforest/BlockForest.h>
#include <mesa_pd/data/STLOverloads.h>
@@ -82,6 +83,7 @@ public:
using force_type = walberla::mesa_pd::Vec3;
using oldForce_type = walberla::mesa_pd::Vec3;
using shapeID_type = size_t;
+ using baseShape_type = std::shared_ptr<walberla::mesa_pd::data::BaseShape>;
using rotation_type = walberla::mesa_pd::Rot3;
using angularVelocity_type = walberla::mesa_pd::Vec3;
using torque_type = walberla::mesa_pd::Vec3;
@@ -94,6 +96,7 @@ public:
using newContactHistory_type = std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory>;
using temperature_type = walberla::real_t;
using heatFlux_type = walberla::real_t;
+ using numContacts_type = uint_t;
using dv_type = walberla::mesa_pd::Vec3;
using dw_type = walberla::mesa_pd::Vec3;
using hydrodynamicForce_type = walberla::mesa_pd::Vec3;
@@ -156,6 +159,10 @@ public:
shapeID_type& getShapeIDRef() {return storage_.getShapeIDRef(i_);}
void setShapeID(shapeID_type const & v) { storage_.setShapeID(i_, v);}
+ baseShape_type const & getBaseShape() const {return storage_.getBaseShape(i_);}
+ baseShape_type& getBaseShapeRef() {return storage_.getBaseShapeRef(i_);}
+ void setBaseShape(baseShape_type const & v) { storage_.setBaseShape(i_, v);}
+
rotation_type const & getRotation() const {return storage_.getRotation(i_);}
rotation_type& getRotationRef() {return storage_.getRotationRef(i_);}
void setRotation(rotation_type const & v) { storage_.setRotation(i_, v);}
@@ -204,6 +211,10 @@ public:
heatFlux_type& getHeatFluxRef() {return storage_.getHeatFluxRef(i_);}
void setHeatFlux(heatFlux_type const & v) { storage_.setHeatFlux(i_, v);}
+ numContacts_type const & getNumContacts() const {return storage_.getNumContacts(i_);}
+ numContacts_type& getNumContactsRef() {return storage_.getNumContactsRef(i_);}
+ void setNumContacts(numContacts_type const & v) { storage_.setNumContacts(i_, v);}
+
dv_type const & getDv() const {return storage_.getDv(i_);}
dv_type& getDvRef() {return storage_.getDvRef(i_);}
void setDv(dv_type const & v) { storage_.setDv(i_, v);}
@@ -339,6 +350,7 @@ public:
using force_type = walberla::mesa_pd::Vec3;
using oldForce_type = walberla::mesa_pd::Vec3;
using shapeID_type = size_t;
+ using baseShape_type = std::shared_ptr<walberla::mesa_pd::data::BaseShape>;
using rotation_type = walberla::mesa_pd::Rot3;
using angularVelocity_type = walberla::mesa_pd::Vec3;
using torque_type = walberla::mesa_pd::Vec3;
@@ -351,6 +363,7 @@ public:
using newContactHistory_type = std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory>;
using temperature_type = walberla::real_t;
using heatFlux_type = walberla::real_t;
+ using numContacts_type = uint_t;
using dv_type = walberla::mesa_pd::Vec3;
using dw_type = walberla::mesa_pd::Vec3;
using hydrodynamicForce_type = walberla::mesa_pd::Vec3;
@@ -413,6 +426,10 @@ public:
shapeID_type& getShapeIDRef(const size_t idx) {return shapeID_[idx];}
void setShapeID(const size_t idx, shapeID_type const & v) { shapeID_[idx] = v; }
+ baseShape_type const & getBaseShape(const size_t idx) const {return baseShape_[idx];}
+ baseShape_type& getBaseShapeRef(const size_t idx) {return baseShape_[idx];}
+ void setBaseShape(const size_t idx, baseShape_type const & v) { baseShape_[idx] = v; }
+
rotation_type const & getRotation(const size_t idx) const {return rotation_[idx];}
rotation_type& getRotationRef(const size_t idx) {return rotation_[idx];}
void setRotation(const size_t idx, rotation_type const & v) { rotation_[idx] = v; }
@@ -461,6 +478,10 @@ public:
heatFlux_type& getHeatFluxRef(const size_t idx) {return heatFlux_[idx];}
void setHeatFlux(const size_t idx, heatFlux_type const & v) { heatFlux_[idx] = v; }
+ numContacts_type const & getNumContacts(const size_t idx) const {return numContacts_[idx];}
+ numContacts_type& getNumContactsRef(const size_t idx) {return numContacts_[idx];}
+ void setNumContacts(const size_t idx, numContacts_type const & v) { numContacts_[idx] = v; }
+
dv_type const & getDv(const size_t idx) const {return dv_[idx];}
dv_type& getDvRef(const size_t idx) {return dv_[idx];}
void setDv(const size_t idx, dv_type const & v) { dv_[idx] = v; }
@@ -627,6 +648,7 @@ public:
std::vector<force_type> force_ {};
std::vector<oldForce_type> oldForce_ {};
std::vector<shapeID_type> shapeID_ {};
+ std::vector<baseShape_type> baseShape_ {};
std::vector<rotation_type> rotation_ {};
std::vector<angularVelocity_type> angularVelocity_ {};
std::vector<torque_type> torque_ {};
@@ -639,6 +661,7 @@ public:
std::vector<newContactHistory_type> newContactHistory_ {};
std::vector<temperature_type> temperature_ {};
std::vector<heatFlux_type> heatFlux_ {};
+ std::vector<numContacts_type> numContacts_ {};
std::vector<dv_type> dv_ {};
std::vector<dw_type> dw_ {};
std::vector<hydrodynamicForce_type> hydrodynamicForce_ {};
@@ -675,6 +698,7 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(const ParticleSt
getForceRef() = rhs.getForce();
getOldForceRef() = rhs.getOldForce();
getShapeIDRef() = rhs.getShapeID();
+ getBaseShapeRef() = rhs.getBaseShape();
getRotationRef() = rhs.getRotation();
getAngularVelocityRef() = rhs.getAngularVelocity();
getTorqueRef() = rhs.getTorque();
@@ -687,6 +711,7 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(const ParticleSt
getNewContactHistoryRef() = rhs.getNewContactHistory();
getTemperatureRef() = rhs.getTemperature();
getHeatFluxRef() = rhs.getHeatFlux();
+ getNumContactsRef() = rhs.getNumContacts();
getDvRef() = rhs.getDv();
getDwRef() = rhs.getDw();
getHydrodynamicForceRef() = rhs.getHydrodynamicForce();
@@ -720,6 +745,7 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(ParticleStorage:
getForceRef() = std::move(rhs.getForceRef());
getOldForceRef() = std::move(rhs.getOldForceRef());
getShapeIDRef() = std::move(rhs.getShapeIDRef());
+ getBaseShapeRef() = std::move(rhs.getBaseShapeRef());
getRotationRef() = std::move(rhs.getRotationRef());
getAngularVelocityRef() = std::move(rhs.getAngularVelocityRef());
getTorqueRef() = std::move(rhs.getTorqueRef());
@@ -732,6 +758,7 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(ParticleStorage:
getNewContactHistoryRef() = std::move(rhs.getNewContactHistoryRef());
getTemperatureRef() = std::move(rhs.getTemperatureRef());
getHeatFluxRef() = std::move(rhs.getHeatFluxRef());
+ getNumContactsRef() = std::move(rhs.getNumContactsRef());
getDvRef() = std::move(rhs.getDvRef());
getDwRef() = std::move(rhs.getDwRef());
getHydrodynamicForceRef() = std::move(rhs.getHydrodynamicForceRef());
@@ -766,6 +793,7 @@ void swap(ParticleStorage::Particle lhs, ParticleStorage::Particle rhs)
std::swap(lhs.getForceRef(), rhs.getForceRef());
std::swap(lhs.getOldForceRef(), rhs.getOldForceRef());
std::swap(lhs.getShapeIDRef(), rhs.getShapeIDRef());
+ std::swap(lhs.getBaseShapeRef(), rhs.getBaseShapeRef());
std::swap(lhs.getRotationRef(), rhs.getRotationRef());
std::swap(lhs.getAngularVelocityRef(), rhs.getAngularVelocityRef());
std::swap(lhs.getTorqueRef(), rhs.getTorqueRef());
@@ -778,6 +806,7 @@ void swap(ParticleStorage::Particle lhs, ParticleStorage::Particle rhs)
std::swap(lhs.getNewContactHistoryRef(), rhs.getNewContactHistoryRef());
std::swap(lhs.getTemperatureRef(), rhs.getTemperatureRef());
std::swap(lhs.getHeatFluxRef(), rhs.getHeatFluxRef());
+ std::swap(lhs.getNumContactsRef(), rhs.getNumContactsRef());
std::swap(lhs.getDvRef(), rhs.getDvRef());
std::swap(lhs.getDwRef(), rhs.getDwRef());
std::swap(lhs.getHydrodynamicForceRef(), rhs.getHydrodynamicForceRef());
@@ -812,6 +841,7 @@ std::ostream& operator<<( std::ostream& os, const ParticleStorage::Particle& p )
"force : " << p.getForce() << "\n" <<
"oldForce : " << p.getOldForce() << "\n" <<
"shapeID : " << p.getShapeID() << "\n" <<
+ "baseShape : " << p.getBaseShape() << "\n" <<
"rotation : " << p.getRotation() << "\n" <<
"angularVelocity : " << p.getAngularVelocity() << "\n" <<
"torque : " << p.getTorque() << "\n" <<
@@ -824,6 +854,7 @@ std::ostream& operator<<( std::ostream& os, const ParticleStorage::Particle& p )
"newContactHistory : " << p.getNewContactHistory() << "\n" <<
"temperature : " << p.getTemperature() << "\n" <<
"heatFlux : " << p.getHeatFlux() << "\n" <<
+ "numContacts : " << p.getNumContacts() << "\n" <<
"dv : " << p.getDv() << "\n" <<
"dw : " << p.getDw() << "\n" <<
"hydrodynamicForce : " << p.getHydrodynamicForce() << "\n" <<
@@ -928,6 +959,7 @@ inline ParticleStorage::iterator ParticleStorage::create(const id_t& uid)
force_.emplace_back(real_t(0));
oldForce_.emplace_back(real_t(0));
shapeID_.emplace_back();
+ baseShape_.emplace_back(make_shared<walberla::mesa_pd::data::BaseShape>());
rotation_.emplace_back();
angularVelocity_.emplace_back(real_t(0));
torque_.emplace_back(real_t(0));
@@ -940,6 +972,7 @@ inline ParticleStorage::iterator ParticleStorage::create(const id_t& uid)
newContactHistory_.emplace_back();
temperature_.emplace_back(real_t(0));
heatFlux_.emplace_back(real_t(0));
+ numContacts_.emplace_back(0);
dv_.emplace_back(real_t(0));
dw_.emplace_back(real_t(0));
hydrodynamicForce_.emplace_back(real_t(0));
@@ -999,6 +1032,7 @@ inline ParticleStorage::iterator ParticleStorage::erase(iterator& it)
force_.pop_back();
oldForce_.pop_back();
shapeID_.pop_back();
+ baseShape_.pop_back();
rotation_.pop_back();
angularVelocity_.pop_back();
torque_.pop_back();
@@ -1011,6 +1045,7 @@ inline ParticleStorage::iterator ParticleStorage::erase(iterator& it)
newContactHistory_.pop_back();
temperature_.pop_back();
heatFlux_.pop_back();
+ numContacts_.pop_back();
dv_.pop_back();
dw_.pop_back();
hydrodynamicForce_.pop_back();
@@ -1057,6 +1092,7 @@ inline void ParticleStorage::reserve(const size_t size)
force_.reserve(size);
oldForce_.reserve(size);
shapeID_.reserve(size);
+ baseShape_.reserve(size);
rotation_.reserve(size);
angularVelocity_.reserve(size);
torque_.reserve(size);
@@ -1069,6 +1105,7 @@ inline void ParticleStorage::reserve(const size_t size)
newContactHistory_.reserve(size);
temperature_.reserve(size);
heatFlux_.reserve(size);
+ numContacts_.reserve(size);
dv_.reserve(size);
dw_.reserve(size);
hydrodynamicForce_.reserve(size);
@@ -1100,6 +1137,7 @@ inline void ParticleStorage::clear()
force_.clear();
oldForce_.clear();
shapeID_.clear();
+ baseShape_.clear();
rotation_.clear();
angularVelocity_.clear();
torque_.clear();
@@ -1112,6 +1150,7 @@ inline void ParticleStorage::clear()
newContactHistory_.clear();
temperature_.clear();
heatFlux_.clear();
+ numContacts_.clear();
dv_.clear();
dw_.clear();
hydrodynamicForce_.clear();
@@ -1144,6 +1183,7 @@ inline size_t ParticleStorage::size() const
//WALBERLA_ASSERT_EQUAL( uid_.size(), force.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), oldForce.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), shapeID.size() );
+ //WALBERLA_ASSERT_EQUAL( uid_.size(), baseShape.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), rotation.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), angularVelocity.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), torque.size() );
@@ -1156,6 +1196,7 @@ inline size_t ParticleStorage::size() const
//WALBERLA_ASSERT_EQUAL( uid_.size(), newContactHistory.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), temperature.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), heatFlux.size() );
+ //WALBERLA_ASSERT_EQUAL( uid_.size(), numContacts.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), dv.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), dw.size() );
//WALBERLA_ASSERT_EQUAL( uid_.size(), hydrodynamicForce.size() );
@@ -1437,6 +1478,15 @@ public:
size_t const & operator()(const data::Particle& p) const {return p.getShapeID();}
};
///Predicate that selects a certain property from a Particle
+class SelectParticleBaseShape
+{
+public:
+ using return_type = std::shared_ptr<walberla::mesa_pd::data::BaseShape>;
+ std::shared_ptr<walberla::mesa_pd::data::BaseShape>& operator()(data::Particle& p) const {return p.getBaseShapeRef();}
+ std::shared_ptr<walberla::mesa_pd::data::BaseShape>& operator()(data::Particle&& p) const {return p.getBaseShapeRef();}
+ std::shared_ptr<walberla::mesa_pd::data::BaseShape> const & operator()(const data::Particle& p) const {return p.getBaseShape();}
+};
+///Predicate that selects a certain property from a Particle
class SelectParticleRotation
{
public:
@@ -1545,6 +1595,15 @@ public:
walberla::real_t const & operator()(const data::Particle& p) const {return p.getHeatFlux();}
};
///Predicate that selects a certain property from a Particle
+class SelectParticleNumContacts
+{
+public:
+ using return_type = uint_t;
+ uint_t& operator()(data::Particle& p) const {return p.getNumContactsRef();}
+ uint_t& operator()(data::Particle&& p) const {return p.getNumContactsRef();}
+ uint_t const & operator()(const data::Particle& p) const {return p.getNumContacts();}
+};
+///Predicate that selects a certain property from a Particle
class SelectParticleDv
{
public:
diff --git a/src/mesa_pd/data/shape/BaseShape.h b/src/mesa_pd/data/shape/BaseShape.h
index abd1eb511..4fbb58344 100644
--- a/src/mesa_pd/data/shape/BaseShape.h
+++ b/src/mesa_pd/data/shape/BaseShape.h
@@ -40,13 +40,14 @@ class BaseShape
public:
using ShapeTypeT = int;
+ BaseShape() = default;
explicit BaseShape(const int shapeType) : shapeType_(shapeType) {}
virtual ~BaseShape() = default;
///Updates mass and inertia according to the actual shape.
virtual void updateMassAndInertia(const real_t density);
- virtual real_t getVolume() const = 0;
+ virtual real_t getVolume() const {WALBERLA_ABORT("Not implemented!");}
const real_t& getMass() const {return mass_;}
const real_t& getInvMass() const {return invMass_;}
diff --git a/src/mesa_pd/kernel/InitParticlesForHCSITS.h b/src/mesa_pd/kernel/InitParticlesForHCSITS.h
index a4b62e657..2022a0b80 100644
--- a/src/mesa_pd/kernel/InitParticlesForHCSITS.h
+++ b/src/mesa_pd/kernel/InitParticlesForHCSITS.h
@@ -94,7 +94,7 @@ inline void InitParticlesForHCSITS::operator()(size_t j, Accessor& ac, real_t dt
* \param ac The particle accessor
* \param body The body whose velocities to time integrate
* \param dv On return the initial linear velocity correction.
- * \param dw On return the initial angular velocity correction.
+ * \param w On return the initial angular velocity correction.
* \param dt The time step size.
* \return void
*
diff --git a/src/mesa_pd/kernel/IntegrateParticlesHCSITS.h b/src/mesa_pd/kernel/IntegrateParticlesHCSITS.h
index d598e312a..0f00fc9a0 100644
--- a/src/mesa_pd/kernel/IntegrateParticlesHCSITS.h
+++ b/src/mesa_pd/kernel/IntegrateParticlesHCSITS.h
@@ -102,7 +102,6 @@ inline void IntegrateParticlesHCSITS::operator()(size_t j, PAccessor& ac, real_t
//*************************************************************************************************
/*!\brief Time integration of the position and orientation of a given body.
*
- * \param ac The particle accessor
* \param body The body whose position and orientation to time integrate
* \param v The linear velocity to use for time integration of the position.
* \param w The angular velocity to use for time integration of the orientation.
diff --git a/src/mesa_pd/mpi/ShapePackUnpack.h b/src/mesa_pd/mpi/ShapePackUnpack.h
index 4461369ec..8d4e991aa 100644
--- a/src/mesa_pd/mpi/ShapePackUnpack.h
+++ b/src/mesa_pd/mpi/ShapePackUnpack.h
@@ -64,6 +64,10 @@ namespace mpi {
buf >> shapeType;
switch (shapeType)
{
+ case BaseShape::INVALID_SHAPE :
+ bs = std::make_unique<mesa_pd::data::BaseShape>();
+ bs->unpack(buf);
+ break;
case Sphere::SHAPE_TYPE :
bs = std::make_unique<mesa_pd::data::Sphere>();
bs->unpack(buf);
diff --git a/src/mesa_pd/mpi/notifications/NumContactNotification.h b/src/mesa_pd/mpi/notifications/NumContactNotification.h
new file mode 100644
index 000000000..13eeb2097
--- /dev/null
+++ b/src/mesa_pd/mpi/notifications/NumContactNotification.h
@@ -0,0 +1,114 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//======================================================================================================================
+//
+// THIS FILE IS GENERATED - PLEASE CHANGE THE TEMPLATE !!!
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/mpi/notifications/NotificationType.h>
+#include <mesa_pd/mpi/notifications/reset.h>
+
+#include <core/mpi/Datatype.h>
+#include <core/mpi/RecvBuffer.h>
+#include <core/mpi/SendBuffer.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+/**
+ * Transmits force and torque information.
+ */
+class NumContactNotification
+{
+public:
+ struct Parameters
+ {
+ id_t uid_;
+ uint_t numContacts_;
+ };
+
+ inline explicit NumContactNotification( const data::Particle& p ) : p_(p) {}
+
+ const data::Particle& p_;
+};
+
+template <>
+void reset<NumContactNotification>(data::Particle& p)
+{
+ p.setNumContacts( 0 );
+}
+
+void reduce(data::Particle&& p, const NumContactNotification::Parameters& objparam)
+{
+ p.getNumContactsRef() += objparam.numContacts_;
+}
+
+void update(data::Particle&& p, const NumContactNotification::Parameters& objparam)
+{
+ p.setNumContacts( objparam.numContacts_ );
+}
+
+} // namespace mesa_pd
+} // namespace walberla
+
+//======================================================================================================================
+//
+// Send/Recv Buffer Serialization Specialization
+//
+//======================================================================================================================
+
+namespace walberla {
+namespace mpi {
+
+template< typename T, // Element type of SendBuffer
+ typename G> // Growth policy of SendBuffer
+mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, const mesa_pd::NumContactNotification& obj )
+{
+ buf.addDebugMarker( "pn" );
+ buf << obj.p_.getUid();
+ buf << obj.p_.getNumContacts();
+ return buf;
+}
+
+template< typename T> // Element type of RecvBuffer
+mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd::NumContactNotification::Parameters& objparam )
+{
+ buf.readDebugMarker( "pn" );
+ buf >> objparam.uid_;
+ buf >> objparam.numContacts_;
+ return buf;
+}
+
+template< >
+struct BufferSizeTrait< mesa_pd::NumContactNotification > {
+ static const bool constantSize = true;
+ static const uint_t size = BufferSizeTrait<id_t>::size +
+ BufferSizeTrait<uint_t>::size +
+ mpi::BUFFER_DEBUG_OVERHEAD;
+};
+
+} // mpi
+} // walberla
\ No newline at end of file
diff --git a/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h b/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
index f0f23574c..78de9c6ea 100644
--- a/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleCopyNotification.h
@@ -58,6 +58,7 @@ public:
walberla::real_t invMass {real_t(1)};
walberla::mesa_pd::Vec3 oldForce {real_t(0)};
size_t shapeID {};
+ std::shared_ptr<walberla::mesa_pd::data::BaseShape> baseShape {make_shared<walberla::mesa_pd::data::BaseShape>()};
walberla::mesa_pd::Rot3 rotation {};
walberla::mesa_pd::Vec3 angularVelocity {real_t(0)};
walberla::mesa_pd::Vec3 oldTorque {real_t(0)};
@@ -99,6 +100,7 @@ inline data::ParticleStorage::iterator createNewParticle(data::ParticleStorage&
pIt->setInvMass(data.invMass);
pIt->setOldForce(data.oldForce);
pIt->setShapeID(data.shapeID);
+ pIt->setBaseShape(data.baseShape);
pIt->setRotation(data.rotation);
pIt->setAngularVelocity(data.angularVelocity);
pIt->setOldTorque(data.oldTorque);
@@ -155,6 +157,7 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
buf << obj.particle_.getInvMass();
buf << obj.particle_.getOldForce();
buf << obj.particle_.getShapeID();
+ buf << obj.particle_.getBaseShape();
buf << obj.particle_.getRotation();
buf << obj.particle_.getAngularVelocity();
buf << obj.particle_.getOldTorque();
@@ -192,6 +195,7 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
buf >> objparam.invMass;
buf >> objparam.oldForce;
buf >> objparam.shapeID;
+ buf >> objparam.baseShape;
buf >> objparam.rotation;
buf >> objparam.angularVelocity;
buf >> objparam.oldTorque;
diff --git a/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h b/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h
index b9e6aa966..fa87d8571 100644
--- a/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h
+++ b/src/mesa_pd/mpi/notifications/ParticleGhostCopyNotification.h
@@ -56,6 +56,7 @@ public:
walberla::mesa_pd::Vec3 linearVelocity {real_t(0)};
walberla::real_t invMass {real_t(1)};
size_t shapeID {};
+ std::shared_ptr<walberla::mesa_pd::data::BaseShape> baseShape {make_shared<walberla::mesa_pd::data::BaseShape>()};
walberla::mesa_pd::Rot3 rotation {};
walberla::mesa_pd::Vec3 angularVelocity {real_t(0)};
walberla::real_t radiusAtTemperature {real_t(0)};
@@ -83,6 +84,7 @@ inline data::ParticleStorage::iterator createNewParticle(data::ParticleStorage&
pIt->setLinearVelocity(data.linearVelocity);
pIt->setInvMass(data.invMass);
pIt->setShapeID(data.shapeID);
+ pIt->setBaseShape(data.baseShape);
pIt->setRotation(data.rotation);
pIt->setAngularVelocity(data.angularVelocity);
pIt->setRadiusAtTemperature(data.radiusAtTemperature);
@@ -125,6 +127,7 @@ mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, cons
buf << obj.particle_.getLinearVelocity();
buf << obj.particle_.getInvMass();
buf << obj.particle_.getShapeID();
+ buf << obj.particle_.getBaseShape();
buf << obj.particle_.getRotation();
buf << obj.particle_.getAngularVelocity();
buf << obj.particle_.getRadiusAtTemperature();
@@ -148,6 +151,7 @@ mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, mesa_pd:
buf >> objparam.linearVelocity;
buf >> objparam.invMass;
buf >> objparam.shapeID;
+ buf >> objparam.baseShape;
buf >> objparam.rotation;
buf >> objparam.angularVelocity;
buf >> objparam.radiusAtTemperature;
diff --git a/src/mesa_pd/vtk/ConvexPolyhedron/MeshParticleVTKOutput.h b/src/mesa_pd/vtk/ConvexPolyhedron/MeshParticleVTKOutput.h
index b744ae730..97bd1e5cd 100644
--- a/src/mesa_pd/vtk/ConvexPolyhedron/MeshParticleVTKOutput.h
+++ b/src/mesa_pd/vtk/ConvexPolyhedron/MeshParticleVTKOutput.h
@@ -39,7 +39,7 @@ namespace mesa_pd {
template<typename MeshType>
class MeshParticleVTKOutput {
- static_assert(MeshType::IsPolyMesh == 1, "We need polygonal meshes here!");
+ static_assert(MeshType::IsPolyMesh == 1, "Polygonal meshes required for VTK output!");
public:
using ParticleSelectorFunc = std::function<bool (const walberla::mesa_pd::data::ParticleStorage::iterator& pIt)>;
@@ -47,16 +47,17 @@ public:
typedef typename mesh::DistributedVTKMeshWriter<MeshType>::Vertices Vertices;
MeshParticleVTKOutput( shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps,
- shared_ptr<walberla::mesa_pd::data::ShapeStorage> ss, const std::string & identifier,
+ const std::string & identifier,
const uint_t writeFrequency, const std::string & baseFolder = "vtk_out")
- : ps_(std::move(ps)), ss_(std::move(ss)), mesh_(make_shared<MeshType>()),
+ : ps_(std::move(ps)), mesh_(make_shared<MeshType>()),
faceToParticleIdxManager_(*mesh_, "particle"), vertexToParticleIdxManager_(*mesh_, "particle"),
meshWriter_(mesh_, identifier, writeFrequency, baseFolder) {
}
void setParticleSelector( const ParticleSelectorFunc& func) {particleSelector_ = func;}
- void assembleMesh();
+ template<typename AccessorWithShape_T>
+ void assembleMesh(AccessorWithShape_T & accessor);
template <typename Selector>
void addVertexOutput(const std::string& name);
@@ -68,8 +69,13 @@ public:
template <typename DataSourceType>
void addFaceDataSource(const shared_ptr<DataSourceType> & dataSource);
- void operator()() {
- assembleMesh();
+ template<typename AccessorWithShape_T>
+ void operator()(AccessorWithShape_T & accessor) {
+ static_assert(std::is_base_of<mesa_pd::data::IAccessor, AccessorWithShape_T>::value, "Provide a valid accessor as template");
+
+ if(meshWriter_.isWriteScheduled()) {
+ assembleMesh(accessor);
+ }
// the mesh writer writes the mesh to vtk files and adds properties as defined by the data sources
meshWriter_();
mesh_->clean(); // the output mesh is no longer needed, thus discard its contents
@@ -77,7 +83,6 @@ public:
private:
const shared_ptr<walberla::mesa_pd::data::ParticleStorage> ps_;
- const shared_ptr<walberla::mesa_pd::data::ShapeStorage> ss_;
shared_ptr<MeshType> mesh_; ///< the output mesh
// these "managers" (which are just maps basically) map the faces and vertices of the output mesh to particles
@@ -158,7 +163,10 @@ void MeshParticleVTKOutput<MeshType>::addFaceDataSource(const shared_ptr<DataSou
* \tparam MeshType
*/
template<typename MeshType>
-void MeshParticleVTKOutput<MeshType>::assembleMesh() {
+template<typename AccessorWithShape_T>
+void MeshParticleVTKOutput<MeshType>::assembleMesh(AccessorWithShape_T & accessor) {
+
+ static_assert(std::is_base_of<mesa_pd::data::IAccessor, AccessorWithShape_T>::value, "Provide a valid accessor as template");
// those will save the newly created vertices and faces for each mesh
// to make it possible to map a vertex/face to the corresponding particle
std::vector<typename MeshType::VertexHandle> newVertices;
@@ -171,13 +179,13 @@ void MeshParticleVTKOutput<MeshType>::assembleMesh() {
for (auto pIt = ps_->begin(); pIt != ps_->end(); ++pIt) {
if (!particleSelector_(pIt)) continue;
- auto& shape = ss_->shapes[pIt->getShapeID()];
+ data::BaseShape* shape = accessor.getShape(pIt->getIdx());
newVertices.clear();
newFaces.clear();
if (shape->getShapeType() == walberla::mesa_pd::data::ConvexPolyhedron::SHAPE_TYPE) {
- const auto& convexPolyhedron = *static_cast<walberla::mesa_pd::data::ConvexPolyhedron*>(shape.get());
+ const auto& convexPolyhedron = *static_cast<walberla::mesa_pd::data::ConvexPolyhedron*>(shape);
const auto& particle = *pIt;
// tessellate: add the shape at the particle's position into the output mesh
diff --git a/src/mesh_common/MeshOperations.h b/src/mesh_common/MeshOperations.h
index 979d0be61..159481e6b 100644
--- a/src/mesh_common/MeshOperations.h
+++ b/src/mesh_common/MeshOperations.h
@@ -67,6 +67,9 @@ Matrix3<typename MeshType::Scalar> computeInertiaTensor( const MeshType & mesh )
template< typename MeshType >
typename MeshType::Point computeCentroid( const MeshType & mesh, const typename MeshType::FaceHandle fh );
+template< typename MeshType >
+typename MeshType::Scalar computeBoundingSphereRadius( const MeshType & mesh, const typename MeshType::Point & referencePoint );
+
template< typename MeshType, typename InputIterator >
std::vector< typename MeshType::VertexHandle > findConnectedVertices( const MeshType & mesh, InputIterator facesBegin, InputIterator facesEnd );
@@ -334,6 +337,20 @@ typename MeshType::Point computeCentroid( const MeshType & mesh, const typename
return centroid;
}
+template< typename MeshType >
+typename MeshType::Scalar computeBoundingSphereRadius( const MeshType & mesh, const typename MeshType::Point & referencePoint )
+{
+ typedef typename MeshType::Scalar Scalar;
+
+ Scalar maxSqRadius(0);
+ for(auto vh : mesh.vertices()) {
+ auto v = mesh.point(vh);
+ auto refPointToVSqr = (v-referencePoint).sqrnorm();
+ maxSqRadius = std::max(maxSqRadius, refPointToVSqr );
+ }
+ return std::sqrt(maxSqRadius);
+}
+
/**
* \brief Computes all mass properties (mass, centroid, inertia tensor at once).
*
diff --git a/src/vtk/VTKTrait.h b/src/vtk/VTKTrait.h
index 1f451a622..7113b642d 100644
--- a/src/vtk/VTKTrait.h
+++ b/src/vtk/VTKTrait.h
@@ -126,5 +126,13 @@ struct VTKTrait<math::Rot3<T>>
constexpr static const uint_t components = 3;
};
+template <typename T, uint_t N>
+struct VTKTrait<std::array<T,N>>
+{
+ using type = typename VTKTrait<T>::type;
+ constexpr static char const * const type_string = VTKTrait<T>::type_string;
+ constexpr static const uint_t components = N;
+};
+
} // namespace vtk
} // namespace walberla
diff --git a/tests/mesa_pd/vtk/ConvexPolyhedronVTKOutput.cpp b/tests/mesa_pd/vtk/ConvexPolyhedronVTKOutput.cpp
index f734d0d0a..57af4a2d5 100644
--- a/tests/mesa_pd/vtk/ConvexPolyhedronVTKOutput.cpp
+++ b/tests/mesa_pd/vtk/ConvexPolyhedronVTKOutput.cpp
@@ -110,7 +110,7 @@ int main(int argc, char** argv) {
auto vtkDomainOutput = walberla::vtk::createVTKOutput_DomainDecomposition(forest, "domain_decomposition", 1, vtk_out, "simulation_step");
vtkDomainOutput->write();
- mesa_pd::MeshParticleVTKOutput< mesh::PolyMesh > meshParticleVTK(ps, ss, "mesh", uint_t(1));
+ mesa_pd::MeshParticleVTKOutput< mesh::PolyMesh > meshParticleVTK(ps, "mesh", uint_t(1));
meshParticleVTK.addFaceOutput< data::SelectParticleUid >("UID");
meshParticleVTK.addVertexOutput< data::SelectParticleInteractionRadius >("InteractionRadius");
meshParticleVTK.addFaceOutput< data::SelectParticleLinearVelocity >("LinearVelocity");
@@ -119,7 +119,7 @@ int main(int argc, char** argv) {
mesa_pd::SurfaceVelocityVertexDataSource< mesh::PolyMesh, mesa_pd::data::ParticleAccessorWithShape > >(
"SurfaceVelocity", *ac);
meshParticleVTK.addVertexDataSource(surfaceVelDataSource);
- meshParticleVTK();
+ meshParticleVTK(*ac);
return EXIT_SUCCESS;
}
--
GitLab