diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp index 31db762a1002e3c633afeaabd4f06d6a98dff634..96f36a14f99ce3a3239f94d1b454f35ad82fc22b 100644 --- a/src/mesh/StencilBuilder.cpp +++ b/src/mesh/StencilBuilder.cpp @@ -5,16 +5,23 @@ #include <utils/GlobalVariableManager.hpp> #include <utils/Messenger.hpp> -#include <set> - template <ItemType item_type> class StencilBuilder::Layer { using ItemId = ItemIdT<item_type>; + ItemValue<const int, item_type> m_number_of; + std::vector<ItemId> m_item_id_vector; std::vector<int> m_item_number_vector; public: + void + clear() + { + m_item_id_vector.clear(); + m_item_number_vector.clear(); + } + size_t size() const { @@ -22,9 +29,17 @@ class StencilBuilder::Layer return m_item_id_vector.size(); } + const std::vector<ItemId>& + itemIdList() const + { + return m_item_id_vector; + } + bool - hasItemNumber(const int item_number) const + hasItemNumber(const ItemId item_id) const { + const int item_number = m_number_of[item_id]; + ssize_t begin = 0; ssize_t end = m_item_number_vector.size(); @@ -37,26 +52,22 @@ class StencilBuilder::Layer if (begin == mid) { break; } - begin = mid - 1; + begin = mid; } else { if (end == mid) { break; } - end = mid + 1; + end = mid; } } return false; } - const std::vector<ItemId>& - itemIdList() const - { - return m_item_id_vector; - } - void - add(const ItemId item_id, const int item_number) + add(const ItemId item_id) { + const int item_number = m_number_of[item_id]; + ssize_t begin = 0; ssize_t end = m_item_number_vector.size(); @@ -97,7 +108,14 @@ class StencilBuilder::Layer } } - Layer() = default; + Layer& operator=(const Layer&) = default; + Layer& operator=(Layer&&) = default; + + template <size_t Dimension> + Layer(const Connectivity<Dimension>& connectivity) : m_number_of{connectivity.template number<item_type>()} + {} + + // Layer() = default; Layer(const Layer&) = default; Layer(Layer&&) = default; @@ -117,41 +135,8 @@ StencilBuilder::_buildSymmetryConnectingItemList(const ConnectivityType& connect symmetry_boundary_descriptor_list.size()}; symmetry_connecting_item_list.fill(false); - if constexpr (ConnectivityType::Dimension > 1) { - auto face_to_connecting_item_matrix = - connectivity.template getItemToItemMatrix<ItemType::face, connecting_item_type>(); - size_t i_symmetry_boundary = 0; - for (auto p_boundary_descriptor : symmetry_boundary_descriptor_list) { - const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor; - - bool found = false; - for (size_t i_ref_face_list = 0; i_ref_face_list < connectivity.template numberOfRefItemList<ItemType::face>(); - ++i_ref_face_list) { - const auto& ref_face_list = connectivity.template refItemList<ItemType::face>(i_ref_face_list); - if (ref_face_list.refId() == boundary_descriptor) { - found = true; - for (size_t i_face = 0; i_face < ref_face_list.list().size(); ++i_face) { - const FaceId face_id = ref_face_list.list()[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) { - const ConnectingItemId connecting_item_id = connecting_item_list[i_connecting_item]; - - symmetry_connecting_item_list[connecting_item_id][i_symmetry_boundary] = true; - } - } - break; - } - } - ++i_symmetry_boundary; - if (not found) { - std::ostringstream error_msg; - error_msg << "cannot find boundary '" << rang::fgB::yellow << boundary_descriptor << rang::fg::reset << '\''; - throw NormalError(error_msg.str()); - } - } - - return symmetry_connecting_item_list; - } else { + if constexpr (ItemTypeId<ConnectivityType::Dimension>::itemTId(connecting_item_type) == + ItemTypeId<ConnectivityType::Dimension>::itemTId(ItemType::face)) { size_t i_symmetry_boundary = 0; for (auto p_boundary_descriptor : symmetry_boundary_descriptor_list) { const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor; @@ -180,110 +165,41 @@ StencilBuilder::_buildSymmetryConnectingItemList(const ConnectivityType& connect throw NormalError(error_msg.str()); } } + } else { + auto face_to_connecting_item_matrix = + connectivity.template getItemToItemMatrix<ItemType::face, connecting_item_type>(); + size_t i_symmetry_boundary = 0; + for (auto p_boundary_descriptor : symmetry_boundary_descriptor_list) { + const IBoundaryDescriptor& boundary_descriptor = *p_boundary_descriptor; - return symmetry_connecting_item_list; - } -} - -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) { - const auto 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) { - const size_t i_source = i_destination - 1; + bool found = false; + for (size_t i_ref_face_list = 0; i_ref_face_list < connectivity.template numberOfRefItemList<ItemType::face>(); + ++i_ref_face_list) { + const auto& ref_face_list = connectivity.template refItemList<ItemType::face>(i_ref_face_list); + if (ref_face_list.refId() == boundary_descriptor) { + found = true; + for (size_t i_face = 0; i_face < ref_face_list.list().size(); ++i_face) { + const FaceId face_id = ref_face_list.list()[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) { + const ConnectingItemId connecting_item_id = connecting_item_list[i_connecting_item]; - column_indices[i_destination] = column_indices[i_source]; - } - ++max_index[cell_id]; - column_indices[i_index] = node_cell_id; + symmetry_connecting_item_list[connecting_item_id][i_symmetry_boundary] = true; } } + break; } } + ++i_symmetry_boundary; + if (not found) { + std::ostringstream error_msg; + error_msg << "cannot find boundary '" << rang::fgB::yellow << boundary_descriptor << rang::fg::reset << '\''; + throw NormalError(error_msg.str()); + } } } - return column_indices; + return symmetry_connecting_item_list; } template <ItemType item_type, ItemType connecting_item_type, typename ConnectivityType> @@ -295,9 +211,7 @@ StencilBuilder::_build_for_same_source_and_target(const ConnectivityType& connec if (number_of_layers == 0) { throw NormalError("number of layers must be greater than 0 to build stencils"); } - if (number_of_layers > 2) { - throw NotImplementedError("number of layers too large"); - } + auto item_to_connecting_item_matrix = connectivity.template getItemToItemMatrix<item_type, connecting_item_type>(); auto connecting_item_to_item_matrix = connectivity.template getItemToItemMatrix<connecting_item_type, item_type>(); @@ -308,215 +222,250 @@ StencilBuilder::_build_for_same_source_and_target(const ConnectivityType& connec using ItemId = ItemIdT<item_type>; using ConnectingItemId = ItemIdT<connecting_item_type>; - if (symmetry_boundary_descriptor_list.size() == 0) { - if (number_of_layers == 1) { - Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1}; - row_map[0] = 0; + ItemArray<bool, connecting_item_type> symmetry_item_list = + this->_buildSymmetryConnectingItemList<connecting_item_type>(connectivity, symmetry_boundary_descriptor_list); - std::vector<ItemId> column_indices_vector; + Array<uint32_t> row_map{connectivity.template numberOf<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<item_type>() + 1}; + symmetry_row_map[0] = 0; + } - for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) { - // First layer is a special case + std::vector<ItemId> column_indices_vector; + std::vector<std::vector<uint32_t>> symmetry_column_indices_vector(symmetry_boundary_descriptor_list.size()); - Layer<item_type> item_layer; - Layer<connecting_item_type> connecting_layer; + std::vector<Layer<item_type>> item_layer_list; + std::vector<Layer<connecting_item_type>> connecting_layer_list; - if (item_is_owned[item_id]) { - for (size_t i_connecting_item_1 = 0; i_connecting_item_1 < item_to_connecting_item_matrix[item_id].size(); - ++i_connecting_item_1) { - const ConnectingItemId layer_1_connecting_item_id = - item_to_connecting_item_matrix[item_id][i_connecting_item_1]; - connecting_layer.add(layer_1_connecting_item_id, connecting_item_number[layer_1_connecting_item_id]); - } + std::vector<Layer<item_type>> symmetry_item_layer_list; + std::vector<Layer<connecting_item_type>> symmetry_connecting_layer_list; - for (auto connecting_item_id : connecting_layer.itemIdList()) { - for (size_t i_item_1 = 0; i_item_1 < connecting_item_to_item_matrix[connecting_item_id].size(); - ++i_item_1) { - const ItemId layer_1_item_id = connecting_item_to_item_matrix[connecting_item_id][i_item_1]; - if (layer_1_item_id != item_id) { - item_layer.add(layer_1_item_id, item_number[layer_1_item_id]); - } - } - } - } + for (size_t i = 0; i < number_of_layers; ++i) { + item_layer_list.emplace_back(Layer<item_type>{connectivity}); + connecting_layer_list.emplace_back(Layer<connecting_item_type>{connectivity}); - for (auto layer_item_id : item_layer.itemIdList()) { - column_indices_vector.push_back(layer_item_id); - } - row_map[item_id + 1] = row_map[item_id] + item_layer.itemIdList().size(); - } + if (symmetry_boundary_descriptor_list.size() > 0) { + symmetry_item_layer_list.emplace_back(Layer<item_type>{connectivity}); + symmetry_connecting_layer_list.emplace_back(Layer<connecting_item_type>{connectivity}); + } + } - if (row_map[row_map.size() - 1] != column_indices_vector.size()) { - throw UnexpectedError("incorrect stencil size"); + for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) { + if (item_is_owned[item_id]) { + for (auto&& item_layer : item_layer_list) { + item_layer.clear(); + } + for (auto&& connecting_layer : connecting_layer_list) { + connecting_layer.clear(); } - Array<uint32_t> column_indices(row_map[row_map.size() - 1]); - column_indices.fill(std::numeric_limits<uint32_t>::max()); - for (size_t i = 0; i < column_indices.size(); ++i) { - column_indices[i] = column_indices_vector[i]; + // First layer is a special case + { + auto& item_layer = item_layer_list[0]; + auto& connecting_layer = connecting_layer_list[0]; + + for (size_t i_connecting_item = 0; i_connecting_item < item_to_connecting_item_matrix[item_id].size(); + ++i_connecting_item) { + const ConnectingItemId connecting_item_id = item_to_connecting_item_matrix[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_item_matrix[connecting_item_id].size(); ++i_item) { + const ItemId layer_item_id = connecting_item_to_item_matrix[connecting_item_id][i_item]; + if (layer_item_id != item_id) { + item_layer.add(layer_item_id); + } + } + } } - return {ConnectivityMatrix{row_map, column_indices}, {}}; - } else { - Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1}; - row_map[0] = 0; - - std::vector<ItemId> column_indices_vector; - - for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) { - if (item_is_owned[item_id]) { - std::set<ItemId, std::function<bool(ItemId, ItemId)>> item_set( - [=](ItemId item_0, ItemId item_1) { return item_number[item_0] < item_number[item_1]; }); - - for (size_t i_connecting_item_1 = 0; i_connecting_item_1 < item_to_connecting_item_matrix[item_id].size(); - ++i_connecting_item_1) { - const ConnectingItemId layer_1_connecting_item_id = - item_to_connecting_item_matrix[item_id][i_connecting_item_1]; - - for (size_t i_item_1 = 0; i_item_1 < connecting_item_to_item_matrix[layer_1_connecting_item_id].size(); - ++i_item_1) { - ItemId item_1_id = connecting_item_to_item_matrix[layer_1_connecting_item_id][i_item_1]; - - for (size_t i_connecting_item_2 = 0; - i_connecting_item_2 < item_to_connecting_item_matrix[item_1_id].size(); ++i_connecting_item_2) { - const ConnectingItemId layer_2_connecting_item_id = - item_to_connecting_item_matrix[item_1_id][i_connecting_item_2]; - - for (size_t i_item_2 = 0; i_item_2 < connecting_item_to_item_matrix[layer_2_connecting_item_id].size(); - ++i_item_2) { - ItemId item_2_id = connecting_item_to_item_matrix[layer_2_connecting_item_id][i_item_2]; - - if (item_2_id != item_id) { - item_set.insert(item_2_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; } } - for (auto stencil_item_id : item_set) { - column_indices_vector.push_back(stencil_item_id); + return false; + }; + + for (auto&& previous_layer_item_id_list : item_layer_list[i_layer - 1].itemIdList()) { + const auto item_to_connecting_item_list = item_to_connecting_item_matrix[previous_layer_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)) { + connecting_layer.add(connecting_item_id); + } } - row_map[item_id + 1] = row_map[item_id] + item_set.size(); } - } - if (row_map[row_map.size() - 1] != column_indices_vector.size()) { - throw UnexpectedError("incorrect stencil size"); - } + auto& item_layer = item_layer_list[i_layer]; - Array<uint32_t> column_indices(row_map[row_map.size() - 1]); - column_indices.fill(std::numeric_limits<uint32_t>::max()); + auto has_layer_item = [&i_layer, &item_layer_list](ItemId layer_item_id) -> bool { + for (ssize_t i = i_layer - 1; i >= 0; --i) { + if (item_layer_list[i].hasItemNumber(layer_item_id)) { + return true; + } + } - for (size_t i = 0; i < column_indices.size(); ++i) { - column_indices[i] = column_indices_vector[i]; - } - ConnectivityMatrix primal_stencil{row_map, column_indices}; + return false; + }; - return {primal_stencil, {}}; - } - } else { - 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<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<item_type>() + 1}; - symmetry_row_map[0] = 0; - } + for (auto connecting_item_id : connecting_layer.itemIdList()) { + const auto& connecting_item_to_item_list = connecting_item_to_item_matrix[connecting_item_id]; + for (size_t i_item = 0; i_item < connecting_item_to_item_list.size(); ++i_item) { + const ItemId layer_item_id = connecting_item_to_item_list[i_item]; + if ((layer_item_id != item_id) and (not has_layer_item(layer_item_id))) { + item_layer.add(layer_item_id); + } + } + } + } - std::vector<uint32_t> column_indices_vector; - std::vector<std::vector<uint32_t>> symmetry_column_indices_vector(symmetry_boundary_descriptor_list.size()); + for (size_t i_symmetry = 0; i_symmetry < symmetry_boundary_descriptor_list.size(); ++i_symmetry) { + for (auto&& symmetry_item_layer : symmetry_item_layer_list) { + symmetry_item_layer.clear(); + } + for (auto&& symmetry_connecting_layer : symmetry_connecting_layer_list) { + symmetry_connecting_layer.clear(); + } - for (ItemId item_id = 0; item_id < connectivity.template numberOf<item_type>(); ++item_id) { - std::set<ItemId> item_set; - std::vector<std::set<ItemId>> by_boundary_symmetry_item(symmetry_boundary_descriptor_list.size()); + // 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); + } + } - if (item_is_owned[item_id]) { - auto item_to_connecting_item_list = item_to_connecting_item_matrix[item_id]; - for (size_t i_connecting_item_of_item = 0; i_connecting_item_of_item < item_to_connecting_item_list.size(); - ++i_connecting_item_of_item) { - const ConnectingItemId connecting_item_id_of_item = item_to_connecting_item_list[i_connecting_item_of_item]; - auto connecting_item_to_item_list = connecting_item_to_item_matrix[connecting_item_id_of_item]; - for (size_t i_item_of_connecting_item = 0; i_item_of_connecting_item < connecting_item_to_item_list.size(); + for (auto&& connecting_item_id : symmetry_connecting_layer_list[0].itemIdList()) { + auto item_of_connecting_item_list = connecting_item_to_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 ItemId item_id_of_connecting_item = connecting_item_to_item_list[i_item_of_connecting_item]; - if (item_id != item_id_of_connecting_item) { - item_set.insert(item_id_of_connecting_item); - } + const ItemId item_id_of_connecting_item = item_of_connecting_item_list[i_item_of_connecting_item]; + symmetry_item_layer_list[0].add(item_id_of_connecting_item); } } - { - std::vector<ItemId> item_vector; - for (auto&& set_item_id : item_set) { - item_vector.push_back(set_item_id); + 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); + } + } } - std::sort(item_vector.begin(), item_vector.end(), - [&item_number](const ItemId& item0_id, const ItemId& item1_id) { - return item_number[item0_id] < item_number[item1_id]; - }); - for (auto&& vector_item_id : item_vector) { - column_indices_vector.push_back(vector_item_id); + for (auto&& previous_layer_item_id_list : symmetry_item_layer_list[i_layer - 1].itemIdList()) { + const auto item_to_connecting_item_list = item_to_connecting_item_matrix[previous_layer_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); + } + } } - } - for (size_t i = 0; i < symmetry_boundary_descriptor_list.size(); ++i) { - std::set<ItemId> symmetry_item_set; - for (size_t i_connecting_item_of_item = 0; i_connecting_item_of_item < item_to_connecting_item_list.size(); - ++i_connecting_item_of_item) { - const ConnectingItemId connecting_item_id_of_item = item_to_connecting_item_list[i_connecting_item_of_item]; - if (symmetry_item_list[connecting_item_id_of_item][i]) { - auto item_of_connecting_item_list = connecting_item_to_item_matrix[connecting_item_id_of_item]; - 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 ItemId item_id_of_connecting_item = item_of_connecting_item_list[i_item_of_connecting_item]; - symmetry_item_set.insert(item_id_of_connecting_item); + auto has_symmetry_layer_item = [&i_layer, &symmetry_item_layer_list](ItemId layer_item_id) -> bool { + for (ssize_t i = i_layer - 1; i >= 0; --i) { + if (symmetry_item_layer_list[i].hasItemNumber(layer_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_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 ItemId item_id_of_connecting_item = item_of_connecting_item_list[i_item_of_connecting_item]; + if (not has_symmetry_layer_item(item_id_of_connecting_item)) { + symmetry_item_layer_list[i_layer].add(item_id_of_connecting_item); } } } - by_boundary_symmetry_item[i] = symmetry_item_set; + } - std::vector<ItemId> item_vector; - for (auto&& set_item_id : symmetry_item_set) { - item_vector.push_back(set_item_id); + for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) { + for (auto symmetry_layer_item_id : symmetry_item_layer_list[i_layer].itemIdList()) { + symmetry_column_indices_vector[i_symmetry].push_back(symmetry_layer_item_id); } - std::sort(item_vector.begin(), item_vector.end(), - [&item_number](const ItemId& item0_id, const ItemId& item1_id) { - return item_number[item0_id] < item_number[item1_id]; - }); + } - for (auto&& vector_item_id : item_vector) { - symmetry_column_indices_vector[i].push_back(vector_item_id); + { + size_t stencil_size = 0; + for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) { + auto& item_layer = symmetry_item_layer_list[i_layer]; + stencil_size += item_layer.size(); } + + symmetry_row_map_list[i_symmetry][item_id + 1] = symmetry_row_map_list[i_symmetry][item_id] + stencil_size; } } - row_map[item_id + 1] = row_map[item_id] + item_set.size(); - for (size_t i = 0; i < symmetry_row_map_list.size(); ++i) { - symmetry_row_map_list[i][item_id + 1] = symmetry_row_map_list[i][item_id] + by_boundary_symmetry_item[i].size(); + { + size_t stencil_size = 0; + for (size_t i_layer = 0; i_layer < number_of_layers; ++i_layer) { + auto& item_layer = 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[item_id + 1] = row_map[item_id] + stencil_size; } - } - ConnectivityMatrix primal_stencil{row_map, convert_to_array(column_indices_vector)}; - - typename StencilArray<item_type, 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<item_type, item_type>:: - BoundaryDescriptorStencilArray{p_boundary_descriptor, - ConnectivityMatrix{symmetry_row_map_list[i], - convert_to_array(symmetry_column_indices_vector[i])}}); - ++i; + + } else { + row_map[item_id + 1] = row_map[item_id]; + for (size_t i = 0; i < symmetry_row_map_list.size(); ++i) { + symmetry_row_map_list[i][item_id + 1] = symmetry_row_map_list[i][item_id]; } } + } - return {{primal_stencil}, {symmetry_boundary_stencil_list}}; + 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<item_type, 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<item_type, 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, diff --git a/src/mesh/StencilBuilder.hpp b/src/mesh/StencilBuilder.hpp index 4f7ead00b618db999dd5386258ac058b664c65a7..f504cde9069d3966658b17f51a884acc5fb9cde5 100644 --- a/src/mesh/StencilBuilder.hpp +++ b/src/mesh/StencilBuilder.hpp @@ -19,13 +19,6 @@ class StencilBuilder template <ItemType item_type> class Layer; - 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; - template <ItemType connecting_item, typename ConnectivityType> auto _buildSymmetryConnectingItemList(const ConnectivityType& connectivity, const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; @@ -74,6 +67,9 @@ class StencilBuilder const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; friend class StencilManager; + + public: +#warning REMAKE THIS FUNCTION PRIVATE template <ItemType source_item_type, ItemType target_item_type> StencilArray<source_item_type, target_item_type> build(const IConnectivity& connectivity, diff --git a/tests/NbGhostLayersTester.hpp b/tests/NbGhostLayersTester.hpp new file mode 100644 index 0000000000000000000000000000000000000000..0468fe18305367834c9cbd0959ba5b3f443fe870 --- /dev/null +++ b/tests/NbGhostLayersTester.hpp @@ -0,0 +1,27 @@ +#ifndef NB_GHOST_LAYERS_TESTER_HPP +#define NB_GHOST_LAYERS_TESTER_HPP + +#include <cstddef> +#include <utils/GlobalVariableManager.hpp> + +class NbGhostLayersTester +{ + private: + const size_t m_original_number_of_ghost_layers; + + public: + PUGS_INLINE + NbGhostLayersTester(const size_t number_of_ghost_layers) + : m_original_number_of_ghost_layers{GlobalVariableManager::instance().getNumberOfGhostLayers()} + { + GlobalVariableManager::instance().m_number_of_ghost_layers = number_of_ghost_layers; + } + + PUGS_INLINE + ~NbGhostLayersTester() + { + GlobalVariableManager::instance().m_number_of_ghost_layers = m_original_number_of_ghost_layers; + } +}; + +#endif // NB_GHOST_LAYERS_TESTER_HPP diff --git a/tests/test_ConnectivityDispatcher.cpp b/tests/test_ConnectivityDispatcher.cpp index abe07cc21e21a1b34ae3f20dfe13ff7f52ac83b9..09ec9982e10e08d84c808f157fe4e11add3d756a 100644 --- a/tests/test_ConnectivityDispatcher.cpp +++ b/tests/test_ConnectivityDispatcher.cpp @@ -9,29 +9,12 @@ #include <utils/Messenger.hpp> #include <MeshDataBaseForTests.hpp> +#include <NbGhostLayersTester.hpp> #include <filesystem> // clazy:excludeall=non-pod-global-static -class NbGhostLayersTester -{ - private: - const size_t m_original_number_of_ghost_layers; - - public: - NbGhostLayersTester(const size_t number_of_ghost_layers) - : m_original_number_of_ghost_layers{GlobalVariableManager::instance().getNumberOfGhostLayers()} - { - GlobalVariableManager::instance().m_number_of_ghost_layers = number_of_ghost_layers; - } - - ~NbGhostLayersTester() - { - GlobalVariableManager::instance().m_number_of_ghost_layers = m_original_number_of_ghost_layers; - } -}; - TEST_CASE("ConnectivityDispatcher", "[mesh]") { auto check_number_of_ghost_layers = [](const auto& connectivity, const size_t number_of_layers) { diff --git a/tests/test_StencilBuilder.cpp b/tests/test_StencilBuilder.cpp index b99ba218aa1e3eb5d8ec884e974cfad2f4a157db..8759b36c5c8dbf5e1b1d6c8af9b72ee4b7fe940f 100644 --- a/tests/test_StencilBuilder.cpp +++ b/tests/test_StencilBuilder.cpp @@ -2,138 +2,420 @@ #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/MeshFlatNodeBoundary.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", "[mesh]") { - SECTION("inner stencil") - { - auto is_valid = [](const auto& connectivity, const auto& stencil_array) { - auto cell_to_node_matrix = connectivity.cellToNodeMatrix(); - auto node_to_cell_matrix = connectivity.nodeToCellMatrix(); + auto is_valid = []<ItemType connecting_item_type>(const auto& connectivity, const auto& stencil_array, + const size_t number_of_layers) { + auto cell_to_connecting_item_matrix = + connectivity.template getItemToItemMatrix<ItemType::cell, connecting_item_type>(); + auto connecting_to_cell_matrix = connectivity.template getItemToItemMatrix<connecting_item_type, ItemType::cell>(); + auto cell_is_owned = connectivity.cellIsOwned(); + auto cell_number = connectivity.cellNumber(); + + using ConnectingItemId = ItemIdT<connecting_item_type>; - 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::vector<CellId> expected_stencil; + + std::set<CellId> marked_cell_set; + marked_cell_set.insert(cell_id); + + std::set<ConnectingItemId> marked_connecting_item_set; + std::set<ConnectingItemId> layer_connecting_item_set; + { + auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[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]; + layer_connecting_item_set.insert(connecting_item_id); + marked_connecting_item_set.insert(connecting_item_id); + } + } - for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) { - if (cell_is_owned[cell_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]; }); - 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); + + 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); } } } - auto cell_stencil = stencil_array[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; + 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[cell_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; - }; + } + return true; + }; - SECTION("1D") + auto are_stencil_by_nodes_valid = [&](const auto& connectivity, const auto& stencil_array, + const size_t number_of_layers) { + return is_valid.template operator()<ItemType::node>(connectivity, stencil_array, number_of_layers); + }; + + auto are_stencil_by_faces_valid = [&](const auto& connectivity, const auto& stencil_array, + const size_t number_of_layers) { + return is_valid.template operator()<ItemType::face>(connectivity, stencil_array, number_of_layers); + }; + + auto are_stencil_by_edges_valid = [&](const auto& connectivity, const auto& stencil_array, + const size_t number_of_layers) { + return is_valid.template operator()<ItemType::edge>(connectivity, stencil_array, number_of_layers); + }; + + SECTION("inner stencil") + { + SECTION("1 layer") { - SECTION("cartesian") + SECTION("1D") { - const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>(); + SECTION("cartesian") + { + const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>(); - const Connectivity<1>& connectivity = mesh.connectivity(); + const Connectivity<1>& connectivity = mesh.connectivity(); - auto stencil_array = - StencilManager::instance() - .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}); + REQUIRE( + are_stencil_by_nodes_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); - REQUIRE(is_valid(connectivity, stencil_array)); - } + REQUIRE( + are_stencil_by_edges_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); - SECTION("unordered") - { - const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>(); + REQUIRE( + are_stencil_by_faces_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(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(); - const Connectivity<1>& connectivity = mesh.connectivity(); - auto stencil_array = - StencilManager::instance() - .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}); + REQUIRE( + are_stencil_by_nodes_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); - REQUIRE(is_valid(connectivity, stencil_array)); + REQUIRE( + are_stencil_by_edges_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + are_stencil_by_faces_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_faces}), + 1)); + } } - } - SECTION("2D") - { - SECTION("cartesian") + SECTION("2D") { - const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>(); - - const Connectivity<2>& connectivity = mesh.connectivity(); - REQUIRE( - is_valid(connectivity, - StencilManager::instance() - .getCellToCellStencilArray(connectivity, - StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}))); + SECTION("cartesian") + { + const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>(); + + const Connectivity<2>& connectivity = mesh.connectivity(); + REQUIRE( + are_stencil_by_nodes_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); + + REQUIRE( + are_stencil_by_edges_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + are_stencil_by_faces_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(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( + are_stencil_by_nodes_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); + + REQUIRE( + are_stencil_by_edges_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + are_stencil_by_faces_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_faces}), + 1)); + } } - SECTION("hybrid") + SECTION("3D") { - const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>(); - - const Connectivity<2>& connectivity = mesh.connectivity(); - REQUIRE( - is_valid(connectivity, - StencilManager::instance() - .getCellToCellStencilArray(connectivity, - StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}))); + SECTION("carteian") + { + const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>(); + + const Connectivity<3>& connectivity = mesh.connectivity(); + REQUIRE( + are_stencil_by_nodes_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); + + REQUIRE( + are_stencil_by_edges_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + are_stencil_by_faces_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(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( + are_stencil_by_nodes_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_nodes}), + 1)); + + REQUIRE( + are_stencil_by_edges_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_edges}), + 1)); + + REQUIRE( + are_stencil_by_faces_valid(connectivity, + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor:: + ConnectionType::by_faces}), + 1)); + } } } - SECTION("3D") + SECTION("2 layers") { - SECTION("carteian") + NbGhostLayersTester nb_ghost_layers_tester(2); + + SECTION("1D") { - const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>(); - - const Connectivity<3>& connectivity = mesh.connectivity(); - REQUIRE( - is_valid(connectivity, - StencilManager::instance() - .getCellToCellStencilArray(connectivity, - StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}))); + 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( + are_stencil_by_nodes_valid(mesh.connectivity(), + StencilManager::instance() + .getCellToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_nodes}), + 2)); + + REQUIRE( + are_stencil_by_edges_valid(mesh.connectivity(), + StencilManager::instance() + .getCellToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_edges}), + 2)); + + REQUIRE( + are_stencil_by_faces_valid(mesh.connectivity(), + StencilManager::instance() + .getCellToCellStencilArray(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( + are_stencil_by_nodes_valid(mesh.connectivity(), + StencilManager::instance() + .getCellToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_nodes}), + 2)); + + REQUIRE( + are_stencil_by_edges_valid(mesh.connectivity(), + StencilManager::instance() + .getCellToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_edges}), + 2)); + + REQUIRE( + are_stencil_by_faces_valid(mesh.connectivity(), + StencilManager::instance() + .getCellToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_faces}), + 2)); + } } - SECTION("hybrid") + SECTION("3D") { - const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); - - const Connectivity<3>& connectivity = mesh.connectivity(); - REQUIRE( - is_valid(connectivity, - StencilManager::instance() - .getCellToCellStencilArray(connectivity, - StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}))); + 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( + are_stencil_by_nodes_valid(mesh.connectivity(), + StencilManager::instance() + .getCellToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_nodes}), + 2)); + + REQUIRE( + are_stencil_by_edges_valid(mesh.connectivity(), + StencilManager::instance() + .getCellToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_edges}), + 2)); + + REQUIRE( + are_stencil_by_faces_valid(mesh.connectivity(), + StencilManager::instance() + .getCellToCellStencilArray(mesh.connectivity(), + StencilDescriptor{2, StencilDescriptor:: + ConnectionType::by_faces}), + 2)); + } } } } @@ -160,13 +442,21 @@ TEST_CASE("StencilBuilder", "[mesh]") return is_empty; }; - auto symmetry_stencils_are_valid = [](const auto& stencil_array, const auto& mesh) { - bool is_valid = true; + auto are_symmetry_stencils_valid = []<ItemType connecting_item_type>(const auto& stencil_array, const auto& mesh, + const size_t number_of_layers) { + bool are_valid_symmetries = true; - auto node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); - auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); - auto cell_is_owned = mesh.connectivity().cellIsOwned(); - auto cell_number = mesh.connectivity().cellNumber(); + 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 cell_is_owned = mesh.connectivity().cellIsOwned(); + auto cell_number = mesh.connectivity().cellNumber(); + + auto connecting_number = mesh.connectivity().template number<connecting_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) { @@ -174,175 +464,620 @@ TEST_CASE("StencilBuilder", "[mesh]") stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].boundaryDescriptor(); auto boundary_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry_stencil].stencilArray(); - auto boundary_node_list = getMeshFlatNodeBoundary(mesh, boundary_descriptor); - - CellValue<bool> boundary_cell{mesh.connectivity()}; - boundary_cell.fill(false); - auto node_list = boundary_node_list.nodeList(); - for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - const NodeId node_id = node_list[i_node]; - auto node_cell_list = node_to_cell_matrix[node_id]; - for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { - const CellId cell_id = node_cell_list[i_cell]; - boundary_cell[cell_id] = true; + 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}}); } - } - std::set<NodeId> symmetry_node; - for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - const NodeId node_id = node_list[i_node]; - symmetry_node.insert(node_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 (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) { - if ((not boundary_cell[cell_id]) or (not cell_is_owned[cell_id])) { - is_valid &= (boundary_stencil[cell_id].size() == 0); + if (not cell_is_owned[cell_id]) { + are_valid_symmetries &= (boundary_stencil[cell_id].size() == 0); } else { - auto cell_node_list = cell_to_node_matrix[cell_id]; - std::set<CellId, std::function<bool(CellId, CellId)>> cell_set( - [&](CellId cell0_id, CellId cell1_id) { return cell_number[cell0_id] < cell_number[cell1_id]; }); - for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) { - const NodeId node_id = cell_node_list[i_node]; - if (symmetry_node.contains(node_id)) { - auto node_cell_list = node_to_cell_matrix[node_id]; - for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) { - const CellId node_cell_id = node_cell_list[i_node_cell]; - cell_set.insert(node_cell_id); + std::vector<CellId> expected_stencil; + + std::set<CellId> marked_cell_set; + marked_cell_set.insert(cell_id); + + std::set<ConnectingItemId> marked_connecting_item_set; + std::set<ConnectingItemId> layer_connecting_item_set; + { + auto cell_to_connecting_item_list = cell_to_connecting_item_matrix[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]; + 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); + } } } } - if (cell_set.size() == boundary_stencil[cell_id].size()) { - size_t i = 0; - for (auto&& id : cell_set) { - is_valid &= (id == boundary_stencil[cell_id][i++]); + auto cell_stencil = boundary_stencil[cell_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; } - } else { - is_valid = false; } } } } - return is_valid; + return are_valid_symmetries; + }; + + auto are_symmetry_stencils_by_nodes_valid = [&](const auto& stencil_array, const auto& mesh, + const size_t number_of_layers) { + return are_symmetry_stencils_valid.template operator()<ItemType::node>(stencil_array, mesh, number_of_layers); + }; + + auto are_symmetry_stencils_by_edges_valid = [&](const auto& stencil_array, const auto& mesh, + const size_t number_of_layers) { + return are_symmetry_stencils_valid.template operator()<ItemType::edge>(stencil_array, mesh, number_of_layers); }; - SECTION("1D") + auto are_symmetry_stencils_by_faces_valid = [&](const auto& stencil_array, const auto& mesh, + const size_t number_of_layers) { + return are_symmetry_stencils_valid.template operator()<ItemType::face>(stencil_array, mesh, number_of_layers); + }; + + SECTION("1 layer") { - StencilManager::BoundaryDescriptorList list; - list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN")); - list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX")); + 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() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_nodes_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_edges_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_faces_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_nodes_valid(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() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_nodes_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_edges_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 1)); + } - SECTION("cartesian") + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_faces_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 1)); + } + } + } + + SECTION("2D") { - const auto& mesh = *MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>(); + 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() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_nodes_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_edges_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_edges_valid(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_faces_valid(stencil_array, mesh, 1)); + REQUIRE(not(are_stencil_by_nodes_valid(connectivity, stencil_array, 1))); + REQUIRE(are_stencil_by_faces_valid(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() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_nodes_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 1)); + } - const Connectivity<1>& connectivity = mesh.connectivity(); + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); - auto stencil_array = - StencilManager::instance() - .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, - list); + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_edges_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_edges_valid(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); - REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); - REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); - REQUIRE(symmetry_stencils_are_valid(stencil_array, mesh)); + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_faces_valid(stencil_array, mesh, 1)); + REQUIRE(not(are_stencil_by_nodes_valid(connectivity, stencil_array, 1))); + REQUIRE(are_stencil_by_faces_valid(connectivity, stencil_array, 1)); + } + } } - SECTION("hybrid") + SECTION("3D") { - const auto& mesh = *MeshDataBaseForTests::get().unordered1DMesh()->get<Mesh<1>>(); + 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() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_nodes_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_edges_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_edges_valid(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_faces_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_faces_valid(connectivity, stencil_array, 1)); + } + } - const Connectivity<1>& connectivity = mesh.connectivity(); - auto stencil = - StencilManager::instance() - .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, - list); + SECTION("hybrid") + { + const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); - REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 2); - REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity)); - REQUIRE(symmetry_stencils_are_valid(stencil, mesh)); + const Connectivity<3>& connectivity = mesh.connectivity(); + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_nodes_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_edges_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_edges_valid(connectivity, stencil_array, 1)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_faces_valid(stencil_array, mesh, 1)); + REQUIRE(are_stencil_by_faces_valid(connectivity, stencil_array, 1)); + } + } } } - SECTION("2D") + SECTION("2 layers") { - 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")); + NbGhostLayersTester nb_ghost_layers_tester(2); - SECTION("cartesian") + SECTION("1D") { - const auto& mesh = *MeshDataBaseForTests::get().cartesian2DMesh()->get<Mesh<2>>(); + StencilManager::BoundaryDescriptorList list; + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMIN")); + list.emplace_back(std::make_shared<NamedBoundaryDescriptor>("XMAX")); - const Connectivity<2>& connectivity = mesh.connectivity(); + 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>>(); - auto stencil_array = - StencilManager::instance() - .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, - list); + const Connectivity<1>& connectivity = mesh.connectivity(); - REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); - REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); - REQUIRE(symmetry_stencils_are_valid(stencil_array, mesh)); - } + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list); - SECTION("hybrid") - { - const auto& mesh = *MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>(); + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_nodes_valid(stencil_array, mesh, 2)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 2)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_edges_valid(stencil_array, mesh, 2)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 2)); + } - const Connectivity<2>& connectivity = mesh.connectivity(); - auto stencil = - StencilManager::instance() - .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, - list); + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list); - REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 4); - REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity)); - REQUIRE(symmetry_stencils_are_valid(stencil, mesh)); + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 2); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_faces_valid(stencil_array, mesh, 2)); + REQUIRE(are_stencil_by_nodes_valid(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") + SECTION("2D") { - const auto& mesh = *MeshDataBaseForTests::get().cartesian3DMesh()->get<Mesh<3>>(); + 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() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_nodes_valid(stencil_array, mesh, 2)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 2)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list); - const Connectivity<3>& connectivity = mesh.connectivity(); - auto stencil = - StencilManager::instance() - .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, - list); + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_edges_valid(stencil_array, mesh, 2)); + REQUIRE(are_stencil_by_edges_valid(connectivity, stencil_array, 2)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list); - REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 6); - REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity)); - REQUIRE(symmetry_stencils_are_valid(stencil, mesh)); + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 4); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_faces_valid(stencil_array, mesh, 2)); + REQUIRE(not(are_stencil_by_nodes_valid(connectivity, stencil_array, 2))); + REQUIRE(are_stencil_by_faces_valid(connectivity, stencil_array, 2)); + } + } } - SECTION("hybrid") + SECTION("3D") { - const auto& mesh = *MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); + 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")); - const Connectivity<3>& connectivity = mesh.connectivity(); - auto stencil = - StencilManager::instance() - .getCellToCellStencilArray(connectivity, StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, - list); + 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() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_nodes}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_nodes_valid(stencil_array, mesh, 2)); + REQUIRE(are_stencil_by_nodes_valid(connectivity, stencil_array, 2)); + } - REQUIRE(stencil.symmetryBoundaryStencilArrayList().size() == 6); - REQUIRE(check_ghost_cells_have_empty_stencils(stencil, connectivity)); - REQUIRE(symmetry_stencils_are_valid(stencil, mesh)); + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_edges}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_edges_valid(stencil_array, mesh, 2)); + REQUIRE(are_stencil_by_edges_valid(connectivity, stencil_array, 2)); + } + + { + auto stencil_array = + StencilManager::instance() + .getCellToCellStencilArray(connectivity, + StencilDescriptor{2, StencilDescriptor::ConnectionType::by_faces}, list); + + REQUIRE(stencil_array.symmetryBoundaryStencilArrayList().size() == 6); + REQUIRE(check_ghost_cells_have_empty_stencils(stencil_array, connectivity)); + REQUIRE(are_symmetry_stencils_by_faces_valid(stencil_array, mesh, 2)); + REQUIRE(are_stencil_by_faces_valid(connectivity, stencil_array, 2)); + } + } } } }