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,6 +720,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               }
             }
           }
+        }
       };
 
       Vector<double> T{mesh->numberOfCells()};
@@ -738,7 +739,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
         mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; });
     }
   }
-  }
 };
 
 std::shared_ptr<const IDiscreteFunction>