From 75ec66e35419ce2850021f8e232999e5ad526b4c Mon Sep 17 00:00:00 2001 From: Drouard <axelle.drouard@orange.fr> Date: Thu, 21 Jul 2022 10:37:59 +0200 Subject: [PATCH] Change to Dirichlet BC at the nodes and Neumann BC at the faces --- src/scheme/ScalarNodalScheme.cpp | 112 ++++++++++++++++++++----------- 1 file changed, 74 insertions(+), 38 deletions(-) diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index 4639f58dd..25bce53f9 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -67,18 +67,18 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand const Array<const NodeId> m_node_list; public: - const Array<const FaceId>& - faceList() const - { - return m_face_list; - } - const Array<const NodeId>& nodeList() const { return m_node_list; } + const Array<const FaceId>& + faceList() const + { + return m_face_list; + } + const Array<const double>& valueList() const { @@ -90,7 +90,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand const Array<const double>& value_list) : m_value_list{value_list}, m_face_list{face_list}, m_node_list{node_list} { - Assert(m_value_list.size() == m_node_list.size()); + Assert(m_value_list.size() == m_face_list.size()); } ~NeumannBoundaryCondition() = default; @@ -101,13 +101,13 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand private: const Array<const double> m_coef_list; const Array<const double> m_value_list; - const Array<const NodeId> m_node_list; + const Array<const FaceId> m_face_list; public: - const Array<const NodeId>& - nodeList() const + const Array<const FaceId>& + faceList() const { - return m_node_list; + return m_face_list; } const Array<const double>& @@ -123,13 +123,13 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand } public: - FourierBoundaryCondition(const Array<const NodeId>& node_list, + FourierBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& coef_list, const Array<const double>& value_list) - : m_coef_list{coef_list}, m_value_list{value_list}, m_node_list{node_list} + : m_coef_list{coef_list}, m_value_list{value_list}, m_face_list{face_list} { - Assert(m_coef_list.size() == m_node_list.size()); - Assert(m_value_list.size() == m_node_list.size()); + Assert(m_coef_list.size() == m_face_list.size()); + Assert(m_value_list.size() == m_face_list.size()); } ~FourierBoundaryCondition() = default; @@ -202,12 +202,13 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand getMeshNodeBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor()); const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId(); - // MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); Array<const double> value_list = - InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(), - mesh_node_boundary - .nodeList()); + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, + mesh_data.xl(), + mesh_face_boundary + .faceList()); boundary_condition_list.push_back( NeumannBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list}); @@ -273,12 +274,13 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand getMeshNodeBoundary(*mesh, neumann_bc_prev_descriptor.boundaryDescriptor()); const FunctionSymbolId g_id = neumann_bc_prev_descriptor.rhsSymbolId(); - // MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); Array<const double> value_list = - InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(), - mesh_node_boundary - .nodeList()); + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, + mesh_data.xl(), + mesh_face_boundary + .faceList()); boundary_prev_condition_list.push_back( NeumannBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list}); @@ -311,8 +313,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand const auto is_boundary_node = mesh->connectivity().isBoundaryNode(); const auto is_boundary_face = mesh->connectivity().isBoundaryFace(); - const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix(); - // const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); + const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix(); + const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); const auto& node_local_numbers_in_their_cells = mesh->connectivity().nodeLocalNumbersInTheirCells(); const CellValue<const double> Vj = mesh_data.Vj(); @@ -327,10 +329,14 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand [&](auto&& bc) { using T = std::decay_t<decltype(bc)>; if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) { - const auto& node_list = bc.nodeList(); - for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - const NodeId node_id = node_list[i_node]; - compute_node_is_neumann[node_id] = true; + const auto& face_list = bc.faceList(); + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; + const auto& face_nodes = face_to_node_matrix[face_id]; + for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) { + const NodeId node_id = face_nodes[i_node]; + compute_node_is_neumann[node_id] = true; + } } } }, @@ -471,12 +477,25 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand [&](auto&& bc) { using T = std::decay_t<decltype(bc)>; if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) { - const auto& node_list = bc.nodeList(); + const auto& face_list = bc.faceList(); const auto& value_list = bc.valueList(); - for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - const NodeId node_id = node_list[i_node]; - compute_node_boundary_values[node_id] = value_list[i_node]; + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; + const auto& face_nodes = face_to_node_matrix[face_id]; + for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) { + const NodeId node_id = face_nodes[i_node]; + if (not node_is_dirichlet[node_id]) { + if (node_is_angle[node_id]) { + compute_node_boundary_values[node_id] += + value_list[i_face] * std::abs(dot(nl[face_id], exterior_normal[node_id])); + + } else { + compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id]; + sum_mes_l[node_id] += mes_l[face_id]; + } + } + } } } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) { @@ -492,7 +511,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand }, boundary_condition); } - // std::cout << node_id << " " << compute_node_boundary_values[node_id] << "\n"; + for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { + if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id]) and not node_is_angle[node_id]) { + compute_node_boundary_values[node_id] /= sum_mes_l[node_id]; + } + // std::cout << node_id << " " << compute_node_boundary_values[node_id] << "\n"; + } return compute_node_boundary_values; }(); @@ -506,13 +530,20 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand [&](auto&& bc) { using T = std::decay_t<decltype(bc)>; if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) { - const auto& node_list = bc.nodeList(); + const auto& face_list = bc.faceList(); const auto& value_list = bc.valueList(); - for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - const NodeId node_id = node_list[i_node]; + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; + const auto& face_nodes = face_to_node_matrix[face_id]; - compute_node_boundary_values[node_id] = value_list[i_node]; + for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) { + const NodeId node_id = face_nodes[i_node]; + if (not node_is_dirichlet[node_id]) { + compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id]; + sum_mes_l[node_id] += mes_l[face_id]; + } + } } } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) { const auto& node_list = bc.nodeList(); @@ -527,6 +558,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand }, boundary_condition); } + for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { + if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id])) { + compute_node_boundary_values[node_id] /= sum_mes_l[node_id]; + } + } return compute_node_boundary_values; }(); -- GitLab