diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index e6fee049283821fda3c35a048b3c109966d70392..b650caa5440c0a2fa9768c9d3534dc8a9366cb22 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -457,35 +457,10 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
if (is_boundary_node[node_id]) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
- // const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
if (node_to_cell.size() == 1) {
compute_node_is_corner[node_id] = true;
- //// std::cout << "node is corner : " << node_id << "; " << compute_node_is_corner[node_id] << "\n";
- // for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
- // const unsigned int i_node = node_local_number_in_its_cell[i_cell];
- // const CellId cell_id = node_to_cell[i_cell];
- // const auto& cell_nodes = cell_to_node_matrix[cell_id];
- // NodeId i_node_prev;
- // NodeId i_node_next;
- // if (i_node == 0) {
- // i_node_prev = cell_nodes.size() - 1;
- // } else {
- // i_node_prev = i_node - 1;
- // }
- // if (i_node == cell_nodes.size() - 1) {
- // i_node_next = 0;
- // } else {
- // i_node_next = i_node + 1;
- // }
- // const NodeId prev_node_id = cell_to_node_matrix[cell_id][i_node_prev];
- // const NodeId next_node_id = cell_to_node_matrix[cell_id][i_node_next];
- // if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) {
- // compute_node_is_corner[node_id] = true;
- // }
- // }
} else if ((Dimension == 3) && (node_to_cell.size() == 2)) {
compute_node_is_corner[node_id] = true;
- //// std::cout << " node is corner droit " << node_id << "; " << compute_node_is_corner[node_id] << "\n";
}
}
});
@@ -501,9 +476,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
const auto& node_to_cell = node_to_cell_matrix[node_id];
if ((Dimension == 3) && (node_to_cell.size() == 2)) {
compute_node_is_edge[node_id] = true;
- //// std::cout << "node_is_edge : " << node_id << "; " <<
- // compute_node_is_edge[node_id]
- // << "\n";
}
}
});
@@ -514,7 +486,6 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
const NodeValue<const bool> node_is_angle = [&] {
NodeValue<bool> compute_node_is_angle{mesh->connectivity()};
compute_node_is_angle.fill(false);
- // for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
if (is_boundary_node[node_id] && Dimension == 2) {
@@ -536,7 +507,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
if (l2Norm(n1 - n2) > (1.E-15) and l2Norm(n1 + n2) > (1.E-15) and not node_is_corner[node_id]) {
compute_node_is_angle[node_id] = true;
}
- //// std::cout << node_id << " " << n1 << " " << n2 << " " << compute_node_is_angle[node_id] << "\n";
+ // std ::cout << node_id << " " << n1 << " " << n2 << " " << compute_node_is_angle[node_id] << "\n";
}
});
return compute_node_is_angle;
@@ -1302,11 +1273,9 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
} else if (node_is_edge[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];
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) {
@@ -1314,54 +1283,92 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
const auto& face_to_node = face_to_node_matrix[face_id];
const size_t numberOfNodesPerFace = face_to_node.size();
if (not is_boundary_face[face_id]) {
- bool l1_OK = false;
- FaceId l1;
- FaceId l2;
- TinyVector<Dimension> n1;
- TinyVector<Dimension> n2;
- TinyVector<Dimension> t;
+ bool l1_OK = false;
+ bool l2_OK = false;
+ FaceId l1 = 0;
+ FaceId l2 = 0;
+ TinyVector<Dimension> n1 = zero;
+ TinyVector<Dimension> n2 = zero;
+ TinyVector<Dimension> t = zero;
+ // std::cout << "Test 1.3.01 \n";
for (size_t j_face = 0; j_face < cell_to_face.size(); ++j_face) {
- const FaceId face_id_j = cell_to_face[j_face];
- const auto& face_to_node_j = face_to_node_matrix[face_id_j];
- for (size_t j_face_to_node = 0; j_face_to_node < face_to_node_j.size(); ++j_face_to_node) {
- if (is_boundary_face[face_id_j]) {
- if (not l1_OK) {
- l1 = face_id_j;
- l1_OK = true;
- } else {
- l2 = face_id_j;
- }
- if (dot(nl[l1], xl[l1] - xj[cell_id_j]) <= 0) {
- n1 = -nl[l1];
- } else {
- n1 = nl[l1];
- }
- if (dot(nl[l2], xl[l2] - xj[cell_id_j]) <= 0) {
- n2 = -nl[l2];
- } else {
- n2 = nl[l2];
- }
- t[0] = n1[1] * n2[2] - n1[2] * n2[1];
- t[1] = n1[2] * n2[0] - n1[0] * n2[2];
- t[2] = n1[0] * n2[1] - n1[1] * n2[0];
+ const FaceId face_id_j = cell_to_face[j_face];
+ if (is_boundary_face[face_id_j]) {
+ // const auto& face_to_node_j = face_to_node_matrix[face_id_j];
+ // std::cout << "Test 1.3.02 \n";
+
+ // std::cout << "Test 1.3.03 \n";
+ if (not l1_OK) {
+ l1 = face_id_j;
+ l1_OK = true;
+ } else {
+ l2 = face_id_j;
+ l2_OK = true;
}
}
}
- // ICI
+ if (dot(nl[l1], xl[l1] - xj[cell_id_j]) <= 0 and l1_OK) {
+ n1 = -nl[l1];
+ } else {
+ n1 = nl[l1];
+ }
+ // std::cout << "Test 1.3.04 \n";
+ if (dot(nl[l2], xl[l2] - xj[cell_id_j]) <= 0 and l2_OK) {
+ n2 = -nl[l2];
+ } else {
+ n2 = nl[l2];
+ }
+ // std::cout << "Test 1.3.05 \n";
+ if (l1_OK && l2_OK) {
+ t[0] = n1[1] * n2[2] - n1[2] * n2[1];
+ t[1] = n1[2] * n2[0] - n1[0] * n2[2];
+ t[2] = n1[0] * n2[1] - n1[1] * n2[0];
+ }
+ // std::cout << "Test 1.3.06 \n";
+
+ // std::cout << "Test 1.3.1 \n";
for (size_t i_face_to_node = 0; i_face_to_node < face_to_node.size(); ++i_face_to_node) {
if (face_to_node[i_face_to_node] == node_id) {
- S1(cell_id_j, cell_id_k) +=
- lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] / numberOfNodesPerFace *
- 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]),
- base_xj1_xj2_orth(cell_id_j, i_face));
- S2(cell_id_j, cell_id_k) -=
- lambda * dt * (cn_coeff / 2.) * mes_l[face_id] / numberOfNodesPerFace *
- 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]),
- base_xj1_xj2_orth(cell_id_j, i_face));
+ // std::cout << "Test 1.3.1.1 \n";
+ const auto& face_to_cell = face_to_cell_matrix[face_id];
+ for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+ // std::cout << "Test 1.3.1.2 \n";
+ CellId cell_id = face_to_cell[i_cell];
+ if (cell_id_k == cell_id_j) {
+ // std::cout << "Test 1.3.1.3 \n";
+ S1(cell_id_j, cell_id_k) += lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] /
+ numberOfNodesPerFace *
+ dot(1. /
+ (dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t) *
+ pow(l2Norm(node_kappar[node_id] * t), 2)) *
+ node_kappar[node_id] * t,
+ base_xj1_xj2_orth(cell_id_j, i_face));
+ S2(cell_id_j, cell_id_k) -= lambda * dt * (cn_coeff / 2.) * mes_l[face_id] /
+ numberOfNodesPerFace *
+ dot(1. /
+ (dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t) *
+ pow(l2Norm(node_kappar[node_id] * t), 2)) *
+ node_kappar[node_id] * t,
+ base_xj1_xj2_orth(cell_id_j, i_face));
+ // std::cout << "Test 1.3.2 \n";
+ } else {
+ S1(cell_id_j, cell_id_k) -= lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] /
+ numberOfNodesPerFace *
+ dot(1. /
+ (dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t) *
+ pow(l2Norm(node_kappar[node_id] * t), 2)) *
+ node_kappar[node_id] * t,
+ base_xj1_xj2_orth(cell_id_j, i_face));
+ S2(cell_id_j, cell_id_k) += lambda * dt * (cn_coeff / 2.) * mes_l[face_id] /
+ numberOfNodesPerFace *
+ dot(1. /
+ (dot(xj1_xj2(cell_id_j, i_face), node_kappar[node_id] * t) *
+ pow(l2Norm(node_kappar[node_id] * t), 2)) *
+ node_kappar[node_id] * t,
+ base_xj1_xj2_orth(cell_id_j, i_face));
+ // std::cout << "Test 1.3.3 \n";
+ }
+ }
}
}
}
@@ -1484,7 +1491,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
}
};
-
+ // std::cout << "Test 1.31" << "\n";
//// std::cout << "S1 = " << S1 << "\n";
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
@@ -1568,14 +1575,14 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
} else if (node_is_corner[node_id] && not node_is_edge[node_id]) {
const auto& node_to_face = node_to_face_matrix[node_id];
const CellId cell_id = node_to_cell[0];
- double sum_mes_l = 0;
for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
- FaceId face_id = node_to_face[i_face];
- sum_mes_l += mes_l[face_id];
+ FaceId face_id = node_to_face[i_face];
+ const auto& face_to_node = face_to_node_matrix[face_id];
+ const size_t numberOfNodesPerFace = face_to_node.size();
+ b[cell_id] += (1 - lambda) / numberOfNodesPerFace * face_boundary_values[face_id] * mes_l[face_id];
+ b_prev[cell_id] +=
+ (1 - lambda) / numberOfNodesPerFace * face_boundary_prev_values[face_id] * mes_l[face_id];
}
- b[cell_id] += (1 - lambda) * 0.5 * node_boundary_values[node_id] * sum_mes_l;
- b_prev[cell_id] += (1 - lambda) * 0.5 * node_boundary_prev_values[node_id] * sum_mes_l;
-
} else if (node_is_edge[node_id]) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
@@ -1588,12 +1595,12 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
const auto& face_to_node = face_to_node_matrix[face_id];
const size_t numberOfNodesPerFace = face_to_node.size();
if (not is_boundary_face[face_id]) {
- bool l1_OK = false;
- FaceId l1;
- FaceId l2;
- TinyVector<Dimension> n1;
- TinyVector<Dimension> n2;
- TinyVector<Dimension> t;
+ bool l1_OK = false;
+ FaceId l1 = 0;
+ FaceId l2 = 0;
+ TinyVector<Dimension> n1 = zero;
+ TinyVector<Dimension> n2 = zero;
+ TinyVector<Dimension> t = zero;
for (size_t j_face = 0; j_face < cell_to_face.size(); ++j_face) {
const FaceId face_id_j = cell_to_face[j_face];
const auto& face_to_node_j = face_to_node_matrix[face_id_j];