Fix hybrid matrix and second member for Neumann BC

src/scheme/ScalarHybridScheme.cpp

@@ -1077,7 +1077,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
               const FaceId face_id     = cell_to_face[i_face];
               const auto& face_to_node = face_to_node_matrix[face_id];
-
+              if (not is_boundary_face[face_id]) {
                 if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
                   S1(cell_id_j, cell_id_k) +=
                     lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
@@ -1089,6 +1089,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
               }
             }
           }
+        }
       } else if ((node_is_neumann[node_id]) && (not node_is_corner[node_id]) && (not node_is_dirichlet[node_id])) {
         const TinyMatrix<Dimension> I   = identity;
         const TinyVector<Dimension> n   = exterior_normal[node_id];
@@ -1113,6 +1114,29 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
               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_j, i_node_j));
+
+            const auto& cell_to_face = cell_to_face_matrix[cell_id_j];
+
+            for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+              const FaceId face_id     = cell_to_face[i_face];
+              const auto& face_to_node = face_to_node_matrix[face_id];
+              if (not is_boundary_face[face_id]) {
+                if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
+                  S1(cell_id_j, cell_id_k) +=
+                    lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+                    dot(inverse(node_betar[node_id]) *
+                          (Cjr(cell_id_k, i_node_k) -
+                           theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
+                        xj1_xj2_orth(cell_id_j, i_face));
+                  S2(cell_id_j, cell_id_k) -=
+                    lambda * dt * (cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+                    dot(inverse(node_betar[node_id]) *
+                          (Cjr(cell_id_k, i_node_k) -
+                           theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
+                        xj1_xj2_orth(cell_id_j, i_face));
+                }
+              }
+            }
           }
         }
       } else if ((node_is_dirichlet[node_id]) && (not node_is_corner[node_id])) {
@@ -1136,7 +1160,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
               const FaceId face_id     = cell_to_face[i_face];
               const auto& face_to_node = face_to_node_matrix[face_id];
-              if (face_is_dirichlet[face_id]) {
+              if (face_is_dirichlet[face_id] or not is_boundary_face[face_id]) {
                 if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
                   S1(cell_id_j, cell_id_k) +=
                     lambda * dt * (1. - cn_coeff / 2.) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
@@ -1278,7 +1302,28 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
                                (1. / (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
                                   dot(kappar_invBetar[node_id] * sum_Cjr[node_id], Q * Cjr(cell_id, i_node)) +
                                 dot(Cjr(cell_id, i_node), n));
+
+            const auto& cell_to_face = cell_to_face_matrix[cell_id];
+
+            for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+              const FaceId face_id     = cell_to_face[i_face];
+              const auto& face_to_node = face_to_node_matrix[face_id];
+              if (not is_boundary_face[face_id]) {
+                if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
+                  b[cell_id] += lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id, i_face)[0] *
+                                node_boundary_values[node_id] /
+                                (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
+                                dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id, i_face));
+                  // std::cout << "xl - xj2 = " << xl[face_id] -
+                  b_prev[cell_id] +=
+                    lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id, i_face)[0] *
+                    node_boundary_prev_values[node_id] / (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
+                    dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id, i_face));
+                }
               }
+            }
+          }
+
         } else if (node_is_corner[node_id]) {
           const auto& node_to_face = node_to_face_matrix[node_id];
           const CellId cell_id     = node_to_cell[0];
@@ -1313,7 +1358,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
             for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
               const FaceId face_id     = cell_to_face[i_face];
               const auto& face_to_node = face_to_node_matrix[face_id];
-              if (face_is_dirichlet[face_id]) {
+              if (face_is_dirichlet[face_id] or not is_boundary_face[face_id]) {
                 if (face_to_node[0] == node_id or face_to_node[1] == node_id) {
                   b[cell_id_j] += lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
                                   node_boundary_values[node_id] *