diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index ecf6d84ed3bf4e76c150d325c923c024fc3e8365..187d452b65ef2d56d95c7c1d72e25c78df12d1d4 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -1077,14 +1077,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 (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), 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));
+ }
}
}
}
@@ -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] *