diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 78b72bcafc640ecc888c0997eeb3ff8f3056db58..d03051cae50a4f2c39f95fbd9f9fb8e8cb06cdf9 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -43,6 +43,7 @@ add_library(
MeshTransformer.cpp
MeshUtils.cpp
MeshVariant.cpp
+ StencilBuilder.cpp
SynchronizerManager.cpp
)
diff --git a/src/mesh/Stencil.hpp b/src/mesh/Stencil.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d67aa09ed03252b5085c9238d4e968c5947c2678
--- /dev/null
+++ b/src/mesh/Stencil.hpp
@@ -0,0 +1,26 @@
+#ifndef STENCIL_HPP
+#define STENCIL_HPP
+
+#include <mesh/ConnectivityMatrix.hpp>
+#include <mesh/ItemId.hpp>
+
+class Stencil
+{
+ private:
+ ConnectivityMatrix m_stencil;
+
+ public:
+ PUGS_INLINE
+ auto
+ operator[](CellId cell_id) const
+ {
+ return m_stencil[cell_id];
+ }
+
+ Stencil(const ConnectivityMatrix& stencil) : m_stencil(stencil) {}
+ Stencil(const Stencil&) = default;
+ Stencil(Stencil&&) = default;
+ ~Stencil() = default;
+};
+
+#endif // STENCIL_HPP
diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a43dbef60cd76a280b8b3ab679c2d0c4e5b5406f
--- /dev/null
+++ b/src/mesh/StencilBuilder.cpp
@@ -0,0 +1,117 @@
+#include <mesh/Connectivity.hpp>
+#include <mesh/StencilBuilder.hpp>
+
+template <typename ConnectivityType>
+Array<const uint32_t>
+StencilBuilder::_getRowMap(const ConnectivityType& connectivity) const
+{
+ auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+
+ Array<uint32_t> row_map{connectivity.numberOfCells() + 1};
+ row_map[0] = 0;
+ std::vector<CellId> neighbors;
+ for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+ neighbors.resize(0);
+ // The stencil is not built for ghost cells
+ if (cell_is_owned[cell_id]) {
+ auto cell_nodes = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ const NodeId node_id = cell_nodes[i_node];
+ auto node_cells = node_to_cell_matrix[node_id];
+ for (size_t i_node_cell = 0; i_node_cell < node_cells.size(); ++i_node_cell) {
+ const CellId node_cell_id = node_cells[i_node_cell];
+ if (node_cell_id != cell_id) {
+ neighbors.push_back(node_cells[i_node_cell]);
+ }
+ }
+ }
+ std::sort(neighbors.begin(), neighbors.end());
+ neighbors.erase(std::unique(neighbors.begin(), neighbors.end()), neighbors.end());
+ }
+ // The cell itself is not counted
+ row_map[cell_id + 1] = row_map[cell_id] + neighbors.size();
+ }
+
+ return row_map;
+}
+
+template <typename ConnectivityType>
+Array<const uint32_t>
+StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Array<const uint32_t>& row_map) const
+{
+ auto cell_number = connectivity.cellNumber();
+
+ Array<uint32_t> max_index(row_map.size() - 1);
+ parallel_for(
+ max_index.size(), PUGS_LAMBDA(size_t i) { max_index[i] = row_map[i]; });
+
+ auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+
+ Array<uint32_t> column_indices(row_map[row_map.size() - 1]);
+ column_indices.fill(std::numeric_limits<uint32_t>::max());
+
+ for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+ // The stencil is not built for ghost cells
+ if (cell_is_owned[cell_id]) {
+ auto cell_nodes = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ const NodeId node_id = cell_nodes[i_node];
+ auto node_cells = node_to_cell_matrix[node_id];
+ for (size_t i_node_cell = 0; i_node_cell < node_cells.size(); ++i_node_cell) {
+ const CellId node_cell_id = node_cells[i_node_cell];
+ if (node_cell_id != cell_id) {
+ bool found = false;
+ for (size_t i_index = row_map[cell_id]; i_index < max_index[cell_id]; ++i_index) {
+ if (column_indices[i_index] == node_cell_id) {
+ found = true;
+ break;
+ }
+ }
+ if (not found) {
+ int node_cell_number = cell_number[node_cell_id];
+ size_t i_index = row_map[cell_id];
+ // search for position for index
+ while ((i_index < max_index[cell_id])) {
+ if (node_cell_number > cell_number[CellId(column_indices[i_index])]) {
+ ++i_index;
+ } else {
+ break;
+ }
+ }
+
+ for (size_t i_destination = max_index[cell_id]; i_destination > i_index; --i_destination) {
+ size_t i_source = i_destination - 1;
+
+ column_indices[i_destination] = column_indices[i_source];
+ }
+ ++max_index[cell_id];
+ column_indices[i_index] = node_cell_id;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ return column_indices;
+}
+
+template <typename ConnectivityType>
+Stencil
+StencilBuilder::build(const ConnectivityType& connectivity) const
+{
+ Array<const uint32_t> row_map = this->_getRowMap(connectivity);
+ Array<const uint32_t> column_indices = this->_getColumnIndices(connectivity, row_map);
+
+ return ConnectivityMatrix{row_map, column_indices};
+}
+
+template Stencil StencilBuilder::build(const Connectivity<1>& connectivity) const;
+template Stencil StencilBuilder::build(const Connectivity<2>& connectivity) const;
+template Stencil StencilBuilder::build(const Connectivity<3>& connectivity) const;
diff --git a/src/mesh/StencilBuilder.hpp b/src/mesh/StencilBuilder.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3a42350795f0caf8e749bb6235874ea3fdb65eed
--- /dev/null
+++ b/src/mesh/StencilBuilder.hpp
@@ -0,0 +1,26 @@
+#ifndef STENCIL_BUILDER_HPP
+#define STENCIL_BUILDER_HPP
+
+#include <mesh/Stencil.hpp>
+
+class StencilBuilder
+{
+ private:
+ template <typename ConnectivityType>
+ Array<const uint32_t> _getRowMap(const ConnectivityType& connectivity) const;
+
+ template <typename ConnectivityType>
+ Array<const uint32_t> _getColumnIndices(const ConnectivityType& connectivity,
+ const Array<const uint32_t>& row_map) const;
+
+ public:
+ template <typename ConnectivityType>
+ Stencil build(const ConnectivityType& connectivity) const;
+
+ StencilBuilder() = default;
+ StencilBuilder(const StencilBuilder&) = default;
+ StencilBuilder(StencilBuilder&&) = default;
+ ~StencilBuilder() = default;
+};
+
+#endif // STENCIL_BUILDER_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index e7c75245c6c40f9f23b61d000d379c0211a62dda..68da0dc88f3ab393e8e7623bc4580fff21161b6e 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -204,6 +204,7 @@ add_executable (mpi_unit_tests
test_ParallelChecker_read.cpp
test_Partitioner.cpp
test_RandomEngine.cpp
+ test_StencilBuilder.cpp
test_SubItemArrayPerItemVariant.cpp
test_SubItemValuePerItem.cpp
test_SubItemValuePerItemVariant.cpp
diff --git a/tests/test_StencilBuilder.cpp b/tests/test_StencilBuilder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..af6031ee499253c283683e7e15d7e779327def17
--- /dev/null
+++ b/tests/test_StencilBuilder.cpp
@@ -0,0 +1,113 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <mesh/StencilBuilder.hpp>
+#include <utils/Messenger.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("StencilBuilder", "[mesh]")
+{
+ auto is_valid = [](const auto& connectivity, const Stencil& stencil) {
+ auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ auto cell_number = connectivity.cellNumber();
+
+ for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+ if (cell_is_owned[cell_id]) {
+ std::set<CellId, std::function<bool(CellId, CellId)>> cell_set(
+ [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; });
+ auto cell_nodes = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ const NodeId node_id = cell_nodes[i_node];
+ auto node_cells = node_to_cell_matrix[node_id];
+ for (size_t i_node_cell = 0; i_node_cell < node_cells.size(); ++i_node_cell) {
+ const CellId node_cell_id = node_cells[i_node_cell];
+ if (node_cell_id != cell_id) {
+ cell_set.insert(node_cell_id);
+ }
+ }
+ }
+
+ auto cell_stencil = stencil[cell_id];
+
+ auto i_set_cell = cell_set.begin();
+ for (size_t index = 0; index < cell_stencil.size(); ++index, ++i_set_cell) {
+ if (*i_set_cell != cell_stencil[index]) {
+ return false;
+ }
+ }
+ }
+ }
+ return true;
+ };
+
+ SECTION("1D")
+ {
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+
+ Stencil stencil = StencilBuilder{}.build(connectivity);
+
+ REQUIRE(is_valid(connectivity, stencil));
+ }
+
+ SECTION("unordered")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>();
+
+ const Connectivity<1>& connectivity = mesh.connectivity();
+ REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
+ }
+ }
+
+ SECTION("2D")
+ {
+ SECTION("cartesian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
+
+ const Connectivity<2>& connectivity = mesh.connectivity();
+ REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
+ }
+ }
+
+ SECTION("3D")
+ {
+ SECTION("carteian")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
+ }
+
+ SECTION("hybrid")
+ {
+ const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
+
+ const Connectivity<3>& connectivity = mesh.connectivity();
+ REQUIRE(is_valid(connectivity, StencilBuilder{}.build(connectivity)));
+ }
+ }
+}