From a40b3c5c7481bfe4618e2ec3969f5c1020740d44 Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Mon, 17 Jun 2019 16:34:28 +0200
Subject: [PATCH] added particle sorting to ParticleStorage
---
.../templates/data/ParticleStorage.templ.h | 56 +++++-
.../HeatFluxNotification.templ.h | 7 +
src/mesa_pd/data/ParticleStorage.h | 74 ++++++-
.../mpi/notifications/HeatFluxNotification.h | 7 +
src/mesa_pd/sorting/HilbertCompareFunctor.cpp | 187 ++++++++++++++++++
src/mesa_pd/sorting/HilbertCompareFunctor.h | 72 +++++++
.../sorting/LinearizedCompareFunctor.cpp | 72 +++++++
.../sorting/LinearizedCompareFunctor.h | 71 +++++++
src/mesa_pd/vtk/ParticleVtkOutput.cpp | 4 +-
tests/mesa_pd/CMakeLists.txt | 3 +
tests/mesa_pd/Sorting.cpp | 89 +++++++++
11 files changed, 636 insertions(+), 6 deletions(-)
create mode 100644 src/mesa_pd/sorting/HilbertCompareFunctor.cpp
create mode 100644 src/mesa_pd/sorting/HilbertCompareFunctor.h
create mode 100644 src/mesa_pd/sorting/LinearizedCompareFunctor.cpp
create mode 100644 src/mesa_pd/sorting/LinearizedCompareFunctor.h
create mode 100644 tests/mesa_pd/Sorting.cpp
diff --git a/python/mesa_pd/templates/data/ParticleStorage.templ.h b/python/mesa_pd/templates/data/ParticleStorage.templ.h
index 7d6b1fa4b..c62e04564 100644
--- a/python/mesa_pd/templates/data/ParticleStorage.templ.h
+++ b/python/mesa_pd/templates/data/ParticleStorage.templ.h
@@ -133,7 +133,7 @@ public:
iterator begin() { return iterator(this, 0); }
iterator end() { return iterator(this, size()); }
- iterator operator[](const size_t n) { return iterator(this, n); }
+ Particle operator[](const size_t n) { return *iterator(this, n); }
{% for prop in properties %}
const {{prop.type}}& get{{prop.name | capFirst}}(const size_t idx) const {return {{prop.name}}_[idx];}
@@ -159,6 +159,8 @@ public:
inline void reserve(const size_t size);
inline void clear();
inline size_t size() const;
+ template <class Compare>
+ void sort(Compare comp);
/**
* Calls the provided functor \p func for all Particles.
@@ -256,6 +258,15 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(ParticleStorage:
return *this;
}
+inline
+void swap(ParticleStorage::Particle lhs, ParticleStorage::Particle rhs)
+{
+ if (lhs.i_ == rhs.i_) return;
+ {%- for prop in properties %}
+ std::swap(lhs.get{{prop.name | capFirst}}Ref(), rhs.get{{prop.name | capFirst}}Ref());
+ {%- endfor %}
+}
+
inline
std::ostream& operator<<( std::ostream& os, const ParticleStorage::Particle& p )
{
@@ -392,7 +403,7 @@ inline ParticleStorage::iterator ParticleStorage::find(const id_t& uid)
//use unordered_map for faster lookup
auto it = uidToIdx_.find(uid);
if (it == uidToIdx_.end()) return end();
- WALBERLA_ASSERT_EQUAL(it->first, uid, "Lookup via uidToIdx map is not up to date!!!");
+ WALBERLA_ASSERT_EQUAL(getUid(it->second), uid, "Lookup via uidToIdx map is not up to date!!!");
return iterator(this, it->second);
}
@@ -419,6 +430,47 @@ inline size_t ParticleStorage::size() const
return {{properties[0].name}}_.size();
}
+template <class Compare>
+void ParticleStorage::sort(Compare comp)
+{
+ using WeightPair = std::pair<size_t, double>; //idx, weight
+
+ const size_t length = size();
+ std::vector<size_t> newIdx(length); //where is old idx now?
+ std::vector<size_t> oldIdx(length); //what old idx is at that idx?
+ std::vector<WeightPair> weight(length);
+
+ for (size_t idx = 0; idx < length; ++idx)
+ {
+ newIdx[idx] = idx;
+ oldIdx[idx] = idx;
+ weight[idx] = std::make_pair(idx, comp.getWeight(operator[](idx)));
+ }
+ std::sort(weight.begin(), weight.end(), [](const WeightPair& lhs, const WeightPair& rhs){return lhs.second < rhs.second;});
+ for (size_t idx = 0; idx < length; ++idx)
+ {
+ using std::swap;
+ WALBERLA_ASSERT_IDENTICAL(weight[idx].second, comp.getWeight(operator[](newIdx[weight[idx].first])));
+
+ WALBERLA_ASSERT_LESS_EQUAL(comp.getWeight(operator[](newIdx[weight[idx].first])), comp.getWeight(operator[](idx)));
+ swap( operator[](idx), operator[](newIdx[weight[idx].first]) );
+ const auto lhsIDX = idx;
+ const auto rhsIDX = newIdx[weight[idx].first];
+
+ newIdx[oldIdx[lhsIDX]] = rhsIDX;
+ newIdx[oldIdx[rhsIDX]] = lhsIDX;
+
+ swap(oldIdx[lhsIDX], oldIdx[rhsIDX]);
+ }
+
+ //rebuild lookup table
+ uidToIdx_.clear();
+ for (size_t idx = 0; idx < length; ++idx)
+ {
+ uidToIdx_[getUid(idx)] = idx;
+ }
+}
+
{%- for const in ["", "const"] %}
template <typename Selector, typename Accessor, typename Func, typename... Args>
inline void ParticleStorage::forEachParticle(const bool openmp,
diff --git a/python/mesa_pd/templates/mpi/notifications/HeatFluxNotification.templ.h b/python/mesa_pd/templates/mpi/notifications/HeatFluxNotification.templ.h
index fa5caceda..ce463a214 100644
--- a/python/mesa_pd/templates/mpi/notifications/HeatFluxNotification.templ.h
+++ b/python/mesa_pd/templates/mpi/notifications/HeatFluxNotification.templ.h
@@ -29,6 +29,7 @@
#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>
@@ -54,6 +55,12 @@ public:
const data::Particle& p_;
};
+template <>
+void reset<HeatFluxNotification>(data::Particle& p)
+{
+ p.setHeatFlux( real_t(0) );
+}
+
void reduce(data::Particle&& p, const HeatFluxNotification::Parameters& objparam)
{
p.getHeatFluxRef() += objparam.heatFlux_;
diff --git a/src/mesa_pd/data/ParticleStorage.h b/src/mesa_pd/data/ParticleStorage.h
index e4d766318..aba34e3fe 100644
--- a/src/mesa_pd/data/ParticleStorage.h
+++ b/src/mesa_pd/data/ParticleStorage.h
@@ -223,7 +223,7 @@ public:
iterator begin() { return iterator(this, 0); }
iterator end() { return iterator(this, size()); }
- iterator operator[](const size_t n) { return iterator(this, n); }
+ Particle operator[](const size_t n) { return *iterator(this, n); }
const walberla::id_t& getUid(const size_t idx) const {return uid_[idx];}
@@ -341,6 +341,8 @@ public:
inline void reserve(const size_t size);
inline void clear();
inline size_t size() const;
+ template <class Compare>
+ void sort(Compare comp);
/**
* Calls the provided functor \p func for all Particles.
@@ -501,6 +503,33 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(ParticleStorage:
return *this;
}
+inline
+void swap(ParticleStorage::Particle lhs, ParticleStorage::Particle rhs)
+{
+ if (lhs.i_ == rhs.i_) return;
+ std::swap(lhs.getUidRef(), rhs.getUidRef());
+ std::swap(lhs.getPositionRef(), rhs.getPositionRef());
+ std::swap(lhs.getInteractionRadiusRef(), rhs.getInteractionRadiusRef());
+ std::swap(lhs.getFlagsRef(), rhs.getFlagsRef());
+ std::swap(lhs.getOwnerRef(), rhs.getOwnerRef());
+ std::swap(lhs.getGhostOwnersRef(), rhs.getGhostOwnersRef());
+ std::swap(lhs.getShapeIDRef(), rhs.getShapeIDRef());
+ std::swap(lhs.getRotationRef(), rhs.getRotationRef());
+ std::swap(lhs.getAngularVelocityRef(), rhs.getAngularVelocityRef());
+ std::swap(lhs.getTorqueRef(), rhs.getTorqueRef());
+ std::swap(lhs.getLinearVelocityRef(), rhs.getLinearVelocityRef());
+ std::swap(lhs.getInvMassRef(), rhs.getInvMassRef());
+ std::swap(lhs.getForceRef(), rhs.getForceRef());
+ std::swap(lhs.getOldForceRef(), rhs.getOldForceRef());
+ std::swap(lhs.getOldTorqueRef(), rhs.getOldTorqueRef());
+ std::swap(lhs.getTypeRef(), rhs.getTypeRef());
+ std::swap(lhs.getNextParticleRef(), rhs.getNextParticleRef());
+ std::swap(lhs.getOldContactHistoryRef(), rhs.getOldContactHistoryRef());
+ std::swap(lhs.getNewContactHistoryRef(), rhs.getNewContactHistoryRef());
+ std::swap(lhs.getTemperatureRef(), rhs.getTemperatureRef());
+ std::swap(lhs.getHeatFluxRef(), rhs.getHeatFluxRef());
+}
+
inline
std::ostream& operator<<( std::ostream& os, const ParticleStorage::Particle& p )
{
@@ -700,7 +729,7 @@ inline ParticleStorage::iterator ParticleStorage::find(const id_t& uid)
//use unordered_map for faster lookup
auto it = uidToIdx_.find(uid);
if (it == uidToIdx_.end()) return end();
- WALBERLA_ASSERT_EQUAL(it->first, uid, "Lookup via uidToIdx map is not up to date!!!");
+ WALBERLA_ASSERT_EQUAL(getUid(it->second), uid, "Lookup via uidToIdx map is not up to date!!!");
return iterator(this, it->second);
}
@@ -789,6 +818,47 @@ inline size_t ParticleStorage::size() const
//WALBERLA_ASSERT_EQUAL( uid_.size(), neighborState.size() );
return uid_.size();
}
+
+template <class Compare>
+void ParticleStorage::sort(Compare comp)
+{
+ using WeightPair = std::pair<size_t, double>; //idx, weight
+
+ const size_t length = size();
+ std::vector<size_t> newIdx(length); //where is old idx now?
+ std::vector<size_t> oldIdx(length); //what old idx is at that idx?
+ std::vector<WeightPair> weight(length);
+
+ for (size_t idx = 0; idx < length; ++idx)
+ {
+ newIdx[idx] = idx;
+ oldIdx[idx] = idx;
+ weight[idx] = std::make_pair(idx, comp.getWeight(operator[](idx)));
+ }
+ std::sort(weight.begin(), weight.end(), [](const WeightPair& lhs, const WeightPair& rhs){return lhs.second < rhs.second;});
+ for (size_t idx = 0; idx < length; ++idx)
+ {
+ using std::swap;
+ WALBERLA_ASSERT_IDENTICAL(weight[idx].second, comp.getWeight(operator[](newIdx[weight[idx].first])));
+
+ WALBERLA_ASSERT_LESS_EQUAL(comp.getWeight(operator[](newIdx[weight[idx].first])), comp.getWeight(operator[](idx)));
+ swap( operator[](idx), operator[](newIdx[weight[idx].first]) );
+ const auto lhsIDX = idx;
+ const auto rhsIDX = newIdx[weight[idx].first];
+
+ newIdx[oldIdx[lhsIDX]] = rhsIDX;
+ newIdx[oldIdx[rhsIDX]] = lhsIDX;
+
+ swap(oldIdx[lhsIDX], oldIdx[rhsIDX]);
+ }
+
+ //rebuild lookup table
+ uidToIdx_.clear();
+ for (size_t idx = 0; idx < length; ++idx)
+ {
+ uidToIdx_[getUid(idx)] = idx;
+ }
+}
template <typename Selector, typename Accessor, typename Func, typename... Args>
inline void ParticleStorage::forEachParticle(const bool openmp,
const Selector& selector,
diff --git a/src/mesa_pd/mpi/notifications/HeatFluxNotification.h b/src/mesa_pd/mpi/notifications/HeatFluxNotification.h
index 6fee2c541..884a684e9 100644
--- a/src/mesa_pd/mpi/notifications/HeatFluxNotification.h
+++ b/src/mesa_pd/mpi/notifications/HeatFluxNotification.h
@@ -29,6 +29,7 @@
#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>
@@ -54,6 +55,12 @@ public:
const data::Particle& p_;
};
+template <>
+void reset<HeatFluxNotification>(data::Particle& p)
+{
+ p.setHeatFlux( real_t(0) );
+}
+
void reduce(data::Particle&& p, const HeatFluxNotification::Parameters& objparam)
{
p.getHeatFluxRef() += objparam.heatFlux_;
diff --git a/src/mesa_pd/sorting/HilbertCompareFunctor.cpp b/src/mesa_pd/sorting/HilbertCompareFunctor.cpp
new file mode 100644
index 000000000..2cf2a8466
--- /dev/null
+++ b/src/mesa_pd/sorting/HilbertCompareFunctor.cpp
@@ -0,0 +1,187 @@
+//======================================================================================================================
+//
+// 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 HilbertCompareFunctor.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "HilbertCompareFunctor.h"
+
+#include <blockforest/HilbertCurveConstruction.h>
+#include <core/DataTypes.h>
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+#include <core/math/Vector3.h>
+#include <core/mpi/MPIManager.h>
+
+#include <stack>
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace sorting {
+
+math::Vector3<uint_t> getChildPosition(const uint_t currentSize, const math::Vector3<uint_t>& parent, const uint_t childId)
+{
+ auto size = currentSize >> 1;
+ auto child = parent;
+
+ if (!(childId & 0b001))
+ {
+ child[0] -= size;
+ }
+
+ if (!(childId & 0b010))
+ {
+ child[1] -= size;
+ }
+
+ if (!(childId & 0b100))
+ {
+ child[2] -= size;
+ }
+
+ return child;
+}
+
+std::vector< math::Vector3<uint_t> > generateHilbertCurve(const uint_t size)
+{
+ using Node = std::pair<uint_t, Vector3<uint_t>>; //level, coordinate
+
+ std::vector< math::Vector3<uint_t> > result;
+
+ std::stack< Node > stack;
+ std::stack< uint_t > orientation;
+
+ stack.push( std::make_pair(size, math::Vector3<uint_t>(size)) );
+ orientation.push( uint_t(0) );
+
+ while( !stack.empty() )
+ {
+ const Node node = stack.top();
+ uint_t hilbertIndex = orientation.top();
+
+ stack.pop();
+ orientation.pop();
+
+ if( node.first > 1 )
+ {
+ for( uint_t c = 8; c-- != 0; )
+ {
+ stack.push( std::make_pair( node.first >> 1,
+ getChildPosition(node.first, node.second, blockforest::hilbertOrder[hilbertIndex][c]) ));
+ orientation.push( blockforest::hilbertOrientation[hilbertIndex][c] );
+ }
+ }
+ else
+ {
+ result.push_back( node.second );
+ }
+ }
+ return result;
+}
+
+
+HilbertCompareFunctor::HilbertCompareFunctor(const math::AABB& domain, const uint_t cells)
+ : domain_(domain)
+ , cells_(cells)
+ , hilbertLookup_(cells_ * cells_ * cells_)
+{
+ WALBERLA_CHECK(math::uintIsPowerOfTwo(cells));
+ inverse_dx[0] = real_t(1.0) / (domain_.xSize() / real_c(cells_));
+ inverse_dx[1] = real_t(1.0) / (domain_.ySize() / real_c(cells_));
+ inverse_dx[2] = real_t(1.0) / (domain_.zSize() / real_c(cells_));
+ initializeLookup();
+}
+
+bool HilbertCompareFunctor::operator()(const data::Particle p1, const data::Particle p2) const
+{
+ const auto hash1 = discretize(p1.getPosition() - domain_.minCorner());
+ WALBERLA_ASSERT_LESS( hash1, cells_*cells_*cells_);
+ const auto hash2 = discretize(p2.getPosition() - domain_.minCorner());
+ WALBERLA_ASSERT_LESS( hash2, cells_*cells_*cells_);
+
+ return hilbertLookup_[hash1] < hilbertLookup_[hash2];
+}
+
+uint_t HilbertCompareFunctor::discretize(const Vec3& pos) const
+{
+ const real_t& x = pos[0];
+ const real_t& y = pos[1];
+ const real_t& z = pos[2];
+
+ const uint_t hashMask = (cells_ - 1);
+
+ math::Vector3<uint_t> hash;
+
+ if( x < 0 ) {
+ real_t i = ( -x ) * inverse_dx[0];
+ hash[0] = cells_ - 1 - ( static_cast<uint_t>( i ) & hashMask );
+ }
+ else {
+ real_t i = x * inverse_dx[0];
+ hash[0] = static_cast<uint_t>( i ) & hashMask;
+ }
+
+ if( y < 0 ) {
+ real_t i = ( -y ) * inverse_dx[1];
+ hash[1] = cells_ - 1 - ( static_cast<uint_t>( i ) & hashMask );
+ }
+ else {
+ real_t i = y * inverse_dx[1];
+ hash[1] = static_cast<uint_t>( i ) & hashMask;
+ }
+
+ if( z < 0 ) {
+ real_t i = ( -z ) * inverse_dx[2];
+ hash[2] = cells_ - 1 - ( static_cast<uint_t>( i ) & hashMask );
+ }
+ else {
+ real_t i = z * inverse_dx[2];
+ hash[2] = static_cast<uint_t>( i ) & hashMask;
+ }
+
+ WALBERLA_ASSERT_LESS(hash[0], cells_);
+ WALBERLA_ASSERT_LESS(hash[1], cells_);
+ WALBERLA_ASSERT_LESS(hash[2], cells_);
+
+ return hash[2] * cells_ * cells_ + hash[1] * cells_ + hash[0];
+}
+
+
+void HilbertCompareFunctor::initializeLookup()
+{
+ const auto sfc = generateHilbertCurve(cells_);
+ uint_t counter = 0;
+ for (auto p : sfc)
+ {
+ WALBERLA_ASSERT_GREATER(p[0], 0);
+ WALBERLA_ASSERT_GREATER(p[1], 0);
+ WALBERLA_ASSERT_GREATER(p[2], 0);
+ p[0] -= 1;
+ p[1] -= 1;
+ p[2] -= 1;
+
+ hilbertLookup_[p[2] * cells_ * cells_ + p[1] * cells_ + p[0]] = counter;
+ ++counter;
+ }
+ WALBERLA_ASSERT_EQUAL(counter, sfc.size());
+ WALBERLA_ASSERT_EQUAL(counter, hilbertLookup_.size());
+}
+
+} //namespace sorting
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/src/mesa_pd/sorting/HilbertCompareFunctor.h b/src/mesa_pd/sorting/HilbertCompareFunctor.h
new file mode 100644
index 000000000..70e91665f
--- /dev/null
+++ b/src/mesa_pd/sorting/HilbertCompareFunctor.h
@@ -0,0 +1,72 @@
+//======================================================================================================================
+//
+// 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 HilbertCompareFunctor.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <core/math/AABB.h>
+#include <core/math/Vector3.h>
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/ParticleStorage.h>
+
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace sorting {
+
+/**
+ * Compare functor which sorts particles according to the hilbert space filling curve.
+ */
+class HilbertCompareFunctor
+{
+public:
+ /**
+ * Subdivides a domain into cells and arranges the cells according to the sfc.
+ * The position within the sfc gives the particles their weight.
+ * \param domain domain which should be partitioned into cells
+ * \param cells number of cells in every direction.
+ * \attention cells has to be a power of 2!
+ */
+ HilbertCompareFunctor(const math::AABB& domain, const uint_t cells);
+ bool operator()(const data::Particle bd1, const data::Particle bd2) const;
+ double getWeight(const data::Particle p1) const;
+private:
+ uint_t discretize(const Vec3& pos) const;
+ void initializeLookup();
+
+ math::AABB domain_;
+ const uint_t cells_;
+ Vec3 inverse_dx;
+ std::vector<uint_t> hilbertLookup_;
+};
+
+inline
+double HilbertCompareFunctor::getWeight(const data::Particle p1) const
+{
+ const auto hash1 = discretize(p1.getPosition() - domain_.minCorner());
+ WALBERLA_ASSERT_LESS( hash1, cells_*cells_*cells_);
+
+ return double_c(hilbertLookup_[hash1]);
+}
+
+} //namespace sorting
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/src/mesa_pd/sorting/LinearizedCompareFunctor.cpp b/src/mesa_pd/sorting/LinearizedCompareFunctor.cpp
new file mode 100644
index 000000000..7588285a4
--- /dev/null
+++ b/src/mesa_pd/sorting/LinearizedCompareFunctor.cpp
@@ -0,0 +1,72 @@
+//======================================================================================================================
+//
+// 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 LinearizedCompareFunctor.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "LinearizedCompareFunctor.h"
+
+#include <core/DataTypes.h>
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+#include <core/math/Vector3.h>
+#include <core/mpi/MPIManager.h>
+
+#include <stack>
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace sorting {
+
+LinearizedCompareFunctor::LinearizedCompareFunctor(const math::AABB& domain, const Vector3<uint_t> cells)
+ : domain_(domain)
+ , cells_(cells)
+{
+ inverse_dx[0] = real_c(cells_[0]) / domain_.xSize();
+ inverse_dx[1] = real_c(cells_[1]) / domain_.ySize();
+ inverse_dx[2] = real_c(cells_[2]) / domain_.zSize();
+}
+
+bool LinearizedCompareFunctor::operator()(const data::Particle p1, const data::Particle p2) const
+{
+ const auto hash1 = discretize(p1.getPosition() - domain_.minCorner());
+ WALBERLA_ASSERT_LESS( hash1, cells_[0]*cells_[1]*cells_[2]);
+ const auto hash2 = discretize(p2.getPosition() - domain_.minCorner());
+ WALBERLA_ASSERT_LESS( hash2, cells_[0]*cells_[1]*cells_[2]);
+
+ return hash1 < hash2;
+}
+
+uint_t LinearizedCompareFunctor::discretize(const Vec3& pos) const
+{
+ int x = int_c(pos[0] * inverse_dx[0]);
+ int y = int_c(pos[1] * inverse_dx[1]);
+ int z = int_c(pos[2] * inverse_dx[2]);
+ if (x<0) x=0;
+ if (y<0) y=0;
+ if (z<0) z=0;
+ if (x>=int_c(cells_[0])) x=int_c(cells_[0])-1;
+ if (y>=int_c(cells_[1])) y=int_c(cells_[1])-1;
+ if (z>=int_c(cells_[2])) z=int_c(cells_[2])-1;
+
+ return uint_c(z) * cells_[1] * cells_[0] + uint_c(y) * cells_[0] + uint_c(x);
+}
+
+} //namespace sorting
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/src/mesa_pd/sorting/LinearizedCompareFunctor.h b/src/mesa_pd/sorting/LinearizedCompareFunctor.h
new file mode 100644
index 000000000..9caf02bd0
--- /dev/null
+++ b/src/mesa_pd/sorting/LinearizedCompareFunctor.h
@@ -0,0 +1,71 @@
+//======================================================================================================================
+//
+// 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 LinearizedCompareFunctor.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <core/math/AABB.h>
+#include <core/math/Vector3.h>
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/ParticleStorage.h>
+
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace sorting {
+
+/**
+ * Compare functor which sorts particles along the x, y and z axes.
+ */
+class LinearizedCompareFunctor
+{
+public:
+ /**
+ * Subdivides a domain into cells and arranges the cells in a linearized fashion.
+ * The position within the linearization gives the particles their weight.
+ * \param domain domain which should be partitioned into cells
+ * \param cells number of cells in every direction.
+ */
+ LinearizedCompareFunctor(const math::AABB& domain, const Vector3<uint_t> cells);
+ bool operator()(const data::Particle bd1, const data::Particle bd2) const;
+ double getWeight(const data::Particle p1) const;
+private:
+ uint_t discretize(const Vec3& pos) const;
+ void initializeLookup();
+
+ math::AABB domain_;
+ const Vector3<uint_t> cells_;
+ Vec3 inverse_dx;
+ std::vector<uint_t> hilbertLookup_;
+};
+
+inline
+double LinearizedCompareFunctor::getWeight(const data::Particle p1) const
+{
+ const auto hash = discretize(p1.getPosition() - domain_.minCorner());
+ WALBERLA_ASSERT_LESS( hash, cells_[0]*cells_[1]*cells_[2]);
+
+ return double_c(hash);
+}
+
+} //namespace sorting
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/src/mesa_pd/vtk/ParticleVtkOutput.cpp b/src/mesa_pd/vtk/ParticleVtkOutput.cpp
index 0210b4e29..b63295893 100644
--- a/src/mesa_pd/vtk/ParticleVtkOutput.cpp
+++ b/src/mesa_pd/vtk/ParticleVtkOutput.cpp
@@ -60,7 +60,7 @@ void ParticleVtkOutput::push( std::ostream& os, const uint_t data, const uint_t
WALBERLA_ASSERT_LESS( point, particleIndices_.size() );
WALBERLA_ASSERT_LESS( particleIndices_[point], ps_->size() );
- selectors_[data].second->push(os, *(*ps_)[particleIndices_[point]], component);
+ selectors_[data].second->push(os, (*ps_)[particleIndices_[point]], component);
}
void ParticleVtkOutput::push( walberla::vtk::Base64Writer& b64, const uint_t data, const uint_t point, const uint_t component )
@@ -69,7 +69,7 @@ void ParticleVtkOutput::push( walberla::vtk::Base64Writer& b64, const uint_t dat
WALBERLA_ASSERT_LESS( point, particleIndices_.size() );
WALBERLA_ASSERT_LESS( particleIndices_[point], ps_->size() );
- selectors_[data].second->push(b64, *(*ps_)[particleIndices_[point]], component);
+ selectors_[data].second->push(b64, (*ps_)[particleIndices_[point]], component);
}
void ParticleVtkOutput::addOutput(const std::string& name, std::shared_ptr<IOutputSelector> selector)
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index ae8501523..1938034f1 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -128,6 +128,9 @@ waLBerla_execute_test( NAME MESA_PD_MPI_ReduceContactHistory PROCESSES 8 )
waLBerla_compile_test( NAME MESA_PD_MPI_ReduceProperty FILES mpi/ReduceProperty.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_MPI_ReduceProperty PROCESSES 8 )
+waLBerla_compile_test( NAME MESA_PD_Sorting FILES Sorting.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_Sorting )
+
waLBerla_compile_test( NAME MESA_PD_VTK_Outputs FILES vtk/VTKOutputs.cpp DEPENDS blockforest core vtk )
waLBerla_execute_test( NAME MESA_PD_VTK_Outputs PROCESSES 8 )
diff --git a/tests/mesa_pd/Sorting.cpp b/tests/mesa_pd/Sorting.cpp
new file mode 100644
index 000000000..d4408cd71
--- /dev/null
+++ b/tests/mesa_pd/Sorting.cpp
@@ -0,0 +1,89 @@
+//======================================================================================================================
+//
+// 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 Sorting.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <mesa_pd/sorting/HilbertCompareFunctor.h>
+#include <mesa_pd/sorting/LinearizedCompareFunctor.h>
+#include <mesa_pd/data/ParticleStorage.h>
+
+#include <core/Abort.h>
+#include <core/Environment.h>
+#include <core/math/Random.h>
+#include <core/grid_generator/SCIterator.h>
+#include <core/logging/Logging.h>
+
+namespace walberla {
+namespace mesa_pd {
+
+int main()
+{
+ auto ps = std::make_shared<data::ParticleStorage>(100);
+
+ const math::AABB domain(real_t(0), real_t(0), real_t(0),
+ real_t(2), real_t(2), real_t(2));
+ for (auto pt : grid_generator::SCGrid(domain, Vec3(real_c(0.5)), real_t(1)))
+ {
+ auto p = ps->create();
+ p->getPositionRef() = pt;
+ }
+
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[0].getPosition(), Vec3(real_t(0.5), real_t(0.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[1].getPosition(), Vec3(real_t(1.5), real_t(0.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[2].getPosition(), Vec3(real_t(0.5), real_t(1.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[3].getPosition(), Vec3(real_t(1.5), real_t(1.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[4].getPosition(), Vec3(real_t(0.5), real_t(0.5), real_t(1.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[5].getPosition(), Vec3(real_t(1.5), real_t(0.5), real_t(1.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[6].getPosition(), Vec3(real_t(0.5), real_t(1.5), real_t(1.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[7].getPosition(), Vec3(real_t(1.5), real_t(1.5), real_t(1.5)));
+
+ sorting::HilbertCompareFunctor hilbert(domain, 2);
+ ps->sort(hilbert);
+
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[0].getPosition(), Vec3(real_t(0.5), real_t(0.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[1].getPosition(), Vec3(real_t(1.5), real_t(0.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[2].getPosition(), Vec3(real_t(1.5), real_t(1.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[3].getPosition(), Vec3(real_t(0.5), real_t(1.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[4].getPosition(), Vec3(real_t(0.5), real_t(1.5), real_t(1.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[5].getPosition(), Vec3(real_t(1.5), real_t(1.5), real_t(1.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[6].getPosition(), Vec3(real_t(1.5), real_t(0.5), real_t(1.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[7].getPosition(), Vec3(real_t(0.5), real_t(0.5), real_t(1.5)));
+
+ sorting::LinearizedCompareFunctor linear(domain, Vector3<uint_t>(2,2,2));
+ ps->sort(linear);
+
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[0].getPosition(), Vec3(real_t(0.5), real_t(0.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[1].getPosition(), Vec3(real_t(1.5), real_t(0.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[2].getPosition(), Vec3(real_t(0.5), real_t(1.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[3].getPosition(), Vec3(real_t(1.5), real_t(1.5), real_t(0.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[4].getPosition(), Vec3(real_t(0.5), real_t(0.5), real_t(1.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[5].getPosition(), Vec3(real_t(1.5), real_t(0.5), real_t(1.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[6].getPosition(), Vec3(real_t(0.5), real_t(1.5), real_t(1.5)));
+ WALBERLA_CHECK_FLOAT_EQUAL((*ps)[7].getPosition(), Vec3(real_t(1.5), real_t(1.5), real_t(1.5)));
+
+ return EXIT_SUCCESS;
+}
+
+} // namespace mesa_pd
+} // namespace walberla
+
+int main( int argc, char* argv[] )
+{
+ walberla::Environment env(argc, argv);
+ return walberla::mesa_pd::main();
+}
--
GitLab