Fix Neumann BC matrix on the border.

src/scheme/ScalarNodalScheme.cpp

@@ -586,7 +586,7 @@ 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];
 
-              S(cell_id_j, cell_id_k) +=
+              S(cell_id_j, cell_id_k) -=
                 dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
             }
           }
@@ -604,10 +604,10 @@ 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];
 
-              S(cell_id_j, cell_id_k) +=
+              S(cell_id_j, cell_id_k) -=
                 dot(kappar_invBetar[node_id] *
                       (Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
-                    Q * Cjr(cell_id_k, i_node_j));
+                    Q * Cjr(cell_id_j, i_node_j));
             }
           }
         } else if ((node_is_dirichlet[node_id]) && (not node_is_corner[node_id])) {
@@ -619,7 +619,7 @@ 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];
 
-              S(cell_id_j, cell_id_k) +=
+              S(cell_id_j, cell_id_k) -=
                 dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
             }
           }
@@ -632,7 +632,7 @@ 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];
 
-              S(cell_id_j, cell_id_k) +=
+              S(cell_id_j, cell_id_k) -=
                 dot(corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
             }
           }
@@ -641,7 +641,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
 
       for (CellId cell_id_j = 0; cell_id_j < mesh->numberOfCells(); ++cell_id_j) {
         for (CellId cell_id_k = 0; cell_id_k < mesh->numberOfCells(); ++cell_id_k) {
-          S(cell_id_j, cell_id_k) *= dt / Vj[cell_id_j];
+          S(cell_id_j, cell_id_k) = -dt / Vj[cell_id_j] * S(cell_id_j, cell_id_k);
         }
       };
 
@@ -720,8 +720,38 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               }
             }
           }
+        }
       };
 
+      // TESTS //
+      NodeValue<TinyVector<Dimension>> gradPr{mesh->connectivity()};
+      NodeValue<double> Pr{mesh->connectivity()};
+      for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+        gradPr[node_id]                            = zero;
+        const auto& node_to_cell                   = node_to_cell_matrix[node_id];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+        if (not is_boundary_node[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 size_t i_node  = node_local_number_in_its_cells[i_cell];
+            gradPr[node_id] += Ph[cell_id] * Cjr(cell_id, i_node);
+          }
+          gradPr[node_id] = -inverse(node_betar[node_id]) * gradPr[node_id];
+          std::cout << "node_id : " << node_id << "; gradPr = " << gradPr[node_id] << "\n";
+        } else if (is_boundary_node[node_id] and not node_is_corner[node_id]) {
+          Pr[node_id] = 0;
+          for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+            const CellId cell_id = node_to_cell[i_cell];
+            const size_t i_node  = node_local_number_in_its_cells[i_cell];
+            Pr[node_id] += theta(cell_id, i_node) * Ph[cell_id];
+          }
+          Pr[node_id] += node_boundary_values[node_id] / l2Norm(node_kappar[node_id] * exterior_normal[node_id]);
+          Pr[node_id] *= 1. / sum_theta[node_id];
+          std::cout << "node_id :" << node_id << "; Pr = " << Pr[node_id] << "\n";
+        }
+      };
+      // END TESTS//
+
       Vector<double> T{mesh->numberOfCells()};
       T = zero;
 
@@ -729,6 +759,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
       parallel_for(
         mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { B[cell_id] = dt * b[cell_id] + Ph[cell_id]; });
 
+      Vector<double> Pvec{mesh->numberOfCells()};
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Pvec[cell_id] = Ph[cell_id]; });
+      std::cout << "B = " << B << "\n"
+                << "AP = " << A * Pvec << "\n";
+
       LinearSolver solver;
       solver.solveLocalSystem(A, T, B);
 
@@ -738,7 +774,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; });
     }
   }
-  }
 };
 
 std::shared_ptr<const IDiscreteFunction>