From 3f96ea84602615e8aed28b3b8ae5531279cbdc59 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 18 Feb 2022 16:07:12 +0100
Subject: [PATCH] Add synchronize, isSynchronized, min, max and sum utilities
for SubItemArrayPerItem
---
src/mesh/SubItemArrayPerItem.hpp | 2 +-
src/mesh/SubItemArrayPerItemUtils.hpp | 359 ++++++++++++++++++++
tests/CMakeLists.txt | 1 +
tests/test_SubItemArrayPerItemUtils.cpp | 430 ++++++++++++++++++++++++
4 files changed, 791 insertions(+), 1 deletion(-)
create mode 100644 src/mesh/SubItemArrayPerItemUtils.hpp
create mode 100644 tests/test_SubItemArrayPerItemUtils.cpp
diff --git a/src/mesh/SubItemArrayPerItem.hpp b/src/mesh/SubItemArrayPerItem.hpp
index f2103789e..7b204d219 100644
--- a/src/mesh/SubItemArrayPerItem.hpp
+++ b/src/mesh/SubItemArrayPerItem.hpp
@@ -50,7 +50,7 @@ class SubItemArrayPerItem
public:
friend PUGS_INLINE SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr>
- copy(SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
+ copy(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& source)
{
SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem, ConnectivityPtr> image;
diff --git a/src/mesh/SubItemArrayPerItemUtils.hpp b/src/mesh/SubItemArrayPerItemUtils.hpp
new file mode 100644
index 000000000..5e4f122e3
--- /dev/null
+++ b/src/mesh/SubItemArrayPerItemUtils.hpp
@@ -0,0 +1,359 @@
+#ifndef SUB_ITEM_ARRAY_PER_ITEM_UTILS_HPP
+#define SUB_ITEM_ARRAY_PER_ITEM_UTILS_HPP
+
+#include <utils/Messenger.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/SubItemArrayPerItem.hpp>
+#include <mesh/Synchronizer.hpp>
+#include <mesh/SynchronizerManager.hpp>
+#include <utils/PugsTraits.hpp>
+
+#include <iostream>
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+std::remove_const_t<DataType>
+min(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_array_per_item)
+{
+ using SubItemArrayPerItemType = SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>;
+ constexpr ItemType item_type = ItemOfItem::item_type;
+ using ItemIsOwnedType = ItemValue<const bool, item_type>;
+ using data_type = std::remove_const_t<typename SubItemArrayPerItemType::data_type>;
+ using index_type = typename SubItemArrayPerItemType::index_type;
+
+ static_assert(std::is_arithmetic_v<data_type>, "min cannot be called on non-arithmetic data");
+ static_assert(not std::is_same_v<data_type, bool>, "min cannot be called on boolean data");
+
+ class SubItemArrayPerItemMin
+ {
+ private:
+ const SubItemArrayPerItemType& m_sub_item_array_per_item;
+ const ItemIsOwnedType m_is_owned;
+
+ public:
+ PUGS_INLINE
+ operator data_type()
+ {
+ data_type reduced_array;
+ parallel_reduce(m_sub_item_array_per_item.numberOfItems(), *this, reduced_array);
+ return reduced_array;
+ }
+
+ PUGS_INLINE
+ void
+ operator()(const index_type& item_id, data_type& data) const
+ {
+ if (m_is_owned[item_id]) {
+ Table sub_item_table = m_sub_item_array_per_item.itemTable(item_id);
+ for (size_t i = 0; i < sub_item_table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < sub_item_table.numberOfColumns(); ++j) {
+ if (sub_item_table(i, j) < data) {
+ data = sub_item_table(i, j);
+ }
+ }
+ }
+ }
+ }
+
+ PUGS_INLINE
+ void
+ join(volatile data_type& dst, const volatile data_type& src) const
+ {
+ if (src < dst) {
+ // cannot be reached if initial array is the min
+ dst = src; // LCOV_EXCL_LINE
+ }
+ }
+
+ PUGS_INLINE
+ void
+ init(data_type& array) const
+ {
+ array = std::numeric_limits<data_type>::max();
+ }
+
+ PUGS_INLINE
+ SubItemArrayPerItemMin(const SubItemArrayPerItemType& item_array)
+ : m_sub_item_array_per_item(item_array), m_is_owned([&](const IConnectivity& connectivity) {
+ Assert((connectivity.dimension() > 0) and (connectivity.dimension() <= 3),
+ "unexpected connectivity dimension");
+
+ switch (connectivity.dimension()) {
+ case 1: {
+ const auto& connectivity_1d = static_cast<const Connectivity1D&>(connectivity);
+ return connectivity_1d.isOwned<item_type>();
+ break;
+ }
+ case 2: {
+ const auto& connectivity_2d = static_cast<const Connectivity2D&>(connectivity);
+ return connectivity_2d.isOwned<item_type>();
+ break;
+ }
+ case 3: {
+ const auto& connectivity_3d = static_cast<const Connectivity3D&>(connectivity);
+ return connectivity_3d.isOwned<item_type>();
+ break;
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("unexpected dimension");
+ }
+ // LCOV_EXCL_STOP
+ }
+ }(*item_array.connectivity_ptr()))
+ {
+ ;
+ }
+
+ PUGS_INLINE
+ ~SubItemArrayPerItemMin() = default;
+ };
+
+ const DataType local_min = SubItemArrayPerItemMin{sub_item_array_per_item};
+ return parallel::allReduceMin(local_min);
+}
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+std::remove_const_t<DataType>
+max(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_array_per_item)
+{
+ using SubItemArrayPerItemType = SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>;
+ constexpr ItemType item_type = ItemOfItem::item_type;
+ using ItemIsOwnedType = ItemValue<const bool, item_type>;
+ using data_type = std::remove_const_t<typename SubItemArrayPerItemType::data_type>;
+ using index_type = typename SubItemArrayPerItemType::index_type;
+
+ static_assert(std::is_arithmetic_v<data_type>, "min cannot be called on non-arithmetic data");
+ static_assert(not std::is_same_v<data_type, bool>, "max cannot be called on boolean data");
+
+ class SubItemArrayPerItemMax
+ {
+ private:
+ const SubItemArrayPerItemType& m_sub_item_array_per_item;
+ const ItemIsOwnedType m_is_owned;
+
+ public:
+ PUGS_INLINE
+ operator data_type()
+ {
+ data_type reduced_array;
+ parallel_reduce(m_sub_item_array_per_item.numberOfItems(), *this, reduced_array);
+ return reduced_array;
+ }
+
+ PUGS_INLINE
+ void
+ operator()(const index_type& item_id, data_type& data) const
+ {
+ if (m_is_owned[item_id]) {
+ Table sub_item_table = m_sub_item_array_per_item.itemTable(item_id);
+ for (size_t i = 0; i < sub_item_table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < sub_item_table.numberOfColumns(); ++j) {
+ if (sub_item_table(i, j) > data) {
+ data = sub_item_table(i, j);
+ }
+ }
+ }
+ }
+ }
+
+ PUGS_INLINE
+ void
+ join(volatile data_type& dst, const volatile data_type& src) const
+ {
+ if (src > dst) {
+ // cannot be reached if initial array is the max
+ dst = src; // LCOV_EXCL_LINE
+ }
+ }
+
+ PUGS_INLINE
+ void
+ init(data_type& array) const
+ {
+ array = std::numeric_limits<data_type>::min();
+ }
+
+ PUGS_INLINE
+ SubItemArrayPerItemMax(const SubItemArrayPerItemType& item_array)
+ : m_sub_item_array_per_item(item_array), m_is_owned([&](const IConnectivity& connectivity) {
+ Assert((connectivity.dimension() > 0) and (connectivity.dimension() <= 3),
+ "unexpected connectivity dimension");
+
+ switch (connectivity.dimension()) {
+ case 1: {
+ const auto& connectivity_1d = static_cast<const Connectivity1D&>(connectivity);
+ return connectivity_1d.isOwned<item_type>();
+ break;
+ }
+ case 2: {
+ const auto& connectivity_2d = static_cast<const Connectivity2D&>(connectivity);
+ return connectivity_2d.isOwned<item_type>();
+ break;
+ }
+ case 3: {
+ const auto& connectivity_3d = static_cast<const Connectivity3D&>(connectivity);
+ return connectivity_3d.isOwned<item_type>();
+ break;
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("unexpected dimension");
+ }
+ // LCOV_EXCL_STOP
+ }
+ }(*item_array.connectivity_ptr()))
+ {
+ ;
+ }
+
+ PUGS_INLINE
+ ~SubItemArrayPerItemMax() = default;
+ };
+
+ const DataType local_max = SubItemArrayPerItemMax{sub_item_array_per_item};
+ return parallel::allReduceMax(local_max);
+}
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+std::remove_const_t<DataType>
+sum(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& item_array)
+{
+ using SubItemArrayPerItemType = SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>;
+ constexpr ItemType item_type = ItemOfItem::item_type;
+ using ItemIsOwnedType = ItemValue<const bool, item_type>;
+ using data_type = std::remove_const_t<typename SubItemArrayPerItemType::data_type>;
+ using index_type = typename SubItemArrayPerItemType::index_type;
+
+ static_assert(not std::is_same_v<data_type, bool>, "sum cannot be called on boolean data");
+
+ class SubItemArrayPerItemSum
+ {
+ private:
+ const SubItemArrayPerItemType& m_sub_item_array_per_item;
+ const ItemIsOwnedType m_is_owned;
+
+ public:
+ PUGS_INLINE
+ operator data_type()
+ {
+ data_type reduced_array;
+ parallel_reduce(m_sub_item_array_per_item.numberOfItems(), *this, reduced_array);
+ return reduced_array;
+ }
+
+ PUGS_INLINE
+ void
+ operator()(const index_type& item_id, data_type& data) const
+ {
+ if (m_is_owned[item_id]) {
+ Table sub_item_table = m_sub_item_array_per_item.itemTable(item_id);
+ for (size_t i = 0; i < sub_item_table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < sub_item_table.numberOfColumns(); ++j) {
+ data += sub_item_table(i, j);
+ }
+ }
+ }
+ }
+
+ PUGS_INLINE
+ void
+ join(volatile data_type& dst, const volatile data_type& src) const
+ {
+ dst += src;
+ }
+
+ PUGS_INLINE
+ void
+ init(data_type& array) const
+ {
+ if constexpr (std::is_arithmetic_v<data_type>) {
+ array = 0;
+ } else {
+ static_assert(is_tiny_vector_v<data_type> or is_tiny_matrix_v<data_type>, "invalid data type");
+ array = zero;
+ }
+ }
+
+ PUGS_INLINE
+ SubItemArrayPerItemSum(const SubItemArrayPerItemType& item_array)
+ : m_sub_item_array_per_item(item_array), m_is_owned([&](const IConnectivity& connectivity) {
+ Assert((connectivity.dimension() > 0) and (connectivity.dimension() <= 3),
+ "unexpected connectivity dimension");
+
+ switch (connectivity.dimension()) {
+ case 1: {
+ const auto& connectivity_1d = static_cast<const Connectivity1D&>(connectivity);
+ return connectivity_1d.isOwned<item_type>();
+ break;
+ }
+ case 2: {
+ const auto& connectivity_2d = static_cast<const Connectivity2D&>(connectivity);
+ return connectivity_2d.isOwned<item_type>();
+ break;
+ }
+ case 3: {
+ const auto& connectivity_3d = static_cast<const Connectivity3D&>(connectivity);
+ return connectivity_3d.isOwned<item_type>();
+ break;
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("unexpected dimension");
+ }
+ // LCOV_EXCL_STOP
+ }
+ }(*item_array.connectivity_ptr()))
+ {
+ ;
+ }
+
+ PUGS_INLINE
+ ~SubItemArrayPerItemSum() = default;
+ };
+
+ const DataType local_sum = SubItemArrayPerItemSum{item_array};
+ return parallel::allReduceSum(local_sum);
+}
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+void
+synchronize(SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_array_per_item)
+{
+ static_assert(not std::is_const_v<DataType>, "cannot synchronize SubItemArrayPerItem of const data");
+ if (parallel::size() > 1) {
+ auto& manager = SynchronizerManager::instance();
+ const IConnectivity* connectivity = sub_item_array_per_item.connectivity_ptr().get();
+ Synchronizer& synchronizer = manager.getConnectivitySynchronizer(connectivity);
+ synchronizer.synchronize(sub_item_array_per_item);
+ }
+}
+
+template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
+bool
+isSynchronized(const SubItemArrayPerItem<DataType, ItemOfItem, ConnectivityPtr>& sub_item_array_per_item)
+{
+ bool is_synchronized = true;
+ if (parallel::size() > 1) {
+ SubItemArrayPerItem<std::remove_const_t<DataType>, ItemOfItem> sub_item_array_per_item_copy =
+ copy(sub_item_array_per_item);
+
+ synchronize(sub_item_array_per_item_copy);
+
+ for (size_t i = 0; i < sub_item_array_per_item.numberOfArrays(); ++i) {
+ Array sub_item_array = sub_item_array_per_item[i];
+ Array sub_item_array_copy = sub_item_array_per_item_copy[i];
+ for (size_t j = 0; j < sub_item_array.size(); ++j) {
+ if (sub_item_array_copy[j] != sub_item_array[j]) {
+ is_synchronized = false;
+ }
+ }
+ }
+
+ is_synchronized = parallel::allReduceAnd(is_synchronized);
+ }
+
+ return is_synchronized;
+}
+
+#endif // SUB_ITEM_ARRAY_PER_ITEM_UTILS_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 594384931..cbd6b32bc 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -175,6 +175,7 @@ add_executable (mpi_unit_tests
test_SubItemValuePerItem.cpp
test_SubItemValuePerItemUtils.cpp
test_SubItemArrayPerItem.cpp
+ test_SubItemArrayPerItemUtils.cpp
test_Synchronizer.cpp
)
diff --git a/tests/test_SubItemArrayPerItemUtils.cpp b/tests/test_SubItemArrayPerItemUtils.cpp
new file mode 100644
index 000000000..d3f29a381
--- /dev/null
+++ b/tests/test_SubItemArrayPerItemUtils.cpp
@@ -0,0 +1,430 @@
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_all.hpp>
+
+#include <MeshDataBaseForTests.hpp>
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/SubItemArrayPerItem.hpp>
+#include <mesh/SubItemArrayPerItemUtils.hpp>
+#include <utils/Messenger.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SubItemArrayPerItemUtils", "[mesh]")
+{
+ SECTION("Synchronize")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_2d = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+ NodeArrayPerFace<int> weak_node_array_per_face{connectivity, 3};
+
+ for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+ auto face_table = weak_node_array_per_face.itemTable(face_id);
+ for (size_t i_node = 0; i_node < face_table.numberOfRows(); ++i_node) {
+ for (size_t l = 0; l < face_table.numberOfColumns(); ++l) {
+ face_table(i_node, l) = parallel::rank() + 2 * i_node + l;
+ }
+ }
+ }
+
+ NodeArrayPerFace<const int> node_array_per_face{weak_node_array_per_face};
+
+ REQUIRE(node_array_per_face.connectivity_ptr() == weak_node_array_per_face.connectivity_ptr());
+
+ if (parallel::size() > 1) {
+ // before synchronization
+ auto face_owner = connectivity.faceOwner();
+
+ bool is_synchronized = true;
+ for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+ auto face_table = node_array_per_face.itemTable(face_id);
+ for (size_t i_node = 0; i_node < face_table.numberOfRows(); ++i_node) {
+ for (size_t l = 0; l < face_table.numberOfColumns(); ++l) {
+ if (face_table(i_node, l) != static_cast<int>(face_owner[face_id] + 2 * i_node + l)) {
+ is_synchronized = false;
+ break;
+ }
+ }
+ }
+ }
+
+ REQUIRE(not is_synchronized);
+ REQUIRE(not isSynchronized(node_array_per_face));
+ }
+
+ synchronize(weak_node_array_per_face);
+
+ { // after synchronization
+ auto face_owner = connectivity.faceOwner();
+
+ bool is_synchronized = true;
+ for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+ auto face_table = node_array_per_face.itemTable(face_id);
+ for (size_t i_node = 0; i_node < face_table.numberOfRows(); ++i_node) {
+ for (size_t l = 0; l < face_table.numberOfColumns(); ++l) {
+ if (face_table(i_node, l) != static_cast<int>(face_owner[face_id] + 2 * i_node + l)) {
+ is_synchronized = false;
+ }
+ }
+ }
+ }
+ REQUIRE(is_synchronized);
+ }
+ REQUIRE(isSynchronized(node_array_per_face));
+ }
+ }
+ }
+
+ SECTION("min")
+ {
+ SECTION("1D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_1d = named_mesh.mesh();
+
+ const Connectivity<1>& connectivity = mesh_1d->connectivity();
+
+ NodeArrayPerCell<int> node_array_per_cell{connectivity, 3};
+ node_array_per_cell.fill(-1);
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ parallel_for(
+ connectivity.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (cell_is_owned[cell_id]) {
+ auto cell_table = node_array_per_cell.itemTable(cell_id);
+ for (size_t i = 0; i < cell_table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < cell_table.numberOfColumns(); ++j) {
+ cell_table(i, j) = 10 + parallel::rank() + i + 2 * j;
+ }
+ }
+ }
+ });
+
+ REQUIRE(min(node_array_per_cell) == 10);
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_2d = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+ FaceArrayPerCell<int> face_array_per_cell{connectivity, 3};
+ face_array_per_cell.fill(-1);
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ parallel_for(
+ connectivity.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (cell_is_owned[cell_id]) {
+ auto cell_table = face_array_per_cell.itemTable(cell_id);
+ for (size_t i = 0; i < cell_table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < cell_table.numberOfColumns(); ++j) {
+ cell_table(i, j) = 10 + parallel::rank() + i + j;
+ }
+ }
+ }
+ });
+
+ REQUIRE(min(face_array_per_cell) == 10);
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_3d = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+ EdgeArrayPerFace<int> edge_array_per_face{connectivity, 3};
+ edge_array_per_face.fill(-1);
+
+ auto face_is_owned = connectivity.faceIsOwned();
+ parallel_for(
+ connectivity.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ if (face_is_owned[face_id]) {
+ auto face_table = edge_array_per_face.itemTable(face_id);
+ for (size_t i = 0; i < face_table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < face_table.numberOfColumns(); ++j) {
+ face_table(i, j) = 10 + parallel::rank() + i + j;
+ }
+ }
+ }
+ });
+
+ REQUIRE(min(edge_array_per_face) == 10);
+ }
+ }
+ }
+ }
+
+ SECTION("max")
+ {
+ SECTION("1D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_1d = named_mesh.mesh();
+
+ const Connectivity<1>& connectivity = mesh_1d->connectivity();
+
+ EdgeArrayPerCell<size_t> edge_array_per_cell{connectivity, 3};
+ edge_array_per_cell.fill(std::numeric_limits<size_t>::max());
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ parallel_for(
+ connectivity.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (cell_is_owned[cell_id]) {
+ auto cell_table = edge_array_per_cell.itemTable(cell_id);
+ for (size_t i = 0; i < cell_table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < cell_table.numberOfColumns(); ++j) {
+ cell_table(i, j) = 10 + parallel::rank() - i - j;
+ }
+ }
+ }
+ });
+
+ REQUIRE(max(edge_array_per_cell) == 9 + parallel::size());
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_2d = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+ EdgeArrayPerCell<size_t> edge_array_per_cell{connectivity, 3};
+ edge_array_per_cell.fill(std::numeric_limits<size_t>::max());
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+ parallel_for(
+ connectivity.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (cell_is_owned[cell_id]) {
+ auto cell_table = edge_array_per_cell.itemTable(cell_id);
+ for (size_t i = 0; i < cell_table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < cell_table.numberOfColumns(); ++j) {
+ cell_table(i, j) = 10 + parallel::rank() - i - j;
+ }
+ }
+ }
+ });
+
+ REQUIRE(max(edge_array_per_cell) == 9 + parallel::size());
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_3d = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+ NodeArrayPerEdge<size_t> node_array_per_edge{connectivity, 3};
+ node_array_per_edge.fill(std::numeric_limits<size_t>::max());
+
+ auto edge_is_owned = connectivity.edgeIsOwned();
+
+ parallel_for(
+ connectivity.numberOfEdges(), PUGS_LAMBDA(EdgeId edge_id) {
+ if (edge_is_owned[edge_id]) {
+ auto edge_table = node_array_per_edge.itemTable(edge_id);
+ for (size_t i = 0; i < edge_table.numberOfRows(); ++i) {
+ for (size_t j = 0; j < edge_table.numberOfColumns(); ++j) {
+ edge_table(i, j) = 10 + parallel::rank() - i - j;
+ }
+ }
+ }
+ });
+
+ REQUIRE(max(node_array_per_edge) == 9 + parallel::size());
+ }
+ }
+ }
+ }
+
+ SECTION("sum")
+ {
+ SECTION("1D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ auto mesh_1d = named_mesh.mesh();
+
+ const Connectivity<1>& connectivity = mesh_1d->connectivity();
+
+ NodeArrayPerCell<size_t> node_array_per_cell{connectivity, 3};
+ node_array_per_cell.fill(5);
+
+ auto cell_is_owned = connectivity.cellIsOwned();
+
+ const size_t global_number_of_nodes_per_cells = [&] {
+ size_t number_of_nodes_per_cells = 0;
+ for (CellId cell_id = 0; cell_id < cell_is_owned.numberOfItems(); ++cell_id) {
+ number_of_nodes_per_cells += cell_is_owned[cell_id] * node_array_per_cell.numberOfSubArrays(cell_id);
+ }
+ return parallel::allReduceSum(number_of_nodes_per_cells);
+ }();
+
+ REQUIRE(sum(node_array_per_cell) ==
+ 5 * global_number_of_nodes_per_cells * node_array_per_cell.sizeOfArrays());
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name() + " for size_t data")
+ {
+ auto mesh_2d = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+ NodeArrayPerFace<size_t> node_array_per_face{connectivity, 3};
+ node_array_per_face.fill(2);
+
+ auto face_is_owned = connectivity.faceIsOwned();
+
+ const size_t global_number_of_nodes_per_faces = [&] {
+ size_t number_of_nodes_per_faces = 0;
+ for (FaceId face_id = 0; face_id < face_is_owned.numberOfItems(); ++face_id) {
+ number_of_nodes_per_faces += face_is_owned[face_id] * node_array_per_face.numberOfSubArrays(face_id);
+ }
+ return parallel::allReduceSum(number_of_nodes_per_faces);
+ }();
+
+ REQUIRE(sum(node_array_per_face) ==
+ 2 * global_number_of_nodes_per_faces * node_array_per_face.sizeOfArrays());
+ }
+
+ SECTION(named_mesh.name() + " for N^2 data")
+ {
+ auto mesh_2d = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh_2d->connectivity();
+
+ using N2 = TinyVector<2, size_t>;
+ NodeArrayPerFace<N2> node_array_per_face{connectivity, 3};
+
+ const N2 data(3, 2);
+ node_array_per_face.fill(data);
+
+ auto face_is_owned = connectivity.faceIsOwned();
+
+ const size_t global_number_of_nodes_per_faces = [&] {
+ size_t number_of_nodes_per_faces = 0;
+ for (FaceId face_id = 0; face_id < face_is_owned.numberOfItems(); ++face_id) {
+ number_of_nodes_per_faces += face_is_owned[face_id] * node_array_per_face.numberOfSubArrays(face_id);
+ }
+ return parallel::allReduceSum(number_of_nodes_per_faces);
+ }();
+
+ REQUIRE(sum(node_array_per_face) ==
+ global_number_of_nodes_per_faces * node_array_per_face.sizeOfArrays() * data);
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+ for (auto named_mesh : mesh_list) {
+ SECTION(named_mesh.name() + " for size_t data")
+ {
+ auto mesh_3d = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+ EdgeArrayPerFace<size_t> edge_array_per_face{connectivity, 3};
+ edge_array_per_face.fill(3);
+
+ auto face_is_owned = connectivity.faceIsOwned();
+
+ const size_t global_number_of_edges_per_faces = [&] {
+ size_t number_of_edges_per_faces = 0;
+ for (FaceId face_id = 0; face_id < face_is_owned.numberOfItems(); ++face_id) {
+ number_of_edges_per_faces += face_is_owned[face_id] * edge_array_per_face.numberOfSubArrays(face_id);
+ }
+ return parallel::allReduceSum(number_of_edges_per_faces);
+ }();
+
+ REQUIRE(sum(edge_array_per_face) ==
+ 3 * global_number_of_edges_per_faces * edge_array_per_face.sizeOfArrays());
+ }
+
+ SECTION(named_mesh.name() + " for N^2x3 data")
+ {
+ auto mesh_3d = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh_3d->connectivity();
+
+ using N2x3 = TinyMatrix<2, 3, size_t>;
+ EdgeArrayPerFace<N2x3> edge_array_per_face{connectivity, 3};
+
+ const N2x3 data(1, 3, 4, 6, 2, 5);
+ edge_array_per_face.fill(data);
+
+ auto face_is_owned = connectivity.faceIsOwned();
+
+ const size_t global_number_of_edges_per_faces = [&] {
+ size_t number_of_edges_per_faces = 0;
+ for (FaceId face_id = 0; face_id < face_is_owned.numberOfItems(); ++face_id) {
+ number_of_edges_per_faces += face_is_owned[face_id] * edge_array_per_face.numberOfSubArrays(face_id);
+ }
+ return parallel::allReduceSum(number_of_edges_per_faces);
+ }();
+
+ REQUIRE(sum(edge_array_per_face) ==
+ global_number_of_edges_per_faces * edge_array_per_face.sizeOfArrays() * data);
+ }
+ }
+ }
+ }
+}
--
GitLab