diff --git a/src/mesh/StencilArray.hpp b/src/mesh/StencilArray.hpp index 9571c346bae6e4e345a16b09188dce15ab73b456..ffe8638bf0d8fdb472780150603c7fab21ca1eb9 100644 --- a/src/mesh/StencilArray.hpp +++ b/src/mesh/StencilArray.hpp @@ -12,6 +12,9 @@ class StencilArray public: using ItemToItemMatrixT = ItemToItemMatrix<source_item_type, target_item_type>; + using SourceItemId = ItemIdT<source_item_type>; + using TargetItemId = ItemIdT<target_item_type>; + class BoundaryDescriptorStencilArray { private: @@ -71,9 +74,9 @@ class StencilArray PUGS_INLINE auto - operator[](CellId cell_id) const + operator[](SourceItemId source_item_id) const { - return m_stencil_array[cell_id]; + return m_stencil_array[source_item_id]; } StencilArray(const ConnectivityMatrix& connectivity_matrix, diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp index 96f36a14f99ce3a3239f94d1b454f35ad82fc22b..20fc3bf61dd6a59ef82ff7bdb25dc5e32c7eef24 100644 --- a/src/mesh/StencilBuilder.cpp +++ b/src/mesh/StencilBuilder.cpp @@ -478,8 +478,291 @@ StencilBuilder::_build_for_different_source_and_target( const size_t& number_of_layers, const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const { - static_assert(source_item_type != target_item_type); - throw NotImplementedError("different source target"); + constexpr size_t Dimension = ConnectivityType::Dimension; + + if (number_of_layers == 0) { + throw NormalError("number of layers must be greater than 0 to build stencils"); + } + + auto connecting_item_to_target_item_matrix = + connectivity.template getItemToItemMatrix<connecting_item_type, target_item_type>(); + auto target_item_to_connecting_item_matrix = + connectivity.template getItemToItemMatrix<target_item_type, connecting_item_type>(); + + auto source_item_to_connecting_item_matrix = + [](const auto& given_connectivity) -> ItemToItemMatrix<source_item_type, connecting_item_type> { + if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) == + ItemTypeId<Dimension>::itemTId(connecting_item_type)) { + return {ConnectivityMatrix{}}; + } else { + return given_connectivity.template getItemToItemMatrix<source_item_type, connecting_item_type>(); + } + }(connectivity); + + auto source_item_is_owned = connectivity.template isOwned<source_item_type>(); + auto target_item_number = connectivity.template number<target_item_type>(); + auto connecting_item_number = connectivity.template number<connecting_item_type>(); + + using SourceItemId = ItemIdT<source_item_type>; + using TargetItemId = ItemIdT<target_item_type>; + using ConnectingItemId = ItemIdT<connecting_item_type>; + + ItemArray<bool, connecting_item_type> symmetry_item_list = + this->_buildSymmetryConnectingItemList<connecting_item_type>(connectivity, symmetry_boundary_descriptor_list); + + Array<uint32_t> row_map{connectivity.template numberOf<source_item_type>() + 1}; + row_map[0] = 0; + std::vector<Array<uint32_t>> symmetry_row_map_list(symmetry_boundary_descriptor_list.size()); + for (auto&& symmetry_row_map : symmetry_row_map_list) { + symmetry_row_map = Array<uint32_t>{connectivity.template numberOf<source_item_type>() + 1}; + symmetry_row_map[0] = 0; + } + + std::vector<TargetItemId> column_indices_vector; + std::vector<std::vector<uint32_t>> symmetry_column_indices_vector(symmetry_boundary_descriptor_list.size()); + + std::vector<Layer<target_item_type>> target_item_layer_list; + std::vector<Layer<connecting_item_type>> connecting_layer_list; + + std::vector<Layer<target_item_type>> symmetry_target_item_layer_list; + std::vector<Layer<connecting_item_type>> symmetry_connecting_layer_list; + + for (size_t i = 0; i < number_of_layers; ++i) { + target_item_layer_list.emplace_back(Layer<target_item_type>{connectivity}); + connecting_layer_list.emplace_back(Layer<connecting_item_type>{connectivity}); + + if (symmetry_boundary_descriptor_list.size() > 0) { + symmetry_target_item_layer_list.emplace_back(Layer<target_item_type>{connectivity}); + symmetry_connecting_layer_list.emplace_back(Layer<connecting_item_type>{connectivity}); + } + } + + for (SourceItemId source_item_id = 0; source_item_id < connectivity.template numberOf<source_item_type>(); + ++source_item_id) { + if (source_item_is_owned[source_item_id]) { + for (auto&& target_item_layer : target_item_layer_list) { + target_item_layer.clear(); + } + for (auto&& connecting_layer : connecting_layer_list) { + connecting_layer.clear(); + } + + // First layer is a special case + { + auto& target_item_layer = target_item_layer_list[0]; + auto& connecting_layer = connecting_layer_list[0]; + + if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) == + ItemTypeId<Dimension>::itemTId(connecting_item_type)) { + connecting_layer.add(ConnectingItemId{static_cast<typename ConnectingItemId::base_type>(source_item_id)}); + } else { + for (size_t i_connecting_item = 0; + i_connecting_item < source_item_to_connecting_item_matrix[source_item_id].size(); ++i_connecting_item) { + const ConnectingItemId connecting_item_id = + source_item_to_connecting_item_matrix[source_item_id][i_connecting_item]; + + connecting_layer.add(connecting_item_id); + } + } + + for (auto connecting_item_id : connecting_layer.itemIdList()) { + for (size_t i_item = 0; i_item < connecting_item_to_target_item_matrix[connecting_item_id].size(); ++i_item) { + const TargetItemId layer_item_id = connecting_item_to_target_item_matrix[connecting_item_id][i_item]; + target_item_layer.add(layer_item_id); + } + } + } + + for (size_t i_layer = 1; i_layer < number_of_layers; ++i_layer) { + auto& connecting_layer = connecting_layer_list[i_layer]; + + auto has_connecting_item = [&i_layer, &connecting_layer_list](ConnectingItemId connecting_item_id) -> bool { + for (ssize_t i = i_layer - 1; i >= 0; --i) { + if (connecting_layer_list[i].hasItemNumber(connecting_item_id)) { + return true; + } + } + + return false; + }; + + for (auto&& previous_layer_item_id_list : target_item_layer_list[i_layer - 1].itemIdList()) { + const auto target_item_to_connecting_item_list = + target_item_to_connecting_item_matrix[previous_layer_item_id_list]; + for (size_t i_connecting_item = 0; i_connecting_item < target_item_to_connecting_item_list.size(); + ++i_connecting_item) { + const ConnectingItemId connecting_item_id = target_item_to_connecting_item_list[i_connecting_item]; + + if (not has_connecting_item(connecting_item_id)) { + connecting_layer.add(connecting_item_id); + } + } + } + + auto& target_item_layer = target_item_layer_list[i_layer]; + + auto has_layer_item = [&i_layer, &target_item_layer_list](TargetItemId layer_item_id) -> bool { + for (ssize_t i = i_layer - 1; i >= 0; --i) { + if (target_item_layer_list[i].hasItemNumber(layer_item_id)) { + return true; + } + } + + return false; + }; + + for (auto connecting_item_id : connecting_layer.itemIdList()) { + const auto& connecting_item_to_target_item_list = connecting_item_to_target_item_matrix[connecting_item_id]; + for (size_t i_item = 0; i_item < connecting_item_to_target_item_list.size(); ++i_item) { + const TargetItemId layer_item_id = connecting_item_to_target_item_list[i_item]; + if (not has_layer_item(layer_item_id)) { + target_item_layer.add(layer_item_id); + } + } + } + } + + for (size_t i_symmetry = 0; i_symmetry < symmetry_boundary_descriptor_list.size(); ++i_symmetry) { + for (auto&& symmetry_target_item_layer : symmetry_target_item_layer_list) { + symmetry_target_item_layer.clear(); + } + for (auto&& symmetry_connecting_layer : symmetry_connecting_layer_list) { + symmetry_connecting_layer.clear(); + } + + // First layer is a special case + for (auto&& connecting_item_id : connecting_layer_list[0].itemIdList()) { + if (symmetry_item_list[connecting_item_id][i_symmetry]) { + symmetry_connecting_layer_list[0].add(connecting_item_id); + } + } + + for (auto&& connecting_item_id : symmetry_connecting_layer_list[0].itemIdList()) { + auto item_of_connecting_item_list = connecting_item_to_target_item_matrix[connecting_item_id]; + for (size_t i_item_of_connecting_item = 0; i_item_of_connecting_item < item_of_connecting_item_list.size(); + ++i_item_of_connecting_item) { + const TargetItemId item_id_of_connecting_item = item_of_connecting_item_list[i_item_of_connecting_item]; + symmetry_target_item_layer_list[0].add(item_id_of_connecting_item); + } + } + + for (size_t i_layer = 1; i_layer < number_of_layers; ++i_layer) { + auto has_connecting_item = [&i_layer, + &symmetry_connecting_layer_list](ConnectingItemId connecting_item_id) -> bool { + for (ssize_t i = i_layer - 1; i >= 0; --i) { + if (symmetry_connecting_layer_list[i].hasItemNumber(connecting_item_id)) { + return true; + } + } + + return false; + }; + + for (auto&& symmetry_connecting_item_id : connecting_layer_list[i_layer].itemIdList()) { + if (symmetry_item_list[symmetry_connecting_item_id][i_symmetry]) { + if (not has_connecting_item(symmetry_connecting_item_id)) { + symmetry_connecting_layer_list[i_layer].add(symmetry_connecting_item_id); + } + } + } + + for (auto&& previous_layer_target_item_id_list : symmetry_target_item_layer_list[i_layer - 1].itemIdList()) { + const auto item_to_connecting_item_list = + target_item_to_connecting_item_matrix[previous_layer_target_item_id_list]; + for (size_t i_connecting_item = 0; i_connecting_item < item_to_connecting_item_list.size(); + ++i_connecting_item) { + const ConnectingItemId connecting_item_id = item_to_connecting_item_list[i_connecting_item]; + + if (not has_connecting_item(connecting_item_id)) { + symmetry_connecting_layer_list[i_layer].add(connecting_item_id); + } + } + } + + auto has_symmetry_layer_item = [&i_layer, + &symmetry_target_item_layer_list](TargetItemId layer_target_item_id) -> bool { + for (ssize_t i = i_layer - 1; i >= 0; --i) { + if (symmetry_target_item_layer_list[i].hasItemNumber(layer_target_item_id)) { + return true; + } + } + + return false; + }; + + for (auto&& connecting_item_id : symmetry_connecting_layer_list[i_layer].itemIdList()) { + auto item_of_connecting_item_list = connecting_item_to_target_item_matrix[connecting_item_id]; + for (size_t i_target_item_of_connecting_item = 0; + i_target_item_of_connecting_item < item_of_connecting_item_list.size(); + ++i_target_item_of_connecting_item) { + const TargetItemId target_item_id_of_connecting_item = + item_of_connecting_item_list[i_target_item_of_connecting_item]; + if (not has_symmetry_layer_item(target_item_id_of_connecting_item)) { + symmetry_target_item_layer_list[i_layer].add(target_item_id_of_connecting_item); + } + } + } + } + + for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) { + for (auto symmetry_layer_target_item_id : symmetry_target_item_layer_list[i_layer].itemIdList()) { + symmetry_column_indices_vector[i_symmetry].push_back(symmetry_layer_target_item_id); + } + } + + { + size_t stencil_size = 0; + for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) { + auto& target_item_layer = symmetry_target_item_layer_list[i_layer]; + stencil_size += target_item_layer.size(); + } + + symmetry_row_map_list[i_symmetry][source_item_id + 1] = + symmetry_row_map_list[i_symmetry][source_item_id] + stencil_size; + } + } + + { + size_t stencil_size = 0; + for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) { + auto& item_layer = target_item_layer_list[i_layer]; + for (auto stencil_item_id : item_layer.itemIdList()) { + column_indices_vector.push_back(stencil_item_id); + } + stencil_size += item_layer.size(); + } + row_map[source_item_id + 1] = row_map[source_item_id] + stencil_size; + } + + } else { + row_map[source_item_id + 1] = row_map[source_item_id]; + for (size_t i = 0; i < symmetry_row_map_list.size(); ++i) { + symmetry_row_map_list[i][source_item_id + 1] = symmetry_row_map_list[i][source_item_id]; + } + } + } + + Array<uint32_t> column_indices{column_indices_vector.size()}; + parallel_for( + column_indices_vector.size(), PUGS_LAMBDA(size_t i) { column_indices[i] = column_indices_vector[i]; }); + + ConnectivityMatrix primal_stencil{row_map, column_indices}; + + typename StencilArray<source_item_type, target_item_type>::BoundaryDescriptorStencilArrayList + symmetry_boundary_stencil_list; + { + size_t i = 0; + for (auto&& p_boundary_descriptor : symmetry_boundary_descriptor_list) { + symmetry_boundary_stencil_list.emplace_back( + typename StencilArray<source_item_type, target_item_type>:: + BoundaryDescriptorStencilArray{p_boundary_descriptor, + ConnectivityMatrix{symmetry_row_map_list[i], + convert_to_array(symmetry_column_indices_vector[i])}}); + ++i; + } + } + + return {{primal_stencil}, {symmetry_boundary_stencil_list}}; } template <ItemType source_item_type, @@ -580,7 +863,35 @@ StencilBuilder::buildC2C(const IConnectivity& connectivity, } NodeToCellStencilArray -StencilBuilder::buildN2C(const IConnectivity&, const StencilDescriptor&, const BoundaryDescriptorList&) const +StencilBuilder::buildN2C(const IConnectivity& connectivity, + const StencilDescriptor& stencil_descriptor, + const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const { - throw NotImplementedError("node to cell stencil"); + if ((parallel::size() > 1) and + (stencil_descriptor.numberOfLayers() > GlobalVariableManager::instance().getNumberOfGhostLayers())) { + std::ostringstream error_msg; + error_msg << "Stencil builder requires" << rang::fgB::yellow << stencil_descriptor.numberOfLayers() + << rang::fg::reset << " layers while parallel number of ghost layer is " + << GlobalVariableManager::instance().getNumberOfGhostLayers() << ".\n"; + error_msg << "Increase the number of ghost layers (using the '--number-of-ghost-layers' option)."; + throw NormalError(error_msg.str()); + } + + switch (connectivity.dimension()) { + case 1: { + return this->_build<ItemType::node, ItemType::cell>(dynamic_cast<const Connectivity<1>&>(connectivity), + stencil_descriptor, symmetry_boundary_descriptor_list); + } + case 2: { + return this->_build<ItemType::node, ItemType::cell>(dynamic_cast<const Connectivity<2>&>(connectivity), + stencil_descriptor, symmetry_boundary_descriptor_list); + } + case 3: { + return this->_build<ItemType::node, ItemType::cell>(dynamic_cast<const Connectivity<3>&>(connectivity), + stencil_descriptor, symmetry_boundary_descriptor_list); + } + default: { + throw UnexpectedError("invalid connectivity dimension"); + } + } } diff --git a/src/scheme/DiscreteFunctionDPkVector.hpp b/src/scheme/DiscreteFunctionDPkVector.hpp index 2fefc078b192304ab0a80b63782986867e13c3d8..e8ac729f4b3e002f8c1c3b2c1cdcbd5f4aef3f77 100644 --- a/src/scheme/DiscreteFunctionDPkVector.hpp +++ b/src/scheme/DiscreteFunctionDPkVector.hpp @@ -139,6 +139,16 @@ class DiscreteFunctionDPkVector return m_by_cell_coefficients[cell_id]; } + PUGS_INLINE auto + componentCoefficients(const CellId cell_id, size_t i_component) const noexcept(NO_ASSERT) + { + Assert(m_mesh_v.use_count() > 0, "DiscreteFunctionDPkVector is not built"); + Assert(i_component < m_number_of_components, "incorrect component number"); + + return ComponentCoefficientView{&m_by_cell_coefficients[cell_id][i_component * m_nb_coefficients_per_component], + m_nb_coefficients_per_component}; + } + PUGS_FORCEINLINE BasisView operator()(const CellId cell_id, size_t i_component) const noexcept(NO_ASSERT) { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 874df648b6b6719d8f5b1e79cfcd6a46683e16ea..c5f872d4e9140a823f7a13ec0525be12101e7521 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -264,7 +264,8 @@ add_executable (mpi_unit_tests test_PolynomialReconstruction.cpp test_PolynomialReconstructionDescriptor.cpp test_RandomEngine.cpp - test_StencilBuilder.cpp + test_StencilBuilder_cell2cell.cpp + test_StencilBuilder_node2cell.cpp test_SubItemArrayPerItemVariant.cpp test_SubItemValuePerItem.cpp test_SubItemValuePerItemVariant.cpp diff --git a/tests/test_StencilBuilder.cpp b/tests/test_StencilBuilder_cell2cell.cpp similarity index 99% rename from tests/test_StencilBuilder.cpp rename to tests/test_StencilBuilder_cell2cell.cpp index 0666ca112873073d243fbc6318bf8aec5573e629..42ceb6a5da583a9869035db09f5239d008c6f780 100644 --- a/tests/test_StencilBuilder.cpp +++ b/tests/test_StencilBuilder_cell2cell.cpp @@ -18,7 +18,7 @@ // clazy:excludeall=non-pod-global-static -TEST_CASE("StencilBuilder", "[mesh]") +TEST_CASE("StencilBuilder cell2cell", "[mesh]") { auto is_valid = []<ItemType connecting_item_type>(const auto& connectivity, const auto& stencil_array, const size_t number_of_layers) { diff --git a/tests/test_StencilBuilder_node2cell.cpp b/tests/test_StencilBuilder_node2cell.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6d9c429c6d4ba903185c90a1e6c37d63bf9c9feb --- /dev/null +++ b/tests/test_StencilBuilder_node2cell.cpp @@ -0,0 +1,1184 @@ +#include <catch2/catch_test_macros.hpp> +#include <catch2/matchers/catch_matchers_all.hpp> + +#include <MeshDataBaseForTests.hpp> +#include <mesh/CartesianMeshBuilder.hpp> +#include <mesh/Connectivity.hpp> +#include <mesh/ConnectivityUtils.hpp> +#include <mesh/ItemValue.hpp> +#include <mesh/ItemValueUtils.hpp> +#include <mesh/Mesh.hpp> +#include <mesh/MeshFaceBoundary.hpp> +#include <mesh/MeshVariant.hpp> +#include <mesh/NamedBoundaryDescriptor.hpp> +#include <mesh/StencilManager.hpp> +#include <utils/Messenger.hpp> + +#include <NbGhostLayersTester.hpp> + +// clazy:excludeall=non-pod-global-static + +TEST_CASE("StencilBuilder node2cell", "[mesh]") +{ + auto is_valid = []<ItemType source_item_type, ItemType connecting_item_type>(const auto& connectivity, + const auto& stencil_array, + const size_t number_of_layers) { + constexpr size_t Dimension = std::decay_t<decltype(connectivity)>::Dimension; + + auto cell_to_connecting_item_matrix = + connectivity.template getItemToItemMatrix<ItemType::cell, connecting_item_type>(); + + auto source_item_to_connecting_item_matrix = + [](const auto& given_connectivity) -> ItemToItemMatrix<source_item_type, connecting_item_type> { + if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) == + ItemTypeId<Dimension>::itemTId(connecting_item_type)) { + return ConnectivityMatrix{}; + } else { + return given_connectivity.template getItemToItemMatrix<source_item_type, connecting_item_type>(); + } + }(connectivity); + + auto connecting_to_cell_matrix = connectivity.template getItemToItemMatrix<connecting_item_type, ItemType::cell>(); + auto cell_is_owned = connectivity.cellIsOwned(); + auto cell_number = connectivity.cellNumber(); + auto source_item_is_owned = connectivity.template isOwned<source_item_type>(); + + using SourceItemId = ItemIdT<source_item_type>; + using ConnectingItemId = ItemIdT<connecting_item_type>; + + for (SourceItemId source_item_id = 0; source_item_id < connectivity.template numberOf<source_item_type>(); + ++source_item_id) { + if (source_item_is_owned[source_item_id]) { + std::vector<CellId> expected_stencil; + + std::set<CellId> marked_cell_set; + if constexpr (source_item_type == ItemType::cell) { + marked_cell_set.insert(source_item_id); + } + + std::set<ConnectingItemId> marked_connecting_item_set; + std::set<ConnectingItemId> layer_connecting_item_set; + { + if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) == + ItemTypeId<Dimension>::itemTId(connecting_item_type)) { + ConnectingItemId connecting_id = + ConnectingItemId{static_cast<typename SourceItemId::base_type>(source_item_id)}; + layer_connecting_item_set.insert(connecting_id); + marked_connecting_item_set.insert(connecting_id); + } else { + auto source_item_to_connecting_item_list = source_item_to_connecting_item_matrix[source_item_id]; + for (size_t i_connecting_item = 0; i_connecting_item < source_item_to_connecting_item_list.size(); + ++i_connecting_item) { + const ConnectingItemId connecting_item_id = source_item_to_connecting_item_list[i_connecting_item]; + layer_connecting_item_set.insert(connecting_item_id); + marked_connecting_item_set.insert(connecting_item_id); + } + } + } + + for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) { + 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]; }); + + for (auto&& connecting_item_id : layer_connecting_item_set) { + auto connecting_item_to_cell_list = connecting_to_cell_matrix[connecting_item_id]; + for (size_t i_connecting_item_of_cell = 0; i_connecting_item_of_cell < connecting_item_to_cell_list.size(); + ++i_connecting_item_of_cell) { + const CellId connecting_item_id_of_cell = connecting_item_to_cell_list[i_connecting_item_of_cell]; + if (not marked_cell_set.contains(connecting_item_id_of_cell)) { + cell_set.insert(connecting_item_id_of_cell); + marked_cell_set.insert(connecting_item_id_of_cell); + } + } + } + + layer_connecting_item_set.clear(); + for (auto layer_cell_id : cell_set) { + auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[layer_cell_id]; + for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size(); + ++i_connecting_item) { + const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item]; + if (not marked_connecting_item_set.contains(connecting_item_id)) { + layer_connecting_item_set.insert(connecting_item_id); + marked_connecting_item_set.insert(connecting_item_id); + } + } + } + + for (auto&& set_cell_id : cell_set) { + expected_stencil.push_back(set_cell_id); + } + } + + auto cell_stencil = stencil_array[source_item_id]; + + auto i_set_cell = expected_stencil.begin(); + if (cell_stencil.size() != expected_stencil.size()) { + return false; + } + + 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("inner stencil") + { + SECTION("1 layer") + { + SECTION("1D") + { + SECTION("cartesian") + { + const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>(); + + const Connectivity<1>& connectivity = mesh.connectivity(); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::node>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::edge>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::face>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_faces}), + 1)); + } + + SECTION("unordered") + { + const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>(); + + const Connectivity<1>& connectivity = mesh.connectivity(); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::node>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::edge>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::face>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_faces}), + 1)); + } + } + + SECTION("2D") + { + SECTION("cartesian") + { + const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>(); + + const Connectivity<2>& connectivity = mesh.connectivity(); + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::node>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::edge>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::face>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_faces}), + 1)); + } + + SECTION("hybrid") + { + const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>(); + + const Connectivity<2>& connectivity = mesh.connectivity(); + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::node>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::edge>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::face>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_faces}), + 1)); + } + } + + SECTION("3D") + { + SECTION("carteian") + { + const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>(); + + const Connectivity<3>& connectivity = mesh.connectivity(); + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::node>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::edge>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::face>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_faces}), + 1)); + } + + SECTION("hybrid") + { + const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); + + const Connectivity<3>& connectivity = mesh.connectivity(); + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::node>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::edge>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::face>(connectivity, + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_faces}), + 1)); + } + } + } + + SECTION("2 layers") + { + NbGhostLayersTester nb_ghost_layers_tester(2); + + SECTION("1D") + { + SECTION("cartesian") + { + auto mesh_v = CartesianMeshBuilder(TinyVector<1>{0}, TinyVector<1>{1}, TinyVector<1, uint64_t>(20)).mesh(); + + const auto& mesh = *mesh_v->get<Mesh<1>>(); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::node>(mesh.connectivity(), + StencilManager::instance() + .getNodeToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_nodes}), + 2)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::edge>(mesh.connectivity(), + StencilManager::instance() + .getNodeToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_edges}), + 2)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::face>(mesh.connectivity(), + StencilManager::instance() + .getNodeToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_faces}), + 2)); + } + } + + SECTION("2D") + { + SECTION("cartesian") + { + auto mesh_v = + CartesianMeshBuilder(TinyVector<2>{0, 0}, TinyVector<2>{1, 2}, TinyVector<2, uint64_t>(5, 7)).mesh(); + const auto& mesh = *mesh_v->get<Mesh<2>>(); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::node>(mesh.connectivity(), + StencilManager::instance() + .getNodeToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_nodes}), + 2)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::edge>(mesh.connectivity(), + StencilManager::instance() + .getNodeToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_edges}), + 2)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::face>(mesh.connectivity(), + StencilManager::instance() + .getNodeToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_faces}), + 2)); + } + } + + SECTION("3D") + { + SECTION("carteian") + { + auto mesh_v = + CartesianMeshBuilder(TinyVector<3>{0, 0, 0}, TinyVector<3>{1, 1, 2}, TinyVector<3, uint64_t>(3, 4, 5)) + .mesh(); + const auto& mesh = *mesh_v->get<Mesh<3>>(); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::node>(mesh.connectivity(), + StencilManager::instance() + .getNodeToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_nodes}), + 2)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::edge>(mesh.connectivity(), + StencilManager::instance() + .getNodeToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_edges}), + 2)); + + REQUIRE( + is_valid.template + operator()<ItemType::node, + ItemType::face>(mesh.connectivity(), + StencilManager::instance() + .getNodeToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_faces}), + 2)); + } + } + } + } + + SECTION("Stencil using symmetries") + { + auto check_ghost_nodes_have_empty_stencils = [](const auto& stencil_array, const auto& connecticity) { + bool is_empty = true; + + auto node_is_owned = connecticity.nodeIsOwned(); + + for (NodeId node_id = 0; node_id < connecticity.numberOfNodes(); ++node_id) { + if (not node_is_owned[node_id]) { + is_empty &= (stencil_array[node_id].size() == 0); + for (size_t i_symmetry_stencil = 0; + i_symmetry_stencil < stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry_stencil) { + is_empty &= + (stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray()[node_id].size() == + 0); + } + } + } + + return is_empty; + }; + + auto are_symmetry_stencils_valid = + []<ItemType source_item_type, ItemType connecting_item_type>(const auto& stencil_array, const auto& mesh, + const size_t number_of_layers) { + bool are_valid_symmetries = true; + + constexpr const size_t Dimension = std::decay_t<decltype(mesh)>::Dimension; + + auto connecting_to_cell_matrix = + mesh.connectivity().template getItemToItemMatrix<connecting_item_type, ItemType::cell>(); + auto cell_to_connecting_item_matrix = + mesh.connectivity().template getItemToItemMatrix<ItemType::cell, connecting_item_type>(); + // auto source_item_is_owned = mesh.connectivity().cellIsOwned(); + auto cell_number = mesh.connectivity().cellNumber(); + auto source_item_is_owned = mesh.connectivity().template isOwned<source_item_type>(); + + auto connecting_number = mesh.connectivity().template number<connecting_item_type>(); + + auto source_item_to_connecting_item_matrix = + [](const auto& given_connectivity) -> ItemToItemMatrix<source_item_type, connecting_item_type> { + if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) == + ItemTypeId<Dimension>::itemTId(connecting_item_type)) { + return ConnectivityMatrix{}; + } else { + return given_connectivity.template getItemToItemMatrix<source_item_type, connecting_item_type>(); + } + }(mesh.connectivity()); + + using SourceItemId = ItemIdT<source_item_type>; + using ConnectingItemId = ItemIdT<connecting_item_type>; + using MeshType = std::decay_t<decltype(mesh)>; + + for (size_t i_symmetry_stencil = 0; + i_symmetry_stencil < stencil_array.symmetryBoundaryStencilArrayList().size(); ++i_symmetry_stencil) { + const IBoundaryDescriptor& boundary_descriptor = + stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].boundaryDescriptor(); + + auto boundary_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray(); + auto boundary_face_list = getMeshFaceBoundary(mesh, boundary_descriptor); + + std::set<ConnectingItemId> sym_connecting_item_set; + + if constexpr (ItemTypeId<MeshType::Dimension>::itemTId(connecting_item_type) == + ItemTypeId<MeshType::Dimension>::itemTId(ItemType::face)) { + for (size_t i_face = 0; i_face < boundary_face_list.faceList().size(); ++i_face) { + const FaceId face_id = boundary_face_list.faceList()[i_face]; + sym_connecting_item_set.insert(ConnectingItemId{FaceId::base_type{face_id}}); + } + + } else { + auto face_to_connecting_item_matrix = + mesh.connectivity().template getItemToItemMatrix<ItemType::face, connecting_item_type>(); + + for (size_t i_face = 0; i_face < boundary_face_list.faceList().size(); ++i_face) { + const FaceId face_id = boundary_face_list.faceList()[i_face]; + auto connecting_item_list = face_to_connecting_item_matrix[face_id]; + for (size_t i_connecting_item = 0; i_connecting_item < connecting_item_list.size(); ++i_connecting_item) { + sym_connecting_item_set.insert(connecting_item_list[i_connecting_item]); + } + } + } + + for (SourceItemId source_item_id = 0; source_item_id < mesh.template numberOf<source_item_type>(); + ++source_item_id) { + if (not source_item_is_owned[source_item_id]) { + are_valid_symmetries &= (boundary_stencil[source_item_id].size() == 0); + } else { + std::vector<CellId> expected_stencil; + + std::set<CellId> marked_cell_set; + if constexpr (source_item_type == ItemType::cell) { + marked_cell_set.insert(source_item_id); + } + + std::set<ConnectingItemId> marked_connecting_item_set; + std::set<ConnectingItemId> layer_connecting_item_set; + if constexpr (ItemTypeId<Dimension>::itemTId(source_item_type) == + ItemTypeId<Dimension>::itemTId(connecting_item_type)) { + const ConnectingItemId connecting_item_id = + static_cast<typename SourceItemId::base_type>(source_item_id); + layer_connecting_item_set.insert(connecting_item_id); + marked_connecting_item_set.insert(connecting_item_id); + } else { + auto source_item_to_connecting_item_list = source_item_to_connecting_item_matrix[source_item_id]; + for (size_t i_connecting_item = 0; i_connecting_item < source_item_to_connecting_item_list.size(); + ++i_connecting_item) { + const ConnectingItemId connecting_item_id = source_item_to_connecting_item_list[i_connecting_item]; + layer_connecting_item_set.insert(connecting_item_id); + marked_connecting_item_set.insert(connecting_item_id); + } + } + + std::set<ConnectingItemId> marked_sym_connecting_item_set; + std::set<ConnectingItemId> layer_sym_connecting_item_set; + + for (auto&& connecting_item_id : marked_connecting_item_set) { + if (sym_connecting_item_set.contains(connecting_item_id)) { + marked_sym_connecting_item_set.insert(connecting_item_id); + layer_sym_connecting_item_set.insert(connecting_item_id); + } + } + + std::set<CellId> marked_sym_cell_set; + + for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) { + 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]; }); + + for (auto&& connecting_item_id : layer_connecting_item_set) { + auto connecting_item_to_cell_list = connecting_to_cell_matrix[connecting_item_id]; + for (size_t i_connecting_item_of_cell = 0; + i_connecting_item_of_cell < connecting_item_to_cell_list.size(); ++i_connecting_item_of_cell) { + const CellId connecting_item_id_of_cell = connecting_item_to_cell_list[i_connecting_item_of_cell]; + if (not marked_cell_set.contains(connecting_item_id_of_cell)) { + cell_set.insert(connecting_item_id_of_cell); + marked_cell_set.insert(connecting_item_id_of_cell); + } + } + } + + std::set<CellId, std::function<bool(CellId, CellId)>> sym_cell_set( + [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; }); + + for (auto&& connecting_item_id : layer_sym_connecting_item_set) { + auto connecting_item_to_cell_list = connecting_to_cell_matrix[connecting_item_id]; + for (size_t i_connecting_item_of_cell = 0; + i_connecting_item_of_cell < connecting_item_to_cell_list.size(); ++i_connecting_item_of_cell) { + const CellId connecting_item_id_of_cell = connecting_item_to_cell_list[i_connecting_item_of_cell]; + if (not marked_sym_cell_set.contains(connecting_item_id_of_cell)) { + sym_cell_set.insert(connecting_item_id_of_cell); + marked_sym_cell_set.insert(connecting_item_id_of_cell); + } + } + } + + for (auto&& set_sym_cell_id : sym_cell_set) { + expected_stencil.push_back(set_sym_cell_id); + } + + layer_connecting_item_set.clear(); + for (auto&& layer_cell_id : cell_set) { + auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[layer_cell_id]; + for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size(); + ++i_connecting_item) { + const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item]; + if (not marked_connecting_item_set.contains(connecting_item_id)) { + layer_connecting_item_set.insert(connecting_item_id); + marked_connecting_item_set.insert(connecting_item_id); + } + } + } + + layer_sym_connecting_item_set.clear(); + + for (auto&& connecting_item_id : layer_connecting_item_set) { + if (sym_connecting_item_set.contains(connecting_item_id)) { + if (not marked_sym_connecting_item_set.contains(connecting_item_id)) { + marked_sym_connecting_item_set.insert(connecting_item_id); + layer_sym_connecting_item_set.insert(connecting_item_id); + } + } + } + + for (auto layer_sym_cell_id : sym_cell_set) { + auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[layer_sym_cell_id]; + for (size_t i_connecting_item = 0; i_connecting_item < cell_to_connecting_item_list.size(); + ++i_connecting_item) { + const ConnectingItemId connecting_item_id = cell_to_connecting_item_list[i_connecting_item]; + if (not marked_sym_connecting_item_set.contains(connecting_item_id)) { + marked_sym_connecting_item_set.insert(connecting_item_id); + layer_sym_connecting_item_set.insert(connecting_item_id); + } + } + } + } + + auto cell_stencil = boundary_stencil[source_item_id]; + + if (cell_stencil.size() != expected_stencil.size()) { + are_valid_symmetries = false; + } + + auto i_set_cell = expected_stencil.begin(); + for (size_t index = 0; index < cell_stencil.size(); ++index, ++i_set_cell) { + if (*i_set_cell != cell_stencil[index]) { + are_valid_symmetries = false; + } + } + } + } + } + + return are_valid_symmetries; + }; + + SECTION("1 layer") + { + SECTION("1D") + { + StencilManager::BoundaryDescriptorList list; + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX")); + + SECTION("cartesian") + { + const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>(); + + const Connectivity<1>& connectivity = mesh.connectivity(); + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1)); + } + } + + SECTION("unordered") + { + const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>(); + + const Connectivity<1>& connectivity = mesh.connectivity(); + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1)); + } + } + } + + SECTION("2D") + { + StencilManager::BoundaryDescriptorList list; + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX")); + + SECTION("cartesian") + { + const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>(); + + const Connectivity<2>& connectivity = mesh.connectivity(); + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 1)); + } + } + + SECTION("hybrid") + { + const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>(); + + const Connectivity<2>& connectivity = mesh.connectivity(); + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 1)); + } + } + } + + SECTION("3D") + { + StencilManager::BoundaryDescriptorList list; + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMAX")); + + SECTION("cartesian") + { + const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>(); + + const Connectivity<3>& connectivity = mesh.connectivity(); + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 1)); + } + } + + SECTION("hybrid") + { + const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); + + const Connectivity<3>& connectivity = mesh.connectivity(); + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 1)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 1)); + } + } + } + } + + SECTION("2 layers") + { + NbGhostLayersTester nb_ghost_layers_tester(2); + + SECTION("1D") + { + StencilManager::BoundaryDescriptorList list; + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX")); + + SECTION("cartesian") + { + auto mesh_v = CartesianMeshBuilder(TinyVector<1>{0}, TinyVector<1>{1}, TinyVector<1, uint64_t>(20)).mesh(); + const auto& mesh = *mesh_v->get<Mesh<1>>(); + + const Connectivity<1>& connectivity = mesh.connectivity(); + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 2)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 2)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 2)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2)); + } + } + } + + SECTION("2D") + { + StencilManager::BoundaryDescriptorList list; + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX")); + + SECTION("cartesian") + { + auto mesh_v = + CartesianMeshBuilder(TinyVector<2>{0, 0}, TinyVector<2>{1, 2}, TinyVector<2, uint64_t>(5, 7)).mesh(); + const auto& mesh = *mesh_v->get<Mesh<2>>(); + + const Connectivity<2>& connectivity = mesh.connectivity(); + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 2)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 2)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 2)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 2)); + REQUIRE(not(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2))); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 2)); + } + } + } + + SECTION("3D") + { + StencilManager::BoundaryDescriptorList list; + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("YMAX")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("ZMAX")); + + SECTION("cartesian") + { + auto mesh_v = + CartesianMeshBuilder(TinyVector<3>{0, 0, 0}, TinyVector<3>{1, 1, 2}, TinyVector<3, uint64_t>(3, 4, 5)) + .mesh(); + const auto& mesh = *mesh_v->get<Mesh<3>>(); + + const Connectivity<3>& connectivity = mesh.connectivity(); + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::node>(stencil_array, mesh, 2)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::node>(connectivity, stencil_array, 2)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::edge>(stencil_array, mesh, 2)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::edge>(connectivity, stencil_array, 2)); + } + + { + auto stencil_array = + StencilManager::instance() + .getNodeToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_nodes_have_empty_stencils(stencil_array, connectivity)); + REQUIRE( + are_symmetry_stencils_valid.template operator()<ItemType::node, ItemType::face>(stencil_array, mesh, 2)); + REQUIRE(is_valid.template operator()<ItemType::node, ItemType::face>(connectivity, stencil_array, 2)); + } + } + } + } + } +}