diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index 61285fab57d4cf8d3fb4b8f38c0d6cb117eab015..ecf6d84ed3bf4e76c150d325c923c024fc3e8365 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -1085,16 +1085,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
                 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), xj1_xj2_orth(cell_id_j, i_face));
-                // if (normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //       normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) !=
-                //     nl[face_id]) {
-                //   std::cout << "(xj1,xj2) = (" << xj[cell_id_k] << "," << xj[cell_id_j] << ") \n"
-                //             << normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //                  normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face)
-                //             << " = " << nl[face_id] << "\n"
-                //             << "\n";
-                // }
-                // std::cout << S1(cell_id_j, cell_id_k) << "\n";
               }
             }
           }
@@ -1146,23 +1136,15 @@ 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_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), 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), xj1_xj2_orth(cell_id_j, i_face));
-                // if (normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //       normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) !=
-                //     nl[face_id]) {
-                //   std::cout << "(xj1,xj2) = (" << xj[cell_id_k] << "," << xj[cell_id_j] << ") \n"
-                //             << normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //                  normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face)
-                //             << " = " << nl[face_id] << "\n"
-                //             << "\n";
-                // }
+              if (face_is_dirichlet[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), 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), xj1_xj2_orth(cell_id_j, i_face));
+                }
               }
             }
           }
@@ -1188,25 +1170,16 @@ 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_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(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
 
-              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(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), 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(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
-
-                // if (normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //       normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face) !=
-                //     nl[face_id]) {
-                //   std::cout << "(xj1,xj2) = (" << xj[cell_id_k] << "," << xj[cell_id_j] << ") \n"
-                //             << normal_face_base(cell_id_j, i_face)[0] * xj1_xj2_orth(cell_id_j, i_face) +
-                //                  normal_face_base(cell_id_j, i_face)[1] * xj1_xj2(cell_id_j, i_face)
-                //             << " = " << nl[face_id] << "\n"
-                //             << "\n";
-                // }
+                  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(corner_betar[node_id]) * Cjr(cell_id_k, i_node_k), xj1_xj2_orth(cell_id_j, i_face));
+                }
               }
             }
           }
@@ -1340,7 +1313,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 (is_boundary_face[face_id]) {
+              if (face_is_dirichlet[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] *
@@ -1376,7 +1349,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 (is_boundary_face[face_id]) {
+              if (face_is_dirichlet[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] *