#ifndef MESH_NODE_BOUNDARY_HPP #define MESH_NODE_BOUNDARY_HPP #include <Kokkos_Vector.hpp> #include <algebra/TinyVector.hpp> #include <mesh/ConnectivityMatrix.hpp> #include <mesh/IConnectivity.hpp> #include <mesh/ItemValue.hpp> #include <mesh/RefItemList.hpp> #include <utils/Array.hpp> #include <utils/Exceptions.hpp> #include <utils/Messenger.hpp> template <size_t Dimension> class MeshNodeBoundary // clazy:exclude=copyable-polymorphic { protected: Array<const NodeId> m_node_list; std::string m_boundary_name; public: MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default; MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default; const Array<const NodeId>& nodeList() const { return m_node_list; } template <typename MeshType> MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list) : m_boundary_name(ref_face_list.refId().tagName()) { static_assert(Dimension == MeshType::Dimension); const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); const Array<const FaceId>& face_list = ref_face_list.list(); parallel_for( face_list.size(), PUGS_LAMBDA(int l) { const auto& face_cells = face_to_cell_matrix[face_list[l]]; if (face_cells.size() > 1) { 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"; throw NormalError(ost.str()); } }); Kokkos::vector<unsigned int> node_ids; // not enough but should reduce significantly the number of resizing node_ids.reserve(Dimension * face_list.size()); const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix(); for (size_t l = 0; l < face_list.size(); ++l) { const FaceId face_number = face_list[l]; const auto& face_nodes = face_to_node_matrix[face_number]; for (size_t r = 0; r < face_nodes.size(); ++r) { node_ids.push_back(face_nodes[r]); } } std::sort(node_ids.begin(), node_ids.end()); auto last = std::unique(node_ids.begin(), node_ids.end()); node_ids.resize(std::distance(node_ids.begin(), last)); Array<NodeId> node_list(node_ids.size()); parallel_for( node_ids.size(), PUGS_LAMBDA(int r) { node_list[r] = node_ids[r]; }); m_node_list = node_list; } template <typename MeshType> MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list) : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName()) { static_assert(Dimension == MeshType::Dimension); } MeshNodeBoundary() = default; virtual ~MeshNodeBoundary() = default; protected: MeshNodeBoundary(const MeshNodeBoundary&) = default; MeshNodeBoundary(MeshNodeBoundary&&) = default; }; template <size_t Dimension> class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension> // clazy:exclude=copyable-polymorphic { public: using Rd = TinyVector<Dimension, double>; private: const Rd m_outgoing_normal; template <typename MeshType> PUGS_INLINE Rd _getNormal(const MeshType& mesh); template <typename MeshType> PUGS_INLINE void _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal, const TinyVector<Dimension, double>& origin, const double length, const MeshType& mesh) const { static_assert(MeshType::Dimension == Dimension); const NodeValue<const Rd>& xr = mesh.xr(); parallel_for( this->m_node_list.size(), PUGS_LAMBDA(const size_t r) { const Rd& x = xr[this->m_node_list[r]]; if (dot(x - origin, normal) > 1E-13 * length) { std::ostringstream ost; ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset << ": boundary is not flat!"; throw NormalError(ost.str()); } }); } template <typename MeshType> PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh); public: const Rd& outgoingNormal() const { return m_outgoing_normal; } MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default; MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default; template <typename MeshType> MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list) : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh)) { ; } template <typename MeshType> MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list) : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh)) { ; } MeshFlatNodeBoundary() = default; MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default; MeshFlatNodeBoundary(MeshFlatNodeBoundary&&) = default; virtual ~MeshFlatNodeBoundary() = default; }; template <> template <typename MeshType> PUGS_INLINE TinyVector<1, double> MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh) { static_assert(MeshType::Dimension == 1); using R = TinyVector<1, double>; const size_t number_of_bc_nodes = [&]() { size_t number_of_bc_nodes = 0; auto node_is_owned = mesh.connectivity().nodeIsOwned(); for (size_t i_node = 0; i_node < m_node_list.size(); ++i_node) { number_of_bc_nodes += (node_is_owned[m_node_list[i_node]]); } return parallel::allReduceMax(number_of_bc_nodes); }(); if (number_of_bc_nodes != 1) { std::ostringstream ost; ost << "invalid boundary " << rang::fgB::yellow << m_boundary_name << rang::style::reset << ": node boundaries in 1D require to have exactly 1 node"; throw NormalError(ost.str()); } return R{1}; } template <> template <typename MeshType> PUGS_INLINE TinyVector<2, double> MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh) { static_assert(MeshType::Dimension == 2); using R2 = TinyVector<2, double>; const NodeValue<const R2>& xr = mesh.xr(); R2 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max()); R2 xmax(-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()); auto update_xmin = [](const R2& x, R2& xmin) { if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) { xmin = x; } }; auto update_xmax = [](const R2& x, R2& xmax) { if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) { xmax = x; } }; for (size_t r = 0; r < m_node_list.size(); ++r) { const R2& x = xr[m_node_list[r]]; update_xmin(x, xmin); update_xmax(x, xmax); } if (parallel::size() > 1) { Array<R2> xmin_array = parallel::allGather(xmin); Array<R2> xmax_array = parallel::allGather(xmax); for (size_t i = 0; i < xmin_array.size(); ++i) { update_xmin(xmin_array[i], xmin); } for (size_t i = 0; i < xmax_array.size(); ++i) { update_xmax(xmax_array[i], xmax); } } if (xmin == xmax) { std::ostringstream ost; ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset << ": unable to compute normal"; throw NormalError(ost.str()); } R2 dx = xmax - xmin; dx *= 1. / l2Norm(dx); const R2 normal(-dx[1], dx[0]); this->_checkBoundaryIsFlat(normal, 0.5 * (xmin + xmax), l2Norm(xmax - xmin), mesh); return normal; } template <> template <typename MeshType> PUGS_INLINE TinyVector<3, double> MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh) { static_assert(MeshType::Dimension == 3); using R3 = TinyVector<3, double>; auto update_xmin = [](const R3& x, R3& xmin) { // XMIN: X.xmin X.ymax X.zmax if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] > xmin[1])) or ((x[0] == xmin[0]) and (x[1] == xmin[1]) and (x[2] > xmin[2]))) { xmin = x; } }; auto update_xmax = [](const R3& x, R3& xmax) { // XMAX: X.xmax X.ymin X.zmin if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] < xmax[1])) or ((x[0] == xmax[0]) and (x[1] == xmax[1]) and (x[2] < xmax[2]))) { xmax = x; } }; auto update_ymin = [](const R3& x, R3& ymin) { // YMIN: X.ymin X.zmax X.xmin if ((x[1] < ymin[1]) or ((x[1] == ymin[1]) and (x[2] > ymin[2])) or ((x[1] == ymin[1]) and (x[2] == ymin[2]) and (x[0] < ymin[0]))) { ymin = x; } }; 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]))) { 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]))) { 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]))) { zmax = x; } }; R3 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()); R3 ymin = xmin; R3 zmin = xmin; R3 xmax = -xmin; R3 ymax = xmax; R3 zmax = xmax; 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_zmax(x, zmax); } if (parallel::size() > 0) { 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> zmax_array = parallel::allGather(zmax); for (size_t i = 0; i < xmin_array.size(); ++i) { update_xmin(xmin_array[i], xmin); } for (size_t i = 0; i < ymin_array.size(); ++i) { update_ymin(ymin_array[i], ymin); } for (size_t i = 0; i < zmin_array.size(); ++i) { update_zmin(zmin_array[i], zmin); } for (size_t i = 0; i < xmax_array.size(); ++i) { update_xmax(xmax_array[i], xmax); } for (size_t i = 0; i < ymax_array.size(); ++i) { update_ymax(ymax_array[i], ymax); } for (size_t i = 0; i < zmax_array.size(); ++i) { update_zmax(zmax_array[i], zmax); } } const R3 u = xmax - xmin; const R3 v = ymax - ymin; const R3 w = zmax - zmin; const R3 uv = crossProduct(u, v); const double uv_l2 = dot(uv, uv); R3 normal = uv; double normal_l2 = uv_l2; const R3 uw = crossProduct(u, w); const double uw_l2 = dot(uw, uw); if (uw_l2 > uv_l2) { normal = uw; normal_l2 = uw_l2; } const R3 vw = crossProduct(v, w); const double vw_l2 = dot(vw, vw); if (vw_l2 > normal_l2) { normal = vw; normal_l2 = vw_l2; } if (normal_l2 == 0) { std::ostringstream ost; ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset << ": unable to compute normal"; throw NormalError(ost.str()); } const double length = sqrt(normal_l2); normal *= 1. / length; this->_checkBoundaryIsFlat(normal, 1. / 6. * (xmin + xmax + ymin + ymax + zmin + zmax), length, mesh); return normal; } template <> template <typename MeshType> PUGS_INLINE TinyVector<1, double> MeshFlatNodeBoundary<1>::_getOutgoingNormal(const MeshType& mesh) { static_assert(MeshType::Dimension == 1); using R = TinyVector<1, double>; const R normal = this->_getNormal(mesh); double max_height = 0; if (m_node_list.size() > 0) { const NodeValue<const R>& xr = mesh.xr(); const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); const NodeId r0 = m_node_list[0]; const CellId j0 = node_to_cell_matrix[r0][0]; const auto& j0_nodes = cell_to_node_matrix[j0]; for (size_t r = 0; r < j0_nodes.size(); ++r) { const double height = dot(xr[j0_nodes[r]] - xr[r0], normal); if (std::abs(height) > std::abs(max_height)) { max_height = height; } } } Array<double> max_height_array = parallel::allGather(max_height); for (size_t i = 0; i < max_height_array.size(); ++i) { const double height = max_height_array[i]; if (std::abs(height) > std::abs(max_height)) { max_height = height; } } if (max_height > 0) { return -normal; } else { return normal; } } template <> template <typename MeshType> PUGS_INLINE TinyVector<2, double> MeshFlatNodeBoundary<2>::_getOutgoingNormal(const MeshType& mesh) { static_assert(MeshType::Dimension == 2); using R2 = TinyVector<2, double>; const R2 normal = this->_getNormal(mesh); double max_height = 0; if (m_node_list.size() > 0) { const NodeValue<const R2>& xr = mesh.xr(); const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); const NodeId r0 = m_node_list[0]; const CellId j0 = node_to_cell_matrix[r0][0]; const auto& j0_nodes = cell_to_node_matrix[j0]; for (size_t r = 0; r < j0_nodes.size(); ++r) { const double height = dot(xr[j0_nodes[r]] - xr[r0], normal); if (std::abs(height) > std::abs(max_height)) { max_height = height; } } } Array<double> max_height_array = parallel::allGather(max_height); for (size_t i = 0; i < max_height_array.size(); ++i) { const double height = max_height_array[i]; if (std::abs(height) > std::abs(max_height)) { max_height = height; } } if (max_height > 0) { return -normal; } else { return normal; } } template <> template <typename MeshType> PUGS_INLINE TinyVector<3, double> MeshFlatNodeBoundary<3>::_getOutgoingNormal(const MeshType& mesh) { static_assert(MeshType::Dimension == 3); using R3 = TinyVector<3, double>; const R3 normal = this->_getNormal(mesh); double max_height = 0; if (m_node_list.size() > 0) { const NodeValue<const R3>& xr = mesh.xr(); const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); const NodeId r0 = m_node_list[0]; const CellId j0 = node_to_cell_matrix[r0][0]; const auto& j0_nodes = cell_to_node_matrix[j0]; for (size_t r = 0; r < j0_nodes.size(); ++r) { const double height = dot(xr[j0_nodes[r]] - xr[r0], normal); if (std::abs(height) > std::abs(max_height)) { max_height = height; } } } Array<double> max_height_array = parallel::allGather(max_height); for (size_t i = 0; i < max_height_array.size(); ++i) { const double height = max_height_array[i]; if (std::abs(height) > std::abs(max_height)) { max_height = height; } } if (max_height > 0) { return -normal; } else { return normal; } } #endif // MESH_NODE_BOUNDARY_HPP