diff --git a/src/mesh/MeshNodeBoundary.cpp b/src/mesh/MeshNodeBoundary.cpp index 56ee047545d51e20f11e7b0421dae839615bf2cc..6c18637a09baddb9241b95033e729ad644259f5c 100644 --- a/src/mesh/MeshNodeBoundary.cpp +++ b/src/mesh/MeshNodeBoundary.cpp @@ -85,23 +85,23 @@ MeshNodeBoundary<3>::_getBounds(const Mesh<Connectivity<3>>& mesh) const auto update_ymax = [](const R3& x, R3& ymax) { // YMAX: X.ymax X.zmin X.xmax if ((x[1] > ymax[1]) or ((x[1] == ymax[1]) and (x[2] < ymax[2])) or - ((x[1] == ymax[1]) and (x[2] == ymax[1]) and (x[0] > ymax[0]))) { + ((x[1] == ymax[1]) and (x[2] == ymax[2]) and (x[0] > ymax[0]))) { ymax = x; } }; auto update_zmin = [](const R3& x, R3& zmin) { // ZMIN: X.zmin X.xmin X.ymin - if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[2] < zmin[2])) or - ((x[1] == zmin[1]) and (x[2] == zmin[1]) and (x[0] < zmin[0]))) { + if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[0] < zmin[0])) or + ((x[2] == zmin[2]) and (x[0] == zmin[0]) and (x[1] < zmin[1]))) { zmin = x; } }; auto update_zmax = [](const R3& x, R3& zmax) { // ZMAX: X.zmax X.xmax X.ymax - if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[2] > zmax[2])) or - ((x[1] == zmax[1]) and (x[2] == zmax[1]) and (x[0] > zmax[0]))) { + if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[0] > zmax[0])) or + ((x[2] == zmax[2]) and (x[0] == zmax[0]) and (x[1] > zmax[1]))) { zmax = x; } }; @@ -115,31 +115,34 @@ MeshNodeBoundary<3>::_getBounds(const Mesh<Connectivity<3>>& mesh) const R3& zmax = bounds[5]; xmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()}; - ymin = xmin; - zmin = xmin; + ymin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()}; + zmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()}; - xmax = -xmin; - ymax = xmax; - zmax = xmax; + xmax = + -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()}; + ymax = + -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()}; + zmax = + -R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()}; const NodeValue<const R3>& xr = mesh.xr(); for (size_t r = 0; r < m_node_list.size(); ++r) { const R3& x = xr[m_node_list[r]]; update_xmin(x, xmin); - update_xmax(x, xmax); update_ymin(x, ymin); - update_ymax(x, ymax); update_zmin(x, zmin); + update_xmax(x, xmax); + update_ymax(x, ymax); update_zmax(x, zmax); } - if (parallel::size() > 0) { + if (parallel::size() > 1) { Array<const R3> xmin_array = parallel::allGather(xmin); - Array<const R3> xmax_array = parallel::allGather(xmax); Array<const R3> ymin_array = parallel::allGather(ymin); - Array<const R3> ymax_array = parallel::allGather(ymax); Array<const R3> zmin_array = parallel::allGather(zmin); + Array<const R3> xmax_array = parallel::allGather(xmax); + Array<const R3> ymax_array = parallel::allGather(ymax); Array<const R3> zmax_array = parallel::allGather(zmax); for (size_t i = 0; i < xmin_array.size(); ++i) { @@ -170,19 +173,8 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension> const RefFaceList& ref_face_list) : m_boundary_name(ref_face_list.refId().tagName()) { - const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); - const Array<const FaceId>& face_list = ref_face_list.list(); - - bool is_bad = false; - parallel_for(face_list.size(), [=, &is_bad](int l) { - const auto& face_cells = face_to_cell_matrix[face_list[l]]; - if (face_cells.size() > 1) { - is_bad = true; - } - }); - - if (parallel::allReduceOr(is_bad)) { + if (not ref_face_list.isBoundary()) { std::ostringstream ost; ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset << ": inner faces cannot be used to define mesh boundaries"; @@ -217,6 +209,11 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension> face_list.size(), PUGS_LAMBDA(int r) { node_list[r] = static_cast<FaceId::base_type>(face_list[r]); }); m_node_list = node_list; } + + // This is quite dirty but it allows a non negligible performance + // improvement + const_cast<Connectivity<Dimension>&>(mesh.connectivity()) + .addRefItemList(RefItemList<ItemType::node>(ref_face_list.refId(), m_node_list, ref_face_list.isBoundary())); } template <size_t Dimension> @@ -225,44 +222,7 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension> : m_boundary_name(ref_edge_list.refId().tagName()) { const Array<const EdgeId>& edge_list = ref_edge_list.list(); - const auto& edge_is_owned = mesh.connectivity().edgeIsOwned(); - - bool is_bad = false; - - if constexpr ((Dimension == 1) or (Dimension == 2)) { - const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); - parallel_for(edge_list.size(), [=, &is_bad](int l) { - const EdgeId edge_id = edge_list[l]; - if (edge_is_owned[edge_id]) { - if (edge_to_cell_matrix[edge_id].size() != 1) { - is_bad = true; - } - } - }); - } else { - static_assert(Dimension == 3); - const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix(); - const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); - - parallel_for(edge_list.size(), [=, &is_bad](int l) { - const EdgeId edge_id = edge_list[l]; - if (edge_is_owned[edge_id]) { - const auto& edge_faces = edge_to_face_matrix[edge_id]; - bool is_connected_to_boundary_face = false; - for (size_t i_edge_face = 0; i_edge_face < edge_faces.size(); ++i_edge_face) { - const FaceId edge_face_id = edge_faces[i_edge_face]; - if (face_to_cell_matrix[edge_face_id].size() == 1) { - is_connected_to_boundary_face = true; - } - } - if (not is_connected_to_boundary_face) { - is_bad = true; - } - } - }); - } - - if (parallel::allReduceOr(is_bad)) { + if (not ref_edge_list.isBoundary()) { std::ostringstream ost; ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset << ": inner edges cannot be used to define mesh boundaries"; @@ -296,12 +256,24 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension> edge_list.size(), PUGS_LAMBDA(int r) { node_list[r] = static_cast<EdgeId::base_type>(edge_list[r]); }); m_node_list = node_list; } + + // This is quite dirty but it allows a non negligible performance + // improvement + const_cast<Connectivity<Dimension>&>(mesh.connectivity()) + .addRefItemList(RefItemList<ItemType::node>(ref_edge_list.refId(), m_node_list, ref_edge_list.isBoundary())); } template <size_t Dimension> MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>>&, const RefNodeList& ref_node_list) : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName()) -{} +{ + if (not ref_node_list.isBoundary()) { + std::ostringstream ost; + ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset + << ": inner nodes cannot be used to define mesh boundaries"; + throw NormalError(ost.str()); + } +} template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefFaceList&); template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);