From 3a1cc197c9fcbe1f0bb9e473dbc528eaf52bf54b Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Fri, 23 Sep 2022 16:47:02 +0200
Subject: [PATCH] Beggining of Neumann BC for 3D
---
src/scheme/ScalarHybridScheme.cpp | 251 +++++++++++++++++++++++-------
1 file changed, 199 insertions(+), 52 deletions(-)
diff --git a/src/scheme/ScalarHybridScheme.cpp b/src/scheme/ScalarHybridScheme.cpp
index a9bafc3e2..e6fee0492 100644
--- a/src/scheme/ScalarHybridScheme.cpp
+++ b/src/scheme/ScalarHybridScheme.cpp
@@ -492,22 +492,22 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
return compute_node_is_corner;
}();
- const NodeValue<const bool> node_is_straight_border = [&] {
- NodeValue<bool> compute_node_is_straight_border{mesh->connectivity()};
- compute_node_is_straight_border.fill(false);
+ const NodeValue<const bool> node_is_edge = [&] {
+ NodeValue<bool> compute_node_is_edge{mesh->connectivity()};
+ compute_node_is_edge.fill(false);
parallel_for(
mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
if (is_boundary_node[node_id]) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
if ((Dimension == 3) && (node_to_cell.size() == 2)) {
- compute_node_is_straight_border[node_id] = true;
- //// std::cout << "node_is_straight_border : " << node_id << "; " <<
- // compute_node_is_straight_border[node_id]
+ compute_node_is_edge[node_id] = true;
+ //// std::cout << "node_is_edge : " << node_id << "; " <<
+ // compute_node_is_edge[node_id]
// << "\n";
}
}
});
- return compute_node_is_straight_border;
+ return compute_node_is_edge;
}();
// std::cout << "Test 1.9 \n";
@@ -1088,14 +1088,14 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
Cjr2[0] = -0.5 * (xrp2[1] - xrp1[1]);
Cjr2[1] = 0.5 * (xrp2[0] - xrp1[0]);
beta[node_id] = tensorProduct(Cjr1, (xr[node_id] - xj1)) + tensorProduct(Cjr2, (xr[node_id] - xj2));
- } else if (Dimension == 3 && not node_is_straight_border[node_id]) {
+ } else if (Dimension == 3 && not node_is_edge[node_id]) {
// std::cout << "Test 1.24.2 \n";
// std::cout << "coin : " << node_id << " " << Vj[cell_id] << "\n";
const TinyMatrix<Dimension> I = identity;
beta[node_id] = 0.125 * Vj[cell_id] * I;
// throw NotImplementedError("The scheme is not implemented in this dimension.");
- } else if (Dimension == 3 && node_is_straight_border[node_id]) {
+ } else if (Dimension == 3 && node_is_edge[node_id]) {
// std::cout << "Test 1.24.3 \n";
// std::cout << " bord droit : " << node_id << " " << Vj[cell_id] << "\n";
const TinyMatrix<Dimension> I = identity;
@@ -1246,52 +1246,123 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
}
}
- } else if ((node_is_neumann[node_id]) && (not node_is_corner[node_id]) && (not node_is_dirichlet[node_id])) {
- const TinyMatrix<Dimension> I = identity;
- const TinyVector<Dimension> n = exterior_normal[node_id];
- const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
- const TinyMatrix<Dimension> Q = I - nxn;
-
- 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];
+ } else if ((node_is_neumann[node_id]) && (not node_is_dirichlet[node_id])) {
+ if (not node_is_corner[node_id]) {
+ const TinyMatrix<Dimension> I = identity;
+ const TinyVector<Dimension> n = exterior_normal[node_id];
+ const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
+ const TinyMatrix<Dimension> Q = I - nxn;
- 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];
+ 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];
- S1(cell_id_j, cell_id_k) +=
- (1 - lambda) * dt * (1. - cn_coeff / 2.) *
- 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));
- S2(cell_id_j, cell_id_k) -=
- (1 - lambda) * dt * (cn_coeff / 2.) *
- 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));
+ 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];
+ S1(cell_id_j, cell_id_k) +=
+ (1 - lambda) * dt * (1. - cn_coeff / 2.) *
+ 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));
+ S2(cell_id_j, cell_id_k) -=
+ (1 - lambda) * dt * (cn_coeff / 2.) *
+ 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];
+ const size_t numberOfNodesPerFace = face_to_node.size();
+ if (not is_boundary_face[face_id]) {
+ 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));
+ }
+ }
+ }
+ }
+ }
+ }
+ } 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_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];
- const size_t numberOfNodesPerFace = face_to_node.size();
- if (not is_boundary_face[face_id]) {
- 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));
+ 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) {
+ const FaceId face_id = cell_to_face[i_face];
+ 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;
+ 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];
+ }
+ }
+ }
+ // ICI
+ 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));
+ }
}
}
}
@@ -1494,7 +1565,7 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
}
- } else if (node_is_corner[node_id]) {
+ } 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;
@@ -1504,6 +1575,82 @@ class ScalarHybridSchemeHandler::ScalarHybridScheme : public ScalarHybridSchemeH
}
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) {
+ const CellId cell_id = node_to_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) {
+ const FaceId face_id = cell_to_face[i_face];
+ 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;
+ 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];
+ }
+ }
+ }
+ 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));
+ }
+ }
+ }
}
} else if (node_is_dirichlet[node_id]) {
--
GitLab