Construction of second member for Dirichlet BC.

src/scheme/ScalarNodalScheme.cpp

@@ -176,7 +176,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                                                                                                       mesh_data.xl(),
                                                                                                       mesh_face_boundary
                                                                                                         .faceList());
-
           boundary_condition_list.push_back(
             DirichletBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
 
@@ -254,7 +253,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             using T = std::decay_t<decltype(bc)>;
             if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
               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];
@@ -495,7 +493,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
               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];
 
-              size_t i_cell             = node_to_cell[0];
+              size_t i_cell             = 0;
               const CellId cell_id      = node_to_cell[i_cell];
               const unsigned int i_node = node_local_number_in_its_cell[i_cell];
 
@@ -520,7 +518,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
                 Cjr2[0] = -0.5 * (xrp2[1] - xrp1[1]);
                 Cjr2[1] = 0.5 * (xrp2[0] - xrp1[0]);
               } else {
-                std::cout << "Dimension " << Dimension << " not implemented. \n";
                 throw NotImplementedError("The scheme is not implemented in this dimension.");
               }
 
@@ -721,17 +718,33 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
             b[cell_id] += 0.5 * node_boundary_values[node_id] * sum_mes_l;
           }
         } else if (node_is_dirichlet[node_id]) {
-          if (is_boundary_node[node_id] && (not node_is_corner[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];
+          if (not node_is_corner[node_id]) {
+            for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+              const CellId cell_id_j = node_to_cell[i_cell_j];
+              const size_t i_node_j  = node_local_number_in_its_cells[i_cell_j];
+
+              for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+                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] += 0;
+                b[cell_id_j] += dot(node_boundary_values[node_id] * kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+                                    Cjr(cell_id_j, i_node_j));
               }
+            }
+
           } else if (node_is_corner[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];
+            for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
+              const CellId cell_id_j = node_to_cell[i_cell_j];
+              const size_t i_node_j  = node_local_number_in_its_cells[i_cell_j];
+
+              for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
+                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] += 0;
+                b[cell_id_j] +=
+                  dot(node_boundary_values[node_id] * corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+                      Cjr(cell_id_j, i_node_j));
+              }
             }
           }
         };