Change to Dirichlet BC at the nodes and Neumann BC at the faces

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});
@@ -312,7 +314,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     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& 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,12 +329,16 @@ 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];
+              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;
                 }
               }
+            }
           },
           boundary_condition);
       }
@@ -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);
       }
+      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;
     }();