From d1862a0ad2bb83c0fadacd33a73b2918fca6454b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 19 Jul 2021 19:15:06 +0200
Subject: [PATCH] Reorganize Mesh*NodeBoundary classes
---
src/mesh/CMakeLists.txt | 3 +
src/mesh/MeshFlatNodeBoundary.cpp | 296 ++++++++++
src/mesh/MeshFlatNodeBoundary.hpp | 60 ++
src/mesh/MeshLineNodeBoundary.cpp | 143 +++++
src/mesh/MeshLineNodeBoundary.hpp | 65 +++
src/mesh/MeshNodeBoundary.cpp | 323 +++++++++++
src/mesh/MeshNodeBoundary.hpp | 885 +-----------------------------
src/mesh/MeshRandomizer.hpp | 2 +
src/scheme/AcousticSolver.cpp | 1 +
9 files changed, 903 insertions(+), 875 deletions(-)
create mode 100644 src/mesh/MeshFlatNodeBoundary.cpp
create mode 100644 src/mesh/MeshFlatNodeBoundary.hpp
create mode 100644 src/mesh/MeshLineNodeBoundary.cpp
create mode 100644 src/mesh/MeshLineNodeBoundary.hpp
create mode 100644 src/mesh/MeshNodeBoundary.cpp
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 6c0088831..2d429282e 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -17,6 +17,9 @@ add_library(
LogicalConnectivityBuilder.cpp
MeshBuilderBase.cpp
MeshDataManager.cpp
+ MeshFlatNodeBoundary.cpp
+ MeshLineNodeBoundary.cpp
+ MeshNodeBoundary.cpp
SynchronizerManager.cpp)
# Additional dependencies
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
new file mode 100644
index 000000000..08e0c800d
--- /dev/null
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -0,0 +1,296 @@
+#include <mesh/MeshFlatNodeBoundary.hpp>
+
+template <size_t Dimension>
+void
+MeshFlatNodeBoundary<Dimension>::_checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
+ const TinyVector<Dimension, double>& origin,
+ const double length,
+ const Mesh<Connectivity<Dimension>>& mesh) const
+{
+ const NodeValue<const Rd>& xr = mesh.xr();
+
+ bool is_bad = false;
+
+ parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
+ const Rd& x = xr[this->m_node_list[r]];
+ if (dot(x - origin, normal) > 1E-13 * length) {
+ is_bad = true;
+ }
+ });
+
+ if (parallel::allReduceOr(is_bad)) {
+ 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 <>
+TinyVector<1, double>
+MeshFlatNodeBoundary<1>::_getNormal(const Mesh<Connectivity<1>>& mesh)
+{
+ 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 <>
+TinyVector<2, double>
+MeshFlatNodeBoundary<2>::_getNormal(const Mesh<Connectivity<2>>& mesh)
+{
+ using R2 = TinyVector<2, double>;
+
+ std::array<R2, 2> bounds = this->_getBounds(mesh);
+
+ const R2& xmin = bounds[0];
+ const R2& xmax = bounds[1];
+
+ 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 <>
+TinyVector<3, double>
+MeshFlatNodeBoundary<3>::_getNormal(const Mesh<Connectivity<3>>& mesh)
+{
+ using R3 = TinyVector<3, double>;
+
+ std::array<R3, 6> bounds = this->_getBounds(mesh);
+
+ const R3& xmin = bounds[0];
+ const R3& ymin = bounds[1];
+ const R3& zmin = bounds[2];
+ const R3& xmax = bounds[3];
+ const R3& ymax = bounds[4];
+ const R3& zmax = bounds[5];
+
+ 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 <>
+TinyVector<1, double>
+MeshFlatNodeBoundary<1>::_getOutgoingNormal(const Mesh<Connectivity<1>>& mesh)
+{
+ 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 <>
+TinyVector<2, double>
+MeshFlatNodeBoundary<2>::_getOutgoingNormal(const Mesh<Connectivity<2>>& mesh)
+{
+ 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 <>
+TinyVector<3, double>
+MeshFlatNodeBoundary<3>::_getOutgoingNormal(const Mesh<Connectivity<3>>& mesh)
+{
+ 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;
+ }
+}
+
+template <size_t Dimension>
+MeshFlatNodeBoundary<Dimension>
+getMeshFlatNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+{
+ for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
+ ++i_ref_node_list) {
+ const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
+ const RefId& ref = ref_node_list.refId();
+ if (ref == boundary_descriptor) {
+ return MeshFlatNodeBoundary<Dimension>{mesh, ref_node_list};
+ }
+ }
+ for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
+ ++i_ref_face_list) {
+ const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == boundary_descriptor) {
+ return MeshFlatNodeBoundary<Dimension>{mesh, ref_face_list};
+ }
+ }
+
+ std::ostringstream ost;
+ ost << "cannot find surface with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
+
+ throw NormalError(ost.str());
+}
+
+template MeshFlatNodeBoundary<1> getMeshFlatNodeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshFlatNodeBoundary<2> getMeshFlatNodeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshFlatNodeBoundary<3> getMeshFlatNodeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshFlatNodeBoundary.hpp b/src/mesh/MeshFlatNodeBoundary.hpp
new file mode 100644
index 000000000..334cf03aa
--- /dev/null
+++ b/src/mesh/MeshFlatNodeBoundary.hpp
@@ -0,0 +1,60 @@
+#ifndef MESH_FLAT_NODE_BOUNDARY_HPP
+#define MESH_FLAT_NODE_BOUNDARY_HPP
+
+#include <mesh/MeshNodeBoundary.hpp>
+
+template <size_t Dimension>
+class MeshFlatNodeBoundary final : public MeshNodeBoundary<Dimension> // clazy:exclude=copyable-polymorphic
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const Rd m_outgoing_normal;
+
+ Rd _getNormal(const Mesh<Connectivity<Dimension>>& mesh);
+
+ void _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
+ const TinyVector<Dimension, double>& origin,
+ const double length,
+ const Mesh<Connectivity<Dimension>>& mesh) const;
+
+ Rd _getOutgoingNormal(const Mesh<Connectivity<Dimension>>& mesh);
+
+ public:
+ const Rd&
+ outgoingNormal() const
+ {
+ return m_outgoing_normal;
+ }
+
+ MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
+ MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
+
+ template <size_t MeshDimension>
+ friend MeshFlatNodeBoundary<MeshDimension> getMeshFlatNodeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
+ const IBoundaryDescriptor& boundary_descriptor);
+
+ private:
+ 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))
+ {}
+
+ public:
+ MeshFlatNodeBoundary() = default;
+ MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
+ MeshFlatNodeBoundary(MeshFlatNodeBoundary&&) = default;
+ ~MeshFlatNodeBoundary() = default;
+};
+
+template <size_t Dimension>
+MeshFlatNodeBoundary<Dimension> getMeshFlatNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+ const IBoundaryDescriptor& boundary_descriptor);
+
+#endif // MESH_FLAT_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshLineNodeBoundary.cpp b/src/mesh/MeshLineNodeBoundary.cpp
new file mode 100644
index 000000000..391a7aeca
--- /dev/null
+++ b/src/mesh/MeshLineNodeBoundary.cpp
@@ -0,0 +1,143 @@
+#include <mesh/MeshLineNodeBoundary.hpp>
+
+template <size_t Dimension>
+void
+MeshLineNodeBoundary<Dimension>::_checkBoundaryIsLine(const TinyVector<Dimension, double>& direction,
+ const TinyVector<Dimension, double>& origin,
+ const double length,
+ const Mesh<Connectivity<Dimension>>& mesh) const
+{
+ using Rdxd = TinyMatrix<Dimension>;
+
+ const NodeValue<const Rd>& xr = mesh.xr();
+
+ Rdxd P = Rdxd{identity} - tensorProduct(direction, direction);
+
+ bool is_bad = false;
+ parallel_for(this->m_node_list.size(), [=, &is_bad](int node_id) {
+ const Rd& x = xr[this->m_node_list[node_id]];
+ const Rd delta = x - origin;
+ if (dot(P * delta, direction) > 1E-13 * length) {
+ is_bad = true;
+ }
+ });
+
+ if (parallel::allReduceOr(is_bad)) {
+ std::ostringstream ost;
+ ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+ << ": boundary is not a line!";
+ throw NormalError(ost.str());
+ }
+}
+
+template <>
+template <>
+TinyVector<1>
+MeshLineNodeBoundary<1>::_getDirection(const Mesh<Connectivity<1>>&)
+{
+ throw UnexpectedError("MeshLineNodeBoundary makes no sense in dimension 1");
+}
+
+template <>
+template <>
+TinyVector<2>
+MeshLineNodeBoundary<2>::_getDirection(const Mesh<Connectivity<2>>& mesh)
+{
+ using R2 = TinyVector<2, double>;
+
+ std::array<R2, 2> bounds = this->_getBounds(mesh);
+
+ const R2& xmin = bounds[0];
+ const R2& xmax = bounds[1];
+
+ if (xmin == xmax) {
+ std::ostringstream ost;
+ ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+ << ": unable to compute direction";
+ throw NormalError(ost.str());
+ }
+
+ R2 direction = xmax - xmin;
+ const double length = l2Norm(direction);
+ direction *= 1. / length;
+
+ this->_checkBoundaryIsLine(direction, xmin, length, mesh);
+
+ return direction;
+}
+
+template <>
+template <>
+TinyVector<3>
+MeshLineNodeBoundary<3>::_getDirection(const Mesh<Connectivity<3>>& mesh)
+{
+ using R3 = TinyVector<3, double>;
+
+ std::array<R3, 6> bounds = this->_getBounds(mesh);
+
+ const R3& xmin = bounds[0];
+ const R3& ymin = bounds[1];
+ const R3& zmin = bounds[2];
+ const R3& xmax = bounds[3];
+ const R3& ymax = bounds[4];
+ const R3& zmax = bounds[5];
+
+ const R3 u = xmax - xmin;
+ const R3 v = ymax - ymin;
+ const R3 w = zmax - zmin;
+
+ R3 direction = u;
+ if (dot(v, v) > dot(direction, direction)) {
+ direction = v;
+ }
+
+ if (dot(w, w) > dot(direction, direction)) {
+ direction = w;
+ }
+
+ const double length = l2Norm(direction);
+ direction *= 1. / length;
+
+ this->_checkBoundaryIsLine(direction, xmin, length, mesh);
+
+ return direction;
+}
+
+template <size_t Dimension>
+MeshLineNodeBoundary<Dimension>
+getMeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+{
+ for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
+ ++i_ref_node_list) {
+ const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
+ const RefId& ref = ref_node_list.refId();
+ if (ref == boundary_descriptor) {
+ return MeshLineNodeBoundary<Dimension>{mesh, ref_node_list};
+ }
+ }
+ for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
+ ++i_ref_edge_list) {
+ const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
+ const RefId& ref = ref_edge_list.refId();
+ if (ref == boundary_descriptor) {
+ return MeshLineNodeBoundary<Dimension>{mesh, ref_edge_list};
+ }
+ }
+ for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
+ ++i_ref_face_list) {
+ const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == boundary_descriptor) {
+ return MeshLineNodeBoundary<Dimension>{mesh, ref_face_list};
+ }
+ }
+
+ std::ostringstream ost;
+ ost << "cannot find surface with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
+
+ throw NormalError(ost.str());
+}
+
+template MeshLineNodeBoundary<1> getMeshLineNodeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshLineNodeBoundary<2> getMeshLineNodeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshLineNodeBoundary<3> getMeshLineNodeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshLineNodeBoundary.hpp b/src/mesh/MeshLineNodeBoundary.hpp
new file mode 100644
index 000000000..0ae35d058
--- /dev/null
+++ b/src/mesh/MeshLineNodeBoundary.hpp
@@ -0,0 +1,65 @@
+#ifndef MESH_LINE_NODE_BOUNDARY_HPP
+#define MESH_LINE_NODE_BOUNDARY_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+
+template <size_t Dimension>
+class MeshLineNodeBoundary final : public MeshNodeBoundary<Dimension> // clazy:exclude=copyable-polymorphic
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const Rd m_direction;
+
+ template <size_t MeshDimension>
+ TinyVector<MeshDimension> _getDirection(const Mesh<Connectivity<MeshDimension>>&);
+
+ void _checkBoundaryIsLine(const TinyVector<Dimension, double>& direction,
+ const TinyVector<Dimension, double>& origin,
+ const double length,
+ const Mesh<Connectivity<Dimension>>& mesh) const;
+
+ public:
+ template <size_t MeshDimension>
+ friend MeshLineNodeBoundary<MeshDimension> getMeshLineNodeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
+ const IBoundaryDescriptor& boundary_descriptor);
+
+ PUGS_INLINE
+ const Rd&
+ direction() const
+ {
+ return m_direction;
+ }
+
+ MeshLineNodeBoundary& operator=(const MeshLineNodeBoundary&) = default;
+ MeshLineNodeBoundary& operator=(MeshLineNodeBoundary&&) = default;
+
+ private:
+ MeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefEdgeList& ref_edge_list)
+ : MeshNodeBoundary<Dimension>(mesh, ref_edge_list), m_direction(_getDirection(mesh))
+ {
+ ;
+ }
+
+ MeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list)
+ : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_direction(_getDirection(mesh))
+ {}
+
+ MeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefNodeList& ref_node_list)
+ : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_direction(_getDirection(mesh))
+ {}
+
+ public:
+ MeshLineNodeBoundary() = default;
+ MeshLineNodeBoundary(const MeshLineNodeBoundary&) = default;
+ MeshLineNodeBoundary(MeshLineNodeBoundary&&) = default;
+ ~MeshLineNodeBoundary() = default;
+};
+
+template <size_t Dimension>
+MeshLineNodeBoundary<Dimension> getMeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+ const IBoundaryDescriptor& boundary_descriptor);
+
+#endif // MESH_LINE_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshNodeBoundary.cpp b/src/mesh/MeshNodeBoundary.cpp
new file mode 100644
index 000000000..94fed7f5c
--- /dev/null
+++ b/src/mesh/MeshNodeBoundary.cpp
@@ -0,0 +1,323 @@
+#include <mesh/MeshNodeBoundary.hpp>
+
+template <>
+std::array<TinyVector<2>, 2>
+MeshNodeBoundary<2>::_getBounds(const Mesh<Connectivity<2>>& mesh) const
+{
+ using R2 = TinyVector<2, double>;
+
+ const NodeValue<const R2>& xr = mesh.xr();
+
+ std::array<R2, 2> bounds;
+ R2& xmin = bounds[0];
+ R2& xmax = bounds[1];
+
+ xmin = R2{std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+ xmax = R2{-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);
+ }
+ }
+
+ return bounds;
+}
+
+template <>
+std::array<TinyVector<3>, 6>
+MeshNodeBoundary<3>::_getBounds(const Mesh<Connectivity<3>>& mesh) const
+{
+ 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;
+ }
+ };
+
+ std::array<R3, 6> bounds;
+ R3& xmin = bounds[0];
+ R3& ymin = bounds[1];
+ R3& zmin = bounds[2];
+ R3& xmax = bounds[3];
+ R3& ymax = bounds[4];
+ 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;
+
+ xmax = -xmin;
+ ymax = xmax;
+ 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);
+ }
+ }
+
+ return bounds;
+}
+
+template <size_t Dimension>
+MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+ 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)) {
+ 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 <size_t Dimension>
+MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+ const RefEdgeList& ref_edge_list)
+ : m_boundary_name(ref_edge_list.refId().tagName())
+{
+ const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix();
+ const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+
+ auto edge_is_owned = mesh.connectivity().edgeIsOwned();
+
+ const Array<const EdgeId>& edge_list = ref_edge_list.list();
+
+ bool is_bad = false;
+ 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)) {
+ 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";
+ throw NormalError(ost.str());
+ }
+
+ Kokkos::vector<unsigned int> node_ids;
+ node_ids.reserve(2 * edge_list.size());
+ const auto& edge_to_node_matrix = mesh.connectivity().edgeToNodeMatrix();
+
+ for (size_t l = 0; l < edge_list.size(); ++l) {
+ const EdgeId edge_number = edge_list[l];
+ const auto& edge_nodes = edge_to_node_matrix[edge_number];
+
+ for (size_t r = 0; r < edge_nodes.size(); ++r) {
+ node_ids.push_back(edge_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 <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())
+{}
+
+template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefFaceList&);
+template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);
+template MeshNodeBoundary<3>::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&);
+
+template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefEdgeList&);
+template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefEdgeList&);
+template MeshNodeBoundary<3>::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefEdgeList&);
+
+template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefNodeList&);
+template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefNodeList&);
+template MeshNodeBoundary<3>::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefNodeList&);
+
+template <size_t Dimension>
+MeshNodeBoundary<Dimension>
+getMeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+{
+ for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
+ ++i_ref_node_list) {
+ const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
+ const RefId& ref = ref_node_list.refId();
+ if (ref == boundary_descriptor) {
+ return MeshNodeBoundary<Dimension>{mesh, ref_node_list};
+ }
+ }
+ for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
+ ++i_ref_edge_list) {
+ const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
+ const RefId& ref = ref_edge_list.refId();
+ if (ref == boundary_descriptor) {
+ return MeshNodeBoundary<Dimension>{mesh, ref_edge_list};
+ }
+ }
+ for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
+ ++i_ref_face_list) {
+ const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == boundary_descriptor) {
+ return MeshNodeBoundary<Dimension>{mesh, ref_face_list};
+ }
+ }
+
+ std::ostringstream ost;
+ ost << "cannot find surface with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
+
+ throw NormalError(ost.str());
+}
+
+template MeshNodeBoundary<1> getMeshNodeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshNodeBoundary<2> getMeshNodeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshNodeBoundary<3> getMeshNodeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index ab98ded0e..927e00d91 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -2,10 +2,8 @@
#define MESH_NODE_BOUNDARY_HPP
#include <Kokkos_Vector.hpp>
-#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
#include <mesh/Connectivity.hpp>
-#include <mesh/ConnectivityMatrix.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/Mesh.hpp>
@@ -21,13 +19,13 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
Array<const NodeId> m_node_list;
std::string m_boundary_name;
- template <typename MeshType>
- std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds(const MeshType& mesh) const;
+ std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds(
+ const Mesh<Connectivity<Dimension>>& mesh) const;
public:
- template <typename MeshType>
- friend MeshNodeBoundary<MeshType::Dimension> getMeshNodeBoundary(const MeshType& mesh,
- const IBoundaryDescriptor& boundary_descriptor);
+ template <size_t MeshDimension>
+ friend MeshNodeBoundary<MeshDimension> getMeshNodeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
+ const IBoundaryDescriptor& boundary_descriptor);
MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
@@ -39,119 +37,9 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
}
protected:
- 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();
-
- 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)) {
- 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& mesh, const RefEdgeList& ref_edge_list)
- : m_boundary_name(ref_edge_list.refId().tagName())
- {
- static_assert(Dimension == MeshType::Dimension);
- const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix();
- const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
-
- auto edge_is_owned = mesh.connectivity().edgeIsOwned();
-
- const Array<const EdgeId>& edge_list = ref_edge_list.list();
-
- bool is_bad = false;
- 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)) {
- 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";
- throw NormalError(ost.str());
- }
-
- Kokkos::vector<unsigned int> node_ids;
- node_ids.reserve(2 * edge_list.size());
- const auto& edge_to_node_matrix = mesh.connectivity().edgeToNodeMatrix();
-
- for (size_t l = 0; l < edge_list.size(); ++l) {
- const EdgeId edge_number = edge_list[l];
- const auto& edge_nodes = edge_to_node_matrix[edge_number];
-
- for (size_t r = 0; r < edge_nodes.size(); ++r) {
- node_ids.push_back(edge_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(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list);
+ MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefEdgeList& ref_edge_list);
+ MeshNodeBoundary(const Mesh<Connectivity<Dimension>>&, const RefNodeList& ref_node_list);
public:
MeshNodeBoundary(const MeshNodeBoundary&) = default;
@@ -161,761 +49,8 @@ class MeshNodeBoundary // clazy:exclude=copyable-polymorphic
virtual ~MeshNodeBoundary() = default;
};
-template <typename MeshType>
-MeshNodeBoundary<MeshType::Dimension>
-getMeshNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
-{
- for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
- ++i_ref_node_list) {
- const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
- const RefId& ref = ref_node_list.refId();
- if (ref == boundary_descriptor) {
- return MeshNodeBoundary<MeshType::Dimension>{mesh, ref_node_list};
- }
- }
- for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
- ++i_ref_edge_list) {
- const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
- const RefId& ref = ref_edge_list.refId();
- if (ref == boundary_descriptor) {
- return MeshNodeBoundary<MeshType::Dimension>{mesh, ref_edge_list};
- }
- }
- for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
- ++i_ref_face_list) {
- const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
- const RefId& ref = ref_face_list.refId();
- if (ref == boundary_descriptor) {
- return MeshNodeBoundary<MeshType::Dimension>{mesh, ref_face_list};
- }
- }
-
- std::ostringstream ost;
- ost << "cannot find surface with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
-
- throw NormalError(ost.str());
-}
-
-template <>
-template <typename MeshType>
-std::array<TinyVector<2>, 2>
-MeshNodeBoundary<2>::_getBounds(const MeshType& mesh) const
-{
- static_assert(MeshType::Dimension == 2);
-
- using R2 = TinyVector<2, double>;
-
- const NodeValue<const R2>& xr = mesh.xr();
-
- std::array<R2, 2> bounds;
- R2& xmin = bounds[0];
- R2& xmax = bounds[1];
-
- xmin = R2{std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
- xmax = R2{-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);
- }
- }
-
- return bounds;
-}
-
-template <>
-template <typename MeshType>
-std::array<TinyVector<3>, 6>
-MeshNodeBoundary<3>::_getBounds(const MeshType& mesh) const
-{
- 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;
- }
- };
-
- std::array<R3, 6> bounds;
- R3& xmin = bounds[0];
- R3& ymin = bounds[1];
- R3& zmin = bounds[2];
- R3& xmax = bounds[3];
- R3& ymax = bounds[4];
- 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;
-
- xmax = -xmin;
- ymax = xmax;
- 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);
- }
- }
-
- return bounds;
-}
-
-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();
-
- bool is_bad = false;
-
- parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
- const Rd& x = xr[this->m_node_list[r]];
- if (dot(x - origin, normal) > 1E-13 * length) {
- is_bad = true;
- }
- });
-
- if (parallel::allReduceOr(is_bad)) {
- 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>
- friend MeshFlatNodeBoundary<MeshType::Dimension> getMeshFlatNodeBoundary(
- const MeshType& mesh,
- const IBoundaryDescriptor& boundary_descriptor);
-
- private:
- 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))
- {
- ;
- }
-
- public:
- MeshFlatNodeBoundary() = default;
- MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
- MeshFlatNodeBoundary(MeshFlatNodeBoundary&&) = default;
- virtual ~MeshFlatNodeBoundary() = default;
-};
-
-template <typename MeshType>
-MeshFlatNodeBoundary<MeshType::Dimension>
-getMeshFlatNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
-{
- for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
- ++i_ref_node_list) {
- const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
- const RefId& ref = ref_node_list.refId();
- if (ref == boundary_descriptor) {
- return MeshFlatNodeBoundary<MeshType::Dimension>{mesh, ref_node_list};
- }
- }
- for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
- ++i_ref_face_list) {
- const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
- const RefId& ref = ref_face_list.refId();
- if (ref == boundary_descriptor) {
- return MeshFlatNodeBoundary<MeshType::Dimension>{mesh, ref_face_list};
- }
- }
-
- std::ostringstream ost;
- ost << "cannot find surface with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
-
- throw NormalError(ost.str());
-}
-
-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>;
-
- std::array<R2, 2> bounds = this->_getBounds(mesh);
-
- const R2& xmin = bounds[0];
- const R2& xmax = bounds[1];
-
- 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>;
-
- std::array<R3, 6> bounds = this->_getBounds(mesh);
-
- const R3& xmin = bounds[0];
- const R3& ymin = bounds[1];
- const R3& zmin = bounds[2];
- const R3& xmax = bounds[3];
- const R3& ymax = bounds[4];
- const R3& zmax = bounds[5];
-
- 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;
- }
-}
-
template <size_t Dimension>
-class MeshLineNodeBoundary : public MeshNodeBoundary<Dimension> // clazy:exclude=copyable-polymorphic
-{
- public:
- using Rd = TinyVector<Dimension, double>;
-
- private:
- const Rd m_direction;
-
- template <typename MeshType>
- PUGS_INLINE TinyVector<Dimension> _getDirection(const MeshType&);
-
- template <typename MeshType>
- PUGS_INLINE void
- _checkBoundaryIsLine(const TinyVector<Dimension, double>& direction,
- const TinyVector<Dimension, double>& origin,
- const double length,
- const MeshType& mesh) const
- {
- static_assert(MeshType::Dimension == Dimension);
- using Rdxd = TinyMatrix<Dimension>;
-
- const NodeValue<const Rd>& xr = mesh.xr();
-
- Rdxd P = Rdxd{identity} - tensorProduct(direction, direction);
-
- bool is_bad = false;
- parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
- const Rd& x = xr[this->m_node_list[r]];
- const Rd delta = x - origin;
- if (dot(P * delta, direction) > 1E-13 * length) {
- is_bad = true;
- }
- });
-
- if (parallel::allReduceOr(is_bad)) {
- std::ostringstream ost;
- ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
- << ": boundary is not a line!";
- throw NormalError(ost.str());
- }
- }
-
- public:
- template <typename MeshType>
- friend MeshLineNodeBoundary<MeshType::Dimension> getMeshLineNodeBoundary(
- const MeshType& mesh,
- const IBoundaryDescriptor& boundary_descriptor);
-
- PUGS_INLINE
- const Rd&
- direction() const
- {
- return m_direction;
- }
-
- MeshLineNodeBoundary& operator=(const MeshLineNodeBoundary&) = default;
- MeshLineNodeBoundary& operator=(MeshLineNodeBoundary&&) = default;
-
- private:
- template <typename MeshType>
- MeshLineNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
- : MeshNodeBoundary<Dimension>(mesh, ref_edge_list), m_direction(_getDirection(mesh))
- {
- ;
- }
-
- template <typename MeshType>
- MeshLineNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
- : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_direction(_getDirection(mesh))
- {
- ;
- }
-
- template <typename MeshType>
- MeshLineNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
- : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_direction(_getDirection(mesh))
- {
- ;
- }
-
- public:
- MeshLineNodeBoundary() = default;
- MeshLineNodeBoundary(const MeshLineNodeBoundary&) = default;
- MeshLineNodeBoundary(MeshLineNodeBoundary&&) = default;
- virtual ~MeshLineNodeBoundary() = default;
-};
-
-template <typename MeshType>
-MeshLineNodeBoundary<MeshType::Dimension>
-getMeshLineNodeBoundary(const MeshType& mesh, const IBoundaryDescriptor& boundary_descriptor)
-{
- for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
- ++i_ref_node_list) {
- const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
- const RefId& ref = ref_node_list.refId();
- if (ref == boundary_descriptor) {
- return MeshLineNodeBoundary<MeshType::Dimension>{mesh, ref_node_list};
- }
- }
- for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
- ++i_ref_edge_list) {
- const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
- const RefId& ref = ref_edge_list.refId();
- if (ref == boundary_descriptor) {
- return MeshLineNodeBoundary<MeshType::Dimension>{mesh, ref_edge_list};
- }
- }
- for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
- ++i_ref_face_list) {
- const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
- const RefId& ref = ref_face_list.refId();
- if (ref == boundary_descriptor) {
- return MeshLineNodeBoundary<MeshType::Dimension>{mesh, ref_face_list};
- }
- }
-
- std::ostringstream ost;
- ost << "cannot find surface with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
-
- throw NormalError(ost.str());
-}
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<1>
-MeshLineNodeBoundary<1>::_getDirection(const MeshType&)
-{
- throw UnexpectedError("MeshLineNodeBoundary::_getDirection make no sense in dimension 1");
-}
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<2>
-MeshLineNodeBoundary<2>::_getDirection(const MeshType& mesh)
-{
- using R2 = TinyVector<2, double>;
-
- std::array<R2, 2> bounds = this->_getBounds(mesh);
-
- const R2& xmin = bounds[0];
- const R2& xmax = bounds[1];
-
- if (xmin == xmax) {
- std::ostringstream ost;
- ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
- << ": unable to compute direction";
- throw NormalError(ost.str());
- }
-
- R2 direction = xmax - xmin;
- const double length = l2Norm(direction);
- direction *= 1. / length;
-
- this->_checkBoundaryIsLine(direction, xmin, length, mesh);
-
- return direction;
-}
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<3>
-MeshLineNodeBoundary<3>::_getDirection(const MeshType& mesh)
-{
- static_assert(MeshType::Dimension == 3);
-
- using R3 = TinyVector<3, double>;
-
- std::array<R3, 6> bounds = this->_getBounds(mesh);
-
- const R3& xmin = bounds[0];
- const R3& ymin = bounds[1];
- const R3& zmin = bounds[2];
- const R3& xmax = bounds[3];
- const R3& ymax = bounds[4];
- const R3& zmax = bounds[5];
-
- const R3 u = xmax - xmin;
- const R3 v = ymax - ymin;
- const R3 w = zmax - zmin;
-
- R3 direction = u;
- if (dot(v, v) > dot(direction, direction)) {
- direction = v;
- }
-
- if (dot(w, w) > dot(direction, direction)) {
- direction = w;
- }
-
- const double length = l2Norm(direction);
- direction *= 1. / length;
-
- this->_checkBoundaryIsLine(direction, xmin, length, mesh);
-
- return direction;
-}
+MeshNodeBoundary<Dimension> getMeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+ const IBoundaryDescriptor& boundary_descriptor);
#endif // MESH_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshRandomizer.hpp b/src/mesh/MeshRandomizer.hpp
index aa2258c59..fb29175a8 100644
--- a/src/mesh/MeshRandomizer.hpp
+++ b/src/mesh/MeshRandomizer.hpp
@@ -7,6 +7,8 @@
#include <mesh/Connectivity.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshLineNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 7bd97aef4..dedfbd2af 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -2,6 +2,7 @@
#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/ItemValueUtils.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
--
GitLab