Fix renormalization v

src/scheme/ScalarNodalScheme.cpp

@@ -143,7 +143,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
   }
 
   ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
-                    const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_bound,
+                    const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_b,
                     const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
                     const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
                     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
@@ -409,7 +409,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
 
     {
       CellValue<const TinyMatrix<Dimension>> cell_kappaj  = cell_k->cellValues();
-      CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_bound->cellValues();
+      CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_b->cellValues();
 
       const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
         NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
@@ -568,8 +568,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       const NodeValue<const TinyVector<Dimension>> v = [&] {
         NodeValue<TinyVector<Dimension>> compute_v{mesh->connectivity()};
         parallel_for(
-          mesh->numberOfNodes(),
-          PUGS_LAMBDA(NodeId node_id) { compute_v[node_id] = node_kapparb[node_id] * exterior_normal[node_id]; });
+          mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+            if (is_boundary_node[node_id]) {
+              compute_v[node_id] = 1. / l2Norm(node_kapparb[node_id] * exterior_normal[node_id]) *
+                                   node_kapparb[node_id] * exterior_normal[node_id];
+            }
+          });
         return compute_v;
       }();
 
@@ -704,7 +708,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
 
               b[cell_id] += node_boundary_values[node_id] *
                             (1. / (sum_theta[node_id] * l2Norm(node_kapparb[node_id] * n)) *
-                               dot(kapparb_invBetar[node_id] * sum_Cjr[node_id], P * Cjr(cell_id, i_node)) -
+                               dot(kapparb_invBetar[node_id] * sum_Cjr[node_id], P * Cjr(cell_id, i_node)) +
                              dot(Cjr(cell_id, i_node), n));
             }
           } else if (node_is_corner[node_id]) {
@@ -727,7 +731,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                 const CellId cell_id_k = node_to_cell[i_cell_k];
                 const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
 
-                b[cell_id_j] += dot(node_boundary_values[node_id] * kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+                b[cell_id_j] +=
+                  dot(node_boundary_values[node_id] * kapparb_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
                       Cjr(cell_id_j, i_node_j));
               }
             }
@@ -742,7 +747,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                 const size_t i_node_k  = node_local_number_in_its_cells[i_cell_k];
 
                 b[cell_id_j] +=
-                  dot(node_boundary_values[node_id] * corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+                  dot(node_boundary_values[node_id] * corner_kapparb_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
                       Cjr(cell_id_j, i_node_j));
               }
             }
@@ -771,7 +776,7 @@ ScalarNodalSchemeHandler::solution() const
 }
 
 ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
-  const std::shared_ptr<const IDiscreteFunction>& cell_k_bound,
+  const std::shared_ptr<const IDiscreteFunction>& cell_k_b,
   const std::shared_ptr<const IDiscreteFunction>& cell_k,
   const std::shared_ptr<const IDiscreteFunction>& f,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
@@ -781,7 +786,7 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
     throw NormalError("primal discrete functions are not defined on the same mesh");
   }
 
-  checkDiscretizationType({cell_k_bound, cell_k, f}, DiscreteFunctionType::P0);
+  checkDiscretizationType({cell_k_b, cell_k, f}, DiscreteFunctionType::P0);
 
   switch (i_mesh->dimension()) {
   case 1: {
@@ -797,7 +802,7 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
 
     m_scheme =
       std::make_unique<ScalarNodalScheme<2>>(mesh,
-                                             std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
+                                             std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_b),
                                              std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
                                              std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
                                              bc_descriptor_list);
@@ -812,7 +817,7 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
 
     m_scheme =
       std::make_unique<ScalarNodalScheme<3>>(mesh,
-                                             std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
+                                             std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_b),
                                              std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
                                              std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
                                              bc_descriptor_list);