From 65b43282dc1708b83db7f1ed332c852fd3e488f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com> Date: Tue, 20 Jul 2021 12:00:08 +0200 Subject: [PATCH] Add MeshFaceBoundary similarly to Mesh*NodeBoundary Plug it into AcousticSolver as a test --- src/mesh/CMakeLists.txt | 1 + src/mesh/MeshFaceBoundary.cpp | 36 +++++++++++++++ src/mesh/MeshFaceBoundary.hpp | 55 +++++++++++++++++++++++ src/scheme/AcousticSolver.cpp | 83 ++++++++++++++++++----------------- 4 files changed, 134 insertions(+), 41 deletions(-) create mode 100644 src/mesh/MeshFaceBoundary.cpp create mode 100644 src/mesh/MeshFaceBoundary.hpp diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt index 2d429282e..97b03eda3 100644 --- a/src/mesh/CMakeLists.txt +++ b/src/mesh/CMakeLists.txt @@ -17,6 +17,7 @@ add_library( LogicalConnectivityBuilder.cpp MeshBuilderBase.cpp MeshDataManager.cpp + MeshFaceBoundary.cpp MeshFlatNodeBoundary.cpp MeshLineNodeBoundary.cpp MeshNodeBoundary.cpp diff --git a/src/mesh/MeshFaceBoundary.cpp b/src/mesh/MeshFaceBoundary.cpp new file mode 100644 index 000000000..3690c5c0c --- /dev/null +++ b/src/mesh/MeshFaceBoundary.cpp @@ -0,0 +1,36 @@ +#include <mesh/MeshFaceBoundary.hpp> + +#include <mesh/Connectivity.hpp> +#include <mesh/Mesh.hpp> +#include <utils/Messenger.hpp> + +template <size_t Dimension> +MeshFaceBoundary<Dimension>::MeshFaceBoundary(const Mesh<Connectivity<Dimension>>&, const RefFaceList& ref_face_list) + : m_face_list(ref_face_list.list()), m_boundary_name(ref_face_list.refId().tagName()) +{} + +template MeshFaceBoundary<2>::MeshFaceBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&); +template MeshFaceBoundary<3>::MeshFaceBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&); + +template <size_t Dimension> +MeshFaceBoundary<Dimension> +getMeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor) +{ + 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 MeshFaceBoundary<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 MeshFaceBoundary<1> getMeshFaceBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&); +template MeshFaceBoundary<2> getMeshFaceBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&); +template MeshFaceBoundary<3> getMeshFaceBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&); diff --git a/src/mesh/MeshFaceBoundary.hpp b/src/mesh/MeshFaceBoundary.hpp new file mode 100644 index 000000000..afd409d93 --- /dev/null +++ b/src/mesh/MeshFaceBoundary.hpp @@ -0,0 +1,55 @@ +#ifndef MESH_FACE_BOUNDARY_HPP +#define MESH_FACE_BOUNDARY_HPP + +#include <algebra/TinyVector.hpp> +#include <mesh/IBoundaryDescriptor.hpp> +#include <mesh/ItemValue.hpp> +#include <mesh/RefItemList.hpp> +#include <utils/Array.hpp> + +template <size_t Dimension> +class Connectivity; + +template <typename ConnectivityType> +class Mesh; + +template <size_t Dimension> +class MeshFaceBoundary // clazy:exclude=copyable-polymorphic +{ + protected: + Array<const FaceId> m_face_list; + std::string m_boundary_name; + + std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds( + const Mesh<Connectivity<Dimension>>& mesh) const; + + public: + template <size_t MeshDimension> + friend MeshFaceBoundary<MeshDimension> getMeshFaceBoundary(const Mesh<Connectivity<MeshDimension>>& mesh, + const IBoundaryDescriptor& boundary_descriptor); + + MeshFaceBoundary& operator=(const MeshFaceBoundary&) = default; + MeshFaceBoundary& operator=(MeshFaceBoundary&&) = default; + + const Array<const FaceId>& + faceList() const + { + return m_face_list; + } + + protected: + MeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list); + + public: + MeshFaceBoundary(const MeshFaceBoundary&) = default; + MeshFaceBoundary(MeshFaceBoundary&&) = default; + + MeshFaceBoundary() = default; + virtual ~MeshFaceBoundary() = default; +}; + +template <size_t Dimension> +MeshFaceBoundary<Dimension> getMeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, + const IBoundaryDescriptor& boundary_descriptor); + +#endif // MESH_FACE_BOUNDARY_HPP diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp index dedfbd2af..5b4d9cb57 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/MeshFaceBoundary.hpp> #include <mesh/MeshFlatNodeBoundary.hpp> #include <mesh/MeshNodeBoundary.hpp> #include <scheme/DirichletBoundaryConditionDescriptor.hpp> @@ -91,8 +92,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler { MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh); - const NodeValuePerCell<const TinyVector<Dimension>> Cjr = mesh_data.Cjr(); - const NodeValuePerCell<const TinyVector<Dimension>> njr = mesh_data.njr(); + const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr(); + const NodeValuePerCell<const Rd> njr = mesh_data.njr(); NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()}; @@ -253,12 +254,12 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler switch (bc_descriptor->type()) { case IBoundaryConditionDescriptor::Type::symmetry: { - bc_list.push_back( + bc_list.emplace_back( SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()))); break; } case IBoundaryConditionDescriptor::Type::fixed: { - bc_list.push_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()))); break; } case IBoundaryConditionDescriptor::Type::dirichlet: { @@ -268,36 +269,35 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler MeshNodeBoundary<Dimension> mesh_node_boundary = getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor()); - Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>( - TinyVector<Dimension>)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(), - mesh.xr(), mesh_node_boundary.nodeList()); + Array<const Rd> value_list = + InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(), + mesh.xr(), + mesh_node_boundary.nodeList()); - bc_list.push_back(VelocityBoundaryCondition{mesh_node_boundary, value_list}); + bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, value_list}); } else if (dirichlet_bc_descriptor.name() == "pressure") { - for (size_t i_ref_face_list = 0; - i_ref_face_list < mesh.connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { - const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list); - const RefId& ref = ref_face_list.refId(); - if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) { - const auto& face_list = ref_face_list.list(); - - const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId(); - - Array<const double> face_values = [&] { - if constexpr (Dimension == 1) { - return InterpolateItemValue<double( - TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh.xr(), face_list); - } else { - MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh); - - return InterpolateItemValue<double( - TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), face_list); - } - }(); - - bc_list.push_back(PressureBoundaryCondition{face_list, face_values}); - } + const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId(); + + if constexpr (Dimension == 1) { + MeshNodeBoundary<Dimension> mesh_node_boundary = + getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()); + + Array<const double> node_values = + InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(), + mesh_node_boundary.nodeList()); + + bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values}); + } else { + MeshFaceBoundary<Dimension> mesh_face_boundary = + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()); + + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh); + Array<const double> face_values = + InterpolateItemValue<double(Rd)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), + mesh_face_boundary.faceList()); + bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values}); } + } else { is_valid_boundary_condition = false; } @@ -499,7 +499,7 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyPressureBC(const MeshTyp const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); - const auto& node_list = bc.faceList(); + const auto& node_list = bc.nodeList(); const auto& value_list = bc.valueList(); parallel_for( node_list.size(), PUGS_LAMBDA(size_t i_node) { @@ -634,14 +634,14 @@ template <size_t Dimension> class AcousticSolverHandler::AcousticSolver<Dimension>::PressureBoundaryCondition { private: + const MeshFaceBoundary<Dimension> m_mesh_face_boundary; const Array<const double> m_value_list; - const Array<const FaceId> m_face_list; public: const Array<const FaceId>& faceList() const { - return m_face_list; + return m_mesh_face_boundary.faceList(); } const Array<const double>& @@ -650,8 +650,9 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::PressureBoundaryConditio return m_value_list; } - PressureBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list) - : m_value_list{value_list}, m_face_list{face_list} + PressureBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary, + const Array<const double>& value_list) + : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list} {} ~PressureBoundaryCondition() = default; @@ -661,14 +662,14 @@ template <> class AcousticSolverHandler::AcousticSolver<1>::PressureBoundaryCondition { private: + const MeshNodeBoundary<1> m_mesh_node_boundary; const Array<const double> m_value_list; - const Array<const NodeId> m_face_list; public: const Array<const NodeId>& - faceList() const + nodeList() const { - return m_face_list; + return m_mesh_node_boundary.nodeList(); } const Array<const double>& @@ -677,8 +678,8 @@ class AcousticSolverHandler::AcousticSolver<1>::PressureBoundaryCondition return m_value_list; } - PressureBoundaryCondition(const Array<const NodeId>& face_list, const Array<const double>& value_list) - : m_value_list{value_list}, m_face_list{face_list} + PressureBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const double>& value_list) + : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list} {} ~PressureBoundaryCondition() = default; -- GitLab