diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index 10d41c948a8e3ccacb6b3acca3e28b7e61f96b78..55c38578f5a025327d59b4f380f56dc9583adb38 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -694,6 +694,8 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
compute_normal_base(cell_id, i_face) = inverse(A) * (face_kappal[face_id] * face_normal);
+ // std::cout << "alpha_l = " << compute_normal_base(cell_id, i_face)[0]
+ // << "; delta_l = " << compute_normal_base(cell_id, i_face)[1] << "\n";
}
}
@@ -1076,6 +1078,8 @@ 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_node_j));
+ // std::cout << "alpha_l = " << normal_face_base(cell_id_j, i_face)[0] << "\n";
+ // std::cout << S1(cell_id_j, cell_id_k) << "\n";
}
}
}
@@ -1190,14 +1194,14 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
const CellId cell_id_k = face_to_cell[i_cell_k];
if (cell_id_j == cell_id_k) {
- S1(cell_id_j, cell_id_k) -=
+ S1(cell_id_j, cell_id_k) +=
lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
- S2(cell_id_j, cell_id_k) +=
+ S2(cell_id_j, cell_id_k) -=
lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
} else {
- S1(cell_id_j, cell_id_k) +=
+ S1(cell_id_j, cell_id_k) -=
lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
- S2(cell_id_j, cell_id_k) -=
+ S2(cell_id_j, cell_id_k) +=
lambda * dt * (cn_coeff / 2.) * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[1];
}
}
@@ -1214,6 +1218,8 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
};
+ // std::cout << "S1 = " << S1 << "\n";
+
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
S1(cell_id, cell_id) += Vj[cell_id];
S2(cell_id, cell_id) += Vj[cell_id];
@@ -1229,13 +1235,11 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
Ph[cell_id] = Pj[cell_id];
};
- std::cout << "P = " << Ph << "\n";
+ // std::cout << "P = " << Ph << "\n";
CRSMatrix A1{S1.getCRSMatrix()};
CRSMatrix A2{S2.getCRSMatrix()};
- // std::cout << A1 << "\n";
-
Vector<double> b{mesh->numberOfCells()};
b = zero;
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
@@ -1304,15 +1308,18 @@ 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) {
- b[cell_id_j] += (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
- face_boundary_values[face_id] *
- dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
- b_prev[cell_id_j] +=
- (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
- face_boundary_prev_values[face_id] *
- dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+ if (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] *
+ dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+ // std::cout << "xl - xj2 = " << xl[face_id] -
+ b_prev[cell_id_j] +=
+ lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ node_boundary_prev_values[node_id] *
+ dot(inverse(node_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+ }
}
}
}
@@ -1338,16 +1345,17 @@ 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) {
- b[cell_id_j] +=
- (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
- face_boundary_values[face_id] *
- dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
- b_prev[cell_id_j] +=
- (1 - lambda) * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
- face_boundary_prev_values[face_id] *
- dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+ if (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] *
+ dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+ b_prev[cell_id_j] +=
+ lambda * 0.5 * mes_l[face_id] * normal_face_base(cell_id_j, i_face)[0] *
+ node_boundary_prev_values[node_id] *
+ dot(inverse(corner_betar[node_id]) * sum_Cjr[node_id], xj1_xj2_orth(cell_id_j, i_node_j));
+ }
}
}
}
@@ -1390,9 +1398,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
B[cell_id] = dt * ((1. - cn_coeff / 2.) * b[cell_id] + cn_coeff / 2. * b_prev[cell_id]) + A2P[cell_id];
});
- std::cout << "B = " << B << "\n"
- << "A2P = " << A2P << "\n";
-
NodeValue<TinyVector<Dimension>> ur{mesh->connectivity()};
ur.fill(zero);
for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {