diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index b650caa5440c0a2fa9768c9d3534dc8a9366cb22..84ec83e474941185aa06420b4e536f70c2c4f539 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -807,6 +807,153 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
return compute_w;
}();
+ const NodeValue<const TinyVector<2>> normal_edge_index = [&] {
+ NodeValue<TinyVector<2>> compute_normal_edge_index{mesh->connectivity()};
+ compute_normal_edge_index.fill(zero);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ bool l1_OK = false;
+
+ if (node_is_edge[node_id]) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ CellId cell_id = node_to_cell[0];
+ const auto& cell_to_node = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+ const NodeId node_id_in_cell = cell_to_node[i_node];
+ if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) {
+ cell_id = node_to_cell[1];
+ }
+ }
+
+ 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];
+ if (is_boundary_face[face_id]) {
+ if (not l1_OK) {
+ compute_normal_edge_index[node_id][0] = face_id;
+ l1_OK = true;
+ } else {
+ compute_normal_edge_index[node_id][1] = face_id;
+ }
+ }
+ }
+ }
+ });
+ return compute_normal_edge_index;
+ }();
+
+ const NodeValue<const TinyVector<Dimension>> normal1_edge = [&] {
+ NodeValue<TinyVector<Dimension>> compute_normal1_edge{mesh->connectivity()};
+ compute_normal1_edge.fill(zero);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (node_is_edge[node_id]) {
+ const FaceId l1 = normal_edge_index[node_id][0];
+
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ CellId cell_id = node_to_cell[0];
+ const auto& cell_to_node = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+ const NodeId node_id_in_cell = cell_to_node[i_node];
+ if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) {
+ cell_id = node_to_cell[1];
+ }
+ }
+
+ if (dot(nl[l1], xl[l1] - xj[cell_id]) <= 0) {
+ compute_normal1_edge[node_id] = -nl[l1];
+ } else {
+ compute_normal1_edge[node_id] = nl[l1];
+ }
+ }
+ });
+ return compute_normal1_edge;
+ }();
+
+ const NodeValue<const TinyVector<Dimension>> normal2_edge = [&] {
+ NodeValue<TinyVector<Dimension>> compute_normal2_edge{mesh->connectivity()};
+ compute_normal2_edge.fill(zero);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (node_is_edge[node_id]) {
+ const FaceId l2 = normal_edge_index[node_id][1];
+
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ CellId cell_id = node_to_cell[0];
+ const auto& cell_to_node = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+ const NodeId node_id_in_cell = cell_to_node[i_node];
+ if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) {
+ cell_id = node_to_cell[1];
+ }
+ }
+ if (dot(nl[l2], xl[l2] - xj[cell_id]) <= 0) {
+ compute_normal2_edge[node_id] = -nl[l2];
+ } else {
+ compute_normal2_edge[node_id] = nl[l2];
+ }
+ }
+ });
+ return compute_normal2_edge;
+ }();
+
+ const NodeValue<const TinyVector<Dimension>> normal_orth_edge = [&] {
+ NodeValue<TinyVector<Dimension>> compute_normal_orth_edge{mesh->connectivity()};
+ compute_normal_orth_edge.fill(zero);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ // bool l1_OK = false;
+ // FaceId l1 = 0;
+ // FaceId l2 = 0;
+ // TinyVector<Dimension> n1 = zero;
+ // TinyVector<Dimension> n2 = zero;
+
+ if (node_is_edge[node_id]) {
+ // const auto& node_to_cell = node_to_cell_matrix[node_id];
+ // CellId cell_id = node_to_cell[0];
+ // const auto& cell_to_node = cell_to_node_matrix[cell_id];
+ // for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+ // const NodeId node_id_in_cell = cell_to_node[i_node];
+ // if (node_is_corner[node_id_in_cell] && not node_is_edge[node_id_in_cell]) {
+ // cell_id = node_to_cell[1];
+ // }
+ // }
+
+ // 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];
+ // if (is_boundary_face[face_id]) {
+ // if (not l1_OK) {
+ // l1 = face_id;
+ // if (dot(nl[l1], xl[l1] - xj[cell_id]) <= 0) {
+ // n1 = -nl[l1];
+ // } else {
+ // n1 = nl[l1];
+ // }
+ // l1_OK = true;
+
+ // } else {
+ // l2 = face_id;
+ // if (dot(nl[l2], xl[l2] - xj[cell_id]) <= 0) {
+ // n2 = -nl[l2];
+ // } else {
+ // n2 = nl[l2];
+ // }
+ // }
+ // }
+ // }
+
+ compute_normal_orth_edge[node_id][0] =
+ normal1_edge[node_id][1] * normal2_edge[node_id][2] - normal1_edge[node_id][2] * normal2_edge[node_id][1];
+ compute_normal_orth_edge[node_id][1] =
+ normal1_edge[node_id][2] * normal2_edge[node_id][0] - normal1_edge[node_id][0] * normal2_edge[node_id][2];
+ compute_normal_orth_edge[node_id][2] =
+ normal1_edge[node_id][0] * normal2_edge[node_id][1] - normal1_edge[node_id][1] * normal2_edge[node_id][0];
+ }
+ });
+ return compute_normal_orth_edge;
+ }();
+
const FaceValue<const double> face_boundary_values = [&] {
FaceValue<double> compute_face_boundary_values{mesh->connectivity()};
FaceValue<double> sum_mes_l{mesh->connectivity()};
@@ -1273,9 +1420,11 @@ 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 CellId cell_id_k = node_to_cell[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) {
@@ -1283,59 +1432,27 @@ 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;
- 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];
- 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;
- }
- }
- }
- 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";
+ TinyVector<Dimension> t = normal_orth_edge[node_id];
- // 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) {
- // 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) += (1. - lambda) * dt * (1. - cn_coeff / 2.) *
+ 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,
+ node_kappar[node_id] * Cjr(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -= (1. - lambda) * dt * (cn_coeff / 2.) *
+ 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,
+ node_kappar[node_id] * Cjr(cell_id_j, i_node_j));
+
S1(cell_id_j, cell_id_k) += lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] /
numberOfNodesPerFace *
dot(1. /
@@ -1350,8 +1467,20 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
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) -= (1. - lambda) * dt * (1. - cn_coeff / 2.) *
+ 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,
+ node_kappar[node_id] * Cjr(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) += (1. - lambda) * dt * (cn_coeff / 2.) *
+ 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,
+ node_kappar[node_id] * Cjr(cell_id_j, i_node_j));
+
S1(cell_id_j, cell_id_k) -= lambda * dt * (1. - cn_coeff / 2.) * mes_l[face_id] /
numberOfNodesPerFace *
dot(1. /
@@ -1574,7 +1703,8 @@ 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];
+ // const CellId cell_id = node_to_cell[0];
+ const CellId cell_id = node_to_cell[0];
for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
FaceId face_id = node_to_face[i_face];
const auto& face_to_node = face_to_node_matrix[face_id];
@@ -1583,11 +1713,15 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
b_prev[cell_id] +=
(1 - lambda) / numberOfNodesPerFace * face_boundary_prev_values[face_id] * mes_l[face_id];
}
+
} else if (node_is_edge[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);
+ // 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) {
const CellId cell_id = node_to_cell[i_cell];
+ const size_t i_node = node_local_number_in_its_cell[i_cell];
+
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) {
@@ -1595,66 +1729,74 @@ 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 = 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];
- 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]) <= 0) {
- n1 = -nl[l1];
- } else {
- n1 = nl[l1];
- }
- if (dot(nl[l2], xl[l2] - xj[cell_id]) <= 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];
- }
+ FaceId l1 = normal_edge_index[node_id][0];
+ FaceId l2 = normal_edge_index[node_id][1];
+ TinyVector<Dimension> n1 = normal1_edge[node_id];
+ TinyVector<Dimension> n2 = normal2_edge[node_id];
+ TinyVector<Dimension> t = normal_orth_edge[node_id];
+
+ // std::cout << "node_id : " << node_id << "; cell_id : " << cell_id
+ // << "; xj1-xj2 = " << xj1_xj2(cell_id, i_face) << "; t = " << t << "\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) {
+ b[cell_id] += (1. - lambda) *
+ dot(face_boundary_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) *
+ node_kappar[node_id] * n1 +
+ face_boundary_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
+ node_kappar[node_id] * n2 -
+ 1. / dot(xj1_xj2(cell_id, i_face), node_kappar[node_id] * t) *
+ dot(xj1_xj2(cell_id, i_face),
+ face_boundary_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) *
+ node_kappar[node_id] * n1 +
+ face_boundary_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
+ node_kappar[node_id] * n2) /
+ pow(l2Norm(node_kappar[node_id] * t), 2) * node_kappar[node_id] * t,
+ node_kappar[node_id] * Cjr(cell_id, i_node));
+ b_prev[cell_id] +=
+ (1. - lambda) *
+ dot(face_boundary_prev_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) *
+ node_kappar[node_id] * n1 +
+ face_boundary_prev_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
+ node_kappar[node_id] * n2 -
+ 1. / dot(xj1_xj2(cell_id, i_face), node_kappar[node_id] * t) *
+ dot(xj1_xj2(cell_id, i_face),
+ face_boundary_prev_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) *
+ node_kappar[node_id] * n1 +
+ face_boundary_prev_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
+ node_kappar[node_id] * n2) /
+ pow(l2Norm(node_kappar[node_id] * t), 2) * node_kappar[node_id] * t,
+ node_kappar[node_id] * Cjr(cell_id, i_node));
+
+ b[cell_id] += lambda * mes_l[face_id] / numberOfNodesPerFace *
+ dot(face_boundary_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) *
+ node_kappar[node_id] * n1 +
+ face_boundary_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
+ node_kappar[node_id] * n2 -
+ 1. / dot(xj1_xj2(cell_id, i_face), node_kappar[node_id] * t) *
+ dot(xj1_xj2(cell_id, i_face),
+ face_boundary_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) *
+ node_kappar[node_id] * n1 +
+ face_boundary_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
+ node_kappar[node_id] * n2) /
+ pow(l2Norm(node_kappar[node_id] * t), 2) * node_kappar[node_id] * t,
+ base_xj1_xj2_orth(cell_id, i_face));
+ b_prev[cell_id] +=
+ lambda * mes_l[face_id] / numberOfNodesPerFace *
+ dot(face_boundary_prev_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) *
+ node_kappar[node_id] * n1 +
+ face_boundary_prev_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
+ node_kappar[node_id] * n2 -
+ 1. / dot(xj1_xj2(cell_id, i_face), node_kappar[node_id] * t) *
+ dot(xj1_xj2(cell_id, i_face),
+ face_boundary_prev_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) *
+ node_kappar[node_id] * n1 +
+ face_boundary_prev_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
+ node_kappar[node_id] * n2) /
+ pow(l2Norm(node_kappar[node_id] * t), 2) * node_kappar[node_id] * t,
+ base_xj1_xj2_orth(cell_id, i_face));
}
}
- b[cell_id] +=
- lambda * mes_l[face_id] / numberOfNodesPerFace *
- dot(face_boundary_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) * node_kappar[node_id] * n1 +
- face_boundary_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) * node_kappar[node_id] *
- n2 -
- 1. / dot(xj1_xj2(cell_id, i_face), node_kappar[node_id] * t) *
- dot(xj1_xj2(cell_id, i_face),
- face_boundary_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) *
- node_kappar[node_id] * n1 +
- face_boundary_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
- node_kappar[node_id] * n2) /
- pow(l2Norm(node_kappar[node_id] * t), 2) * node_kappar[node_id] * t,
- base_xj1_xj2_orth(cell_id, i_face));
- b_prev[cell_id] +=
- lambda * mes_l[face_id] / numberOfNodesPerFace *
- dot(face_boundary_prev_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) * node_kappar[node_id] *
- n1 +
- face_boundary_prev_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
- node_kappar[node_id] * n2 -
- 1. / dot(xj1_xj2(cell_id, i_face), node_kappar[node_id] * t) *
- dot(xj1_xj2(cell_id, i_face),
- face_boundary_prev_values[l1] / pow(l2Norm(node_kappar[node_id] * n1), 2) *
- node_kappar[node_id] * n1 +
- face_boundary_prev_values[l2] / pow(l2Norm(node_kappar[node_id] * n2), 2) *
- node_kappar[node_id] * n2) /
- pow(l2Norm(node_kappar[node_id] * t), 2) * node_kappar[node_id] * t,
- base_xj1_xj2_orth(cell_id, i_face));
}
}
}