Fix corner detection

src/scheme/ScalarNodalScheme.cpp

@@ -306,38 +306,38 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
           if (is_boundary_node[node_id]) {
             const auto& node_to_cell                  = node_to_cell_matrix[node_id];
-
+            const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
             for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+              const unsigned int i_node = node_local_number_in_its_cell[i_cell];
               const CellId cell_id      = node_to_cell[i_cell];
               const auto& cell_nodes    = cell_to_node_matrix[cell_id];
-              const NodeId prev_node_id = cell_to_node_matrix[cell_id][(node_id - 1) % cell_nodes.size()];
-              const NodeId next_node_id = cell_to_node_matrix[cell_id][(node_id + 1) % cell_nodes.size()];
-
+              const NodeId prev_node_id = cell_to_node_matrix[cell_id][(i_node - 1) % cell_nodes.size()];
+              const NodeId next_node_id = cell_to_node_matrix[cell_id][(i_node + 1) % cell_nodes.size()];
               if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) {
                 compute_node_is_corner[node_id] = true;
               }
             }
           }
         });
-
       return compute_node_is_corner;
     }();
 
     const NodeValue<const TinyVector<Dimension>> exterior_normal = [&] {
-      NodeValue<TinyVector<Dimension>> compute_exterior_normal;
-      parallel_for(
-        mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
-          compute_exterior_normal[node_id]          = zero;
+      NodeValue<TinyVector<Dimension>> compute_exterior_normal{mesh->connectivity()};
+      compute_exterior_normal.fill(zero);
+      for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+        if (is_boundary_node[node_id]) {
           const auto& node_to_cell                  = node_to_cell_matrix[node_id];
           const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
           for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
             const CellId cell_id      = node_to_cell[i_cell];
-            const unsigned int i_node_loc = node_local_number_in_its_cell[cell_id];
-            compute_exterior_normal[node_id] += Cjr(cell_id, i_node_loc);
+            const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+            compute_exterior_normal[node_id] += Cjr(cell_id, i_node);
+          }
+          const double norm_exterior_normal = l2Norm(compute_exterior_normal[node_id]);
+          compute_exterior_normal[node_id] *= 1. / norm_exterior_normal;
+        }
       }
-          compute_exterior_normal[node_id] =
-            1. / l2Norm(compute_exterior_normal[node_id]) * compute_exterior_normal[node_id];
-        });
       return compute_exterior_normal;
     }();
 
@@ -352,8 +352,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
     }();
 
     const NodeValue<const double> node_boundary_values = [&] {
-      NodeValue<double> compute_node_boundary_values;
-      NodeValue<double> sum_mes_l;
+      NodeValue<double> compute_node_boundary_values{mesh->connectivity()};
+      NodeValue<double> sum_mes_l{mesh->connectivity()};
       compute_node_boundary_values.fill(0);
       sum_mes_l.fill(0);
       for (const auto& boundary_condition : boundary_condition_list) {
@@ -376,11 +376,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                   }
                 }
               }
-              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];
-                }
-              }
+
             } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
               const auto& face_list  = bc.faceList();
               const auto& value_list = bc.valueList();
@@ -399,15 +395,17 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                   }
                 }
               }
-              for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
-                if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
-                  compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
-                }
-              }
             }
           },
           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];
+        } else if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
+          compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
+        }
+      }
       return compute_node_boundary_values;
     }();
 
@@ -416,7 +414,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_bound->cellValues();
 
       const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
-        NodeValue<TinyMatrix<Dimension>> kappa;
+        NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
             kappa[node_id]           = zero;
@@ -434,7 +432,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       }();
 
       const NodeValue<const TinyMatrix<Dimension>> node_kapparb = [&] {
-        NodeValue<TinyMatrix<Dimension>> kappa;
+        NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
             kappa[node_id]           = zero;
@@ -452,41 +450,45 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       }();
 
       const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
-        NodeValue<TinyMatrix<Dimension>> beta;
+        NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()};
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
             const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
             const auto& node_to_cell                  = node_to_cell_matrix[node_id];
             beta[node_id]                             = zero;
-            for (size_t j = 0; j < node_to_cell.size(); ++j) {
-              const CellId J       = node_to_cell[j];
-              const unsigned int R = node_local_number_in_its_cell[j];
-              beta[node_id] += tensorProduct(Cjr(J, R), xr[node_id] - xj[J]);
+            for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+              const CellId cell_id      = node_to_cell[i_cell];
+              const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+              beta[node_id] += tensorProduct(Cjr(cell_id, i_node), xr[node_id] - xj[cell_id]);
             }
           });
         return beta;
       }();
 
       const NodeValue<const TinyMatrix<Dimension>> kappar_invBetar = [&] {
-        NodeValue<TinyMatrix<Dimension>> kappa_invBeta;
+        NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+            if (not node_is_corner[node_id]) {
               kappa_invBeta[node_id] = node_kappar[node_id] * inverse(node_betar[node_id]);
+            }
           });
         return kappa_invBeta;
       }();
 
       const NodeValue<const TinyMatrix<Dimension>> kapparb_invBetar = [&] {
-        NodeValue<TinyMatrix<Dimension>> kappa_invBeta;
+        NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()};
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+            if (not node_is_corner[node_id]) {
               kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(node_betar[node_id]);
+            }
           });
         return kappa_invBeta;
       }();
 
       const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] {
-        NodeValue<TinyVector<Dimension>> compute_sum_Cjr;
+        NodeValue<TinyVector<Dimension>> compute_sum_Cjr{mesh->connectivity()};
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
             const auto& node_to_cell                  = node_to_cell_matrix[node_id];
@@ -502,7 +504,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       }();
 
       const NodeValue<const TinyVector<Dimension>> v = [&] {
-        NodeValue<TinyVector<Dimension>> compute_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]; });
@@ -510,20 +512,23 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       }();
 
       const NodeValuePerCell<const double> theta = [&] {
-        NodeValuePerCell<double> compute_theta;
-        parallel_for(
-          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        NodeValuePerCell<double> compute_theta{mesh->connectivity()};
+        // parallel_for(
+        //  mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
           const auto& cell_nodes = cell_to_node_matrix[cell_id];
           for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
             const NodeId node_id = cell_nodes[i_node];
+            std::cout << "node_id = " << node_id << "\n";
             compute_theta(cell_id, i_node) = dot(inverse(node_betar[node_id]) * Cjr(cell_id, i_node), v[node_id]);
           }
-          });
+        }
+        //          });
         return compute_theta;
       }();
 
       const NodeValue<const double> sum_theta = [&] {
-        NodeValue<double> compute_sum_theta;
+        NodeValue<double> compute_sum_theta{mesh->connectivity()};
         parallel_for(
           mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
             const auto& node_to_cell                  = node_to_cell_matrix[node_id];