diff --git a/src/mesh/StencilBuilder.cpp b/src/mesh/StencilBuilder.cpp index 1880c049371f7e04303ae4ce82207c96c4d18c10..1fa54e7023246bb91bf5ca0d63197e2a14e6c8e7 100644 --- a/src/mesh/StencilBuilder.cpp +++ b/src/mesh/StencilBuilder.cpp @@ -7,85 +7,84 @@ #include <set> -template <typename ConnectivityType> +template <ItemType connecting_item_type, typename ConnectivityType> auto -StencilBuilder::_buildSymmetryNodeList(const ConnectivityType& connectivity, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const +StencilBuilder::_buildSymmetryConnectingItemList(const ConnectivityType& connectivity, + const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const { - static_assert(ConnectivityType::Dimension > 1); - NodeArray<bool> symmetry_node_list{connectivity, symmetry_boundary_descriptor_list.size()}; - symmetry_node_list.fill(false); - - auto face_to_node_matrix = connectivity.faceToNodeMatrix(); - 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_node_list = 0; i_ref_node_list < connectivity.template numberOfRefItemList<ItemType::face>(); - ++i_ref_node_list) { - const auto& ref_face_list = connectivity.template refItemList<ItemType::face>(i_ref_node_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 node_list = face_to_node_matrix[face_id]; - for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - const NodeId node_id = node_list[i_node]; - - symmetry_node_list[node_id][i_symmetry_boundary] = true; + static_assert(connecting_item_type != ItemType::cell, "cells cannot be used to define symmetry boundaries"); + + using ConnectingItemId = ItemIdT<connecting_item_type>; + ItemArray<bool, connecting_item_type> symmetry_connecting_item_list{connectivity, + 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; } - 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()); } } - ++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_node_list; -} -template <> -auto -StencilBuilder::_buildSymmetryNodeList(const Connectivity<1>& connectivity, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const -{ - NodeArray<bool> symmetry_node_list{connectivity, symmetry_boundary_descriptor_list.size()}; - symmetry_node_list.fill(false); - - 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_node_list = 0; i_ref_node_list < connectivity.template numberOfRefItemList<ItemType::node>(); - ++i_ref_node_list) { - const auto& ref_node_list = connectivity.template refItemList<ItemType::node>(i_ref_node_list); - if (ref_node_list.refId() == boundary_descriptor) { - found = true; - auto node_list = ref_node_list.list(); - for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - const NodeId node_id = node_list[i_node]; - - symmetry_node_list[node_id][i_symmetry_boundary] = true; + return symmetry_connecting_item_list; + } else { + 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_connecting_item_list = 0; + i_ref_connecting_item_list < connectivity.template numberOfRefItemList<connecting_item_type>(); + ++i_ref_connecting_item_list) { + const auto& ref_connecting_item_list = + connectivity.template refItemList<connecting_item_type>(i_ref_connecting_item_list); + if (ref_connecting_item_list.refId() == boundary_descriptor) { + found = true; + auto connecting_item_list = ref_connecting_item_list.list(); + 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; } - 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()); } } - ++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_node_list; + return symmetry_connecting_item_list; + } } template <typename ConnectivityType> @@ -161,8 +160,8 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar } } if (not found) { - int node_cell_number = cell_number[node_cell_id]; - size_t i_index = row_map[cell_id]; + 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])]) { @@ -173,7 +172,7 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar } for (size_t i_destination = max_index[cell_id]; i_destination > i_index; --i_destination) { - size_t i_source = i_destination - 1; + const size_t i_source = i_destination - 1; column_indices[i_destination] = column_indices[i_source]; } @@ -189,24 +188,26 @@ StencilBuilder::_getColumnIndices(const ConnectivityType& connectivity, const Ar return column_indices; } -template <typename ConnectivityType> -CellToCellStencilArray -StencilBuilder::_buildC2C(const ConnectivityType& connectivity, - size_t number_of_layers, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const +template <ItemType item_type, ItemType connecting_item_type, typename ConnectivityType> +StencilArray<item_type, item_type> +StencilBuilder::_build_for_same_source_and_target(const ConnectivityType& connectivity, + const size_t& number_of_layers, + const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const { - if ((parallel::size() > 1) and (number_of_layers > GlobalVariableManager::instance().getNumberOfGhostLayers())) { - std::ostringstream error_msg; - error_msg << "Stencil builder requires" << rang::fgB::yellow << number_of_layers << 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()); + 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>(); + + auto item_is_owned = connectivity.template isOwned<item_type>(); + auto item_number = connectivity.template number<item_type>(); + + using ItemId = ItemIdT<item_type>; + using ConnectingItemId = ItemIdT<connecting_item_type>; if (symmetry_boundary_descriptor_list.size() == 0) { if (number_of_layers == 1) { @@ -215,46 +216,46 @@ StencilBuilder::_buildC2C(const ConnectivityType& connectivity, return {ConnectivityMatrix{row_map, column_indices}, {}}; } else { - auto cell_to_node_matrix = connectivity.cellToNodeMatrix(); - auto node_to_cell_matrix = connectivity.nodeToCellMatrix(); - - auto cell_is_owned = connectivity.cellIsOwned(); - auto cell_number = connectivity.cellNumber(); - - Array<uint32_t> row_map{connectivity.numberOfCells() + 1}; + Array<uint32_t> row_map{connectivity.template numberOf<item_type>() + 1}; row_map[0] = 0; - std::vector<CellId> column_indices_vector; + std::vector<ItemId> column_indices_vector; - for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) { - if (cell_is_owned[cell_id]) { - std::set<CellId, std::function<bool(CellId, CellId)>> cell_set( - [=](CellId cell_0, CellId cell_1) { return cell_number[cell_0] < cell_number[cell_1]; }); + 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_node_1 = 0; i_node_1 < cell_to_node_matrix[cell_id].size(); ++i_node_1) { - const NodeId layer_1_node_id = cell_to_node_matrix[cell_id][i_node_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_cell_1 = 0; i_cell_1 < node_to_cell_matrix[layer_1_node_id].size(); ++i_cell_1) { - CellId cell_1_id = node_to_cell_matrix[layer_1_node_id][i_cell_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_node_2 = 0; i_node_2 < cell_to_node_matrix[cell_1_id].size(); ++i_node_2) { - const NodeId layer_2_node_id = cell_to_node_matrix[cell_1_id][i_node_2]; + 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_cell_2 = 0; i_cell_2 < node_to_cell_matrix[layer_2_node_id].size(); ++i_cell_2) { - CellId cell_2_id = node_to_cell_matrix[layer_2_node_id][i_cell_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 (cell_2_id != cell_id) { - cell_set.insert(cell_2_id); + if (item_2_id != item_id) { + item_set.insert(item_2_id); } } } } } - for (auto stencil_cell_id : cell_set) { - column_indices_vector.push_back(stencil_cell_id); + for (auto stencil_item_id : item_set) { + column_indices_vector.push_back(stencil_item_id); } - row_map[cell_id + 1] = row_map[cell_id] + cell_set.size(); + row_map[item_id + 1] = row_map[item_id] + item_set.size(); } } @@ -273,99 +274,98 @@ StencilBuilder::_buildC2C(const ConnectivityType& connectivity, return {primal_stencil, {}}; } } else { - NodeArray<bool> symmetry_node_list = this->_buildSymmetryNodeList(connectivity, symmetry_boundary_descriptor_list); - - auto cell_to_node_matrix = connectivity.cellToNodeMatrix(); - auto node_to_cell_matrix = connectivity.nodeToCellMatrix(); + ItemArray<bool, connecting_item_type> symmetry_item_list = + this->_buildSymmetryConnectingItemList<connecting_item_type>(connectivity, symmetry_boundary_descriptor_list); - auto cell_is_owned = connectivity.cellIsOwned(); - auto cell_number = connectivity.cellNumber(); - - Array<uint32_t> row_map{connectivity.numberOfCells() + 1}; + 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.numberOfCells() + 1}; + symmetry_row_map = Array<uint32_t>{connectivity.template numberOf<item_type>() + 1}; symmetry_row_map[0] = 0; } std::vector<uint32_t> column_indices_vector; std::vector<std::vector<uint32_t>> symmetry_column_indices_vector(symmetry_boundary_descriptor_list.size()); - for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) { - std::set<CellId> cell_set; - std::vector<std::set<CellId>> by_boundary_symmetry_cell(symmetry_boundary_descriptor_list.size()); - - if (cell_is_owned[cell_id]) { - auto cell_node_list = cell_to_node_matrix[cell_id]; - for (size_t i_cell_node = 0; i_cell_node < cell_node_list.size(); ++i_cell_node) { - const NodeId cell_node_id = cell_node_list[i_cell_node]; - auto node_cell_list = node_to_cell_matrix[cell_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]; - if (cell_id != node_cell_id) { - cell_set.insert(node_cell_id); + 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()); + + 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(); + ++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); } } } { - std::vector<CellId> cell_vector; - for (auto&& set_cell_id : cell_set) { - cell_vector.push_back(set_cell_id); + std::vector<ItemId> item_vector; + for (auto&& set_item_id : item_set) { + item_vector.push_back(set_item_id); } - std::sort(cell_vector.begin(), cell_vector.end(), - [&cell_number](const CellId& cell0_id, const CellId& cell1_id) { - return cell_number[cell0_id] < cell_number[cell1_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_cell_id : cell_vector) { - column_indices_vector.push_back(vector_cell_id); + for (auto&& vector_item_id : item_vector) { + column_indices_vector.push_back(vector_item_id); } } for (size_t i = 0; i < symmetry_boundary_descriptor_list.size(); ++i) { - std::set<CellId> symmetry_cell_set; - for (size_t i_cell_node = 0; i_cell_node < cell_node_list.size(); ++i_cell_node) { - const NodeId cell_node_id = cell_node_list[i_cell_node]; - if (symmetry_node_list[cell_node_id][i]) { - auto node_cell_list = node_to_cell_matrix[cell_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]; - symmetry_cell_set.insert(node_cell_id); + 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); } } } - by_boundary_symmetry_cell[i] = symmetry_cell_set; + by_boundary_symmetry_item[i] = symmetry_item_set; - std::vector<CellId> cell_vector; - for (auto&& set_cell_id : symmetry_cell_set) { - cell_vector.push_back(set_cell_id); + std::vector<ItemId> item_vector; + for (auto&& set_item_id : symmetry_item_set) { + item_vector.push_back(set_item_id); } - std::sort(cell_vector.begin(), cell_vector.end(), - [&cell_number](const CellId& cell0_id, const CellId& cell1_id) { - return cell_number[cell0_id] < cell_number[cell1_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_cell_id : cell_vector) { - symmetry_column_indices_vector[i].push_back(vector_cell_id); + for (auto&& vector_item_id : item_vector) { + symmetry_column_indices_vector[i].push_back(vector_item_id); } } } - row_map[cell_id + 1] = row_map[cell_id] + cell_set.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][cell_id + 1] = symmetry_row_map_list[i][cell_id] + by_boundary_symmetry_cell[i].size(); + symmetry_row_map_list[i][item_id + 1] = symmetry_row_map_list[i][item_id] + by_boundary_symmetry_item[i].size(); } } ConnectivityMatrix primal_stencil{row_map, convert_to_array(column_indices_vector)}; - CellToCellStencilArray::BoundaryDescriptorStencilArrayList symmetry_boundary_stencil_list; + 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( - CellToCellStencilArray:: + typename StencilArray<item_type, item_type>:: BoundaryDescriptorStencilArray{p_boundary_descriptor, ConnectivityMatrix{symmetry_row_map_list[i], convert_to_array(symmetry_column_indices_vector[i])}}); @@ -377,63 +377,74 @@ StencilBuilder::_buildC2C(const ConnectivityType& connectivity, } } -template <> -CellToCellStencilArray -StencilBuilder::_buildC2C(const Connectivity<1>& connectivity, - const StencilDescriptor& stencil_descriptor, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const +template <ItemType source_item_type, + ItemType connecting_item_type, + ItemType target_item_type, + typename ConnectivityType> +StencilArray<source_item_type, target_item_type> +StencilBuilder::_build_for_different_source_and_target( + const ConnectivityType& connectivity, + const size_t& number_of_layers, + const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const { - if (stencil_descriptor.connectionType() != StencilDescriptor::ConnectionType::by_cells) { - return StencilBuilder::_buildC2C(connectivity, stencil_descriptor.numberOfLayers(), - symmetry_boundary_descriptor_list); - } else { - throw UnexpectedError("Cannot build cell to cell stencil using cell connectivity"); - } + static_assert(source_item_type != target_item_type); + throw NotImplementedError("different source target"); } -template <> -CellToCellStencilArray -StencilBuilder::_buildC2C(const Connectivity<2>& connectivity, - const StencilDescriptor& stencil_descriptor, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const +template <ItemType source_item_type, + ItemType connecting_item_type, + ItemType target_item_type, + typename ConnectivityType> +StencilArray<source_item_type, target_item_type> +StencilBuilder::_build(const ConnectivityType& connectivity, + const size_t& number_of_layers, + const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const { - switch (stencil_descriptor.connectionType()) { - case StencilDescriptor::ConnectionType::by_nodes: { - return StencilBuilder::_buildC2C(connectivity, stencil_descriptor.numberOfLayers(), - symmetry_boundary_descriptor_list); - } - case StencilDescriptor::ConnectionType::by_edges: - case StencilDescriptor::ConnectionType::by_faces: { - throw NotImplementedError("non by_node stencil"); - } - case StencilDescriptor::ConnectionType::by_cells: { - throw UnexpectedError("Cannot build cell to cell stencil using cell connectivity"); - } - // LCOV_EXCL_START - default: { - throw UnexpectedError("invalid connection type"); - } - // LCOV_EXCL_STOP + if constexpr (connecting_item_type != target_item_type) { + if constexpr (source_item_type == target_item_type) { + return this + ->_build_for_same_source_and_target<source_item_type, connecting_item_type>(connectivity, number_of_layers, + symmetry_boundary_descriptor_list); + } else { + return this->_build_for_different_source_and_target<source_item_type, connecting_item_type, + target_item_type>(connectivity, number_of_layers, + symmetry_boundary_descriptor_list); + } + } else { + std::ostringstream error_msg; + error_msg << "cannot build stencil of " << rang::fgB::yellow << itemName(target_item_type) << rang::fg::reset + << " using " << rang::fgB::yellow << itemName(connecting_item_type) << rang::fg::reset + << " for connectivity"; + throw UnexpectedError(error_msg.str()); } } -template <> -CellToCellStencilArray -StencilBuilder::_buildC2C(const Connectivity<3>& connectivity, - const StencilDescriptor& stencil_descriptor, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const +template <ItemType source_item_type, ItemType target_item_type, typename ConnectivityType> +StencilArray<source_item_type, target_item_type> +StencilBuilder::_build(const ConnectivityType& connectivity, + const StencilDescriptor& stencil_descriptor, + const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const { switch (stencil_descriptor.connectionType()) { case StencilDescriptor::ConnectionType::by_nodes: { - return StencilBuilder::_buildC2C(connectivity, stencil_descriptor.numberOfLayers(), - symmetry_boundary_descriptor_list); + return this->_build<source_item_type, ItemType::node, target_item_type>(connectivity, + stencil_descriptor.numberOfLayers(), + symmetry_boundary_descriptor_list); + } + case StencilDescriptor::ConnectionType::by_edges: { + return this->_build<source_item_type, ItemType::edge, target_item_type>(connectivity, + stencil_descriptor.numberOfLayers(), + symmetry_boundary_descriptor_list); } - case StencilDescriptor::ConnectionType::by_edges: case StencilDescriptor::ConnectionType::by_faces: { - throw NotImplementedError("non by_node stencil"); + return this->_build<source_item_type, ItemType::face, target_item_type>(connectivity, + stencil_descriptor.numberOfLayers(), + symmetry_boundary_descriptor_list); } case StencilDescriptor::ConnectionType::by_cells: { - throw UnexpectedError("Cannot build cell to cell stencil using cell connectivity"); + return this->_build<source_item_type, ItemType::cell, target_item_type>(connectivity, + stencil_descriptor.numberOfLayers(), + symmetry_boundary_descriptor_list); } // LCOV_EXCL_START default: { @@ -448,18 +459,28 @@ StencilBuilder::buildC2C(const IConnectivity& connectivity, const StencilDescriptor& stencil_descriptor, const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const { + 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 StencilBuilder::_buildC2C(dynamic_cast<const Connectivity<1>&>(connectivity), stencil_descriptor, - symmetry_boundary_descriptor_list); + return this->_build<ItemType::cell, ItemType::cell>(dynamic_cast<const Connectivity<1>&>(connectivity), + stencil_descriptor, symmetry_boundary_descriptor_list); } case 2: { - return StencilBuilder::_buildC2C(dynamic_cast<const Connectivity<2>&>(connectivity), stencil_descriptor, - symmetry_boundary_descriptor_list); + return this->_build<ItemType::cell, ItemType::cell>(dynamic_cast<const Connectivity<2>&>(connectivity), + stencil_descriptor, symmetry_boundary_descriptor_list); } case 3: { - return StencilBuilder::_buildC2C(dynamic_cast<const Connectivity<3>&>(connectivity), stencil_descriptor, - symmetry_boundary_descriptor_list); + return this->_build<ItemType::cell, ItemType::cell>(dynamic_cast<const Connectivity<3>&>(connectivity), + stencil_descriptor, symmetry_boundary_descriptor_list); } default: { throw UnexpectedError("invalid connectivity dimension"); @@ -467,12 +488,6 @@ StencilBuilder::buildC2C(const IConnectivity& connectivity, } } -CellToFaceStencilArray -StencilBuilder::buildC2F(const IConnectivity&, const StencilDescriptor&, const BoundaryDescriptorList&) const -{ - throw NotImplementedError("cell to face stencil"); -} - NodeToCellStencilArray StencilBuilder::buildN2C(const IConnectivity&, const StencilDescriptor&, const BoundaryDescriptorList&) const { diff --git a/src/mesh/StencilBuilder.hpp b/src/mesh/StencilBuilder.hpp index 464a7fc2732a34f9365ea07ca904e6da54a53054..5bfd821ca49eabd0131f8895477d2058f3eaaff0 100644 --- a/src/mesh/StencilBuilder.hpp +++ b/src/mesh/StencilBuilder.hpp @@ -23,33 +23,39 @@ class StencilBuilder Array<const uint32_t> _getColumnIndices(const ConnectivityType& connectivity, const Array<const uint32_t>& row_map) const; - template <typename ConnectivityType> - auto _buildSymmetryNodeList(const ConnectivityType& connectivity, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; + template <ItemType connecting_item, typename ConnectivityType> + auto _buildSymmetryConnectingItemList(const ConnectivityType& connectivity, + const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; - template <typename ConnectivityType, - ItemType source_item_type, - ItemType connecting_item_type, - ItemType target_item_type> - StencilArray<source_item_type, target_item_type> _build( + template <ItemType item_type, ItemType connecting_item_type, typename ConnectivityType> + StencilArray<item_type, item_type> _build_for_same_source_and_target( const ConnectivityType& connectivity, - const StencilDescriptor& stencil_descriptor, + const size_t& number_of_layers, const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; - template <typename ConnectivityType> - CellToCellStencilArray _buildC2C(const ConnectivityType& connectivity, - const StencilDescriptor& stencil_descriptor, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; + template <ItemType source_item_type, + ItemType connecting_item_type, + ItemType target_item_type, + typename ConnectivityType> + StencilArray<source_item_type, target_item_type> _build_for_different_source_and_target( + const ConnectivityType& connectivity, + const size_t& number_of_layers, + const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; - template <typename ConnectivityType> - CellToCellStencilArray _buildC2C(const ConnectivityType& connectivity, - size_t number_of_layers, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; + template <ItemType source_item_type, + ItemType connecting_item_type, + ItemType target_item_type, + typename ConnectivityType> + StencilArray<source_item_type, target_item_type> _build( + const ConnectivityType& connectivity, + const size_t& number_of_layers, + const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; - template <typename ConnectivityType> - CellToFaceStencilArray _buildC2F(const ConnectivityType& connectivity, - size_t number_of_layers, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; + template <ItemType source_item_type, ItemType target_item_type, typename ConnectivityType> + StencilArray<source_item_type, target_item_type> _build( + const ConnectivityType& connectivity, + const StencilDescriptor& stencil_descriptor, + const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; template <typename ConnectivityType> NodeToCellStencilArray _buildN2C(const ConnectivityType& connectivity, @@ -60,10 +66,6 @@ class StencilBuilder const StencilDescriptor& stencil_descriptor, const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; - CellToFaceStencilArray buildC2F(const IConnectivity& connectivity, - const StencilDescriptor& stencil_descriptor, - const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; - NodeToCellStencilArray buildN2C(const IConnectivity& connectivity, const StencilDescriptor& stencil_descriptor, const BoundaryDescriptorList& symmetry_boundary_descriptor_list) const; @@ -77,8 +79,6 @@ class StencilBuilder { if constexpr ((source_item_type == ItemType::cell) and (target_item_type == ItemType::cell)) { return buildC2C(connectivity, stencil_descriptor, symmetry_boundary_descriptor_list); - } else if constexpr ((source_item_type == ItemType::cell) and (target_item_type == ItemType::face)) { - return buildC2F(connectivity, stencil_descriptor, symmetry_boundary_descriptor_list); } else if constexpr ((source_item_type == ItemType::node) and (target_item_type == ItemType::cell)) { return buildN2C(connectivity, stencil_descriptor, symmetry_boundary_descriptor_list); } else {