From 2bbdd33e35979d902b26884ee43c3402afcee130 Mon Sep 17 00:00:00 2001
From: Julie PATELA <julie.patela.ocre@cea.fr>
Date: Fri, 7 Aug 2020 16:35:18 +0200
Subject: [PATCH] Fix Neumann boundary conditions in 2d
---
.../algorithms/HeatDiamondAlgorithm.cpp | 169 +++++-------------
1 file changed, 46 insertions(+), 123 deletions(-)
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index 04df2737b..74382e9be 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -354,27 +354,6 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
w_rl(i_node, i_face) = x[cpt_face++];
}
}
- std::cout << i_node << ":";
- double s = 0;
- for (size_t i = 0; i < x.size(); ++i) {
- std::cout << ' ' << x[i];
- s += x[i];
- }
- std::cout << " -> " << s << '\n';
- }
- }
-
- for (NodeId i_node = 0; i_node < mesh->numberOfNodes(); ++i_node) {
- if (not is_dirichlet[i_node]) {
- std::cout << "node = " << i_node << '\n';
- for (size_t i_cell = 0; i_cell < node_to_cell_matrix[i_node].size(); ++i_cell) {
- std::cout << "w_rj(" << i_node << ',' << i_cell << ") = " << w_rj(i_node, i_cell) << '\n';
- }
- if (primal_node_is_on_boundary[i_node]) {
- for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
- std::cout << "w_rl(" << i_node << ',' << i_face << ") = " << w_rl(i_node, i_face) << '\n';
- }
- }
}
}
@@ -462,17 +441,19 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
SparseMatrixDescriptor S{number_of_dof};
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
- const double beta_l = 0.5 * alpha_l[face_id] * mes_l[face_id];
-
- for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
- const CellId cell_id1 = primal_face_to_cell[i_cell];
- for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
- const CellId cell_id2 = primal_face_to_cell[j_cell];
- if (i_cell == j_cell) {
- S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
- } else {
- S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
+ if (not primal_face_is_neumann[face_id]) {
+ const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+ const double beta_l = 1. / Dimension * alpha_l[face_id] * mes_l[face_id];
+
+ for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
+ const CellId cell_id1 = primal_face_to_cell[i_cell];
+ for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
+ const CellId cell_id2 = primal_face_to_cell[j_cell];
+ if (i_cell == j_cell) {
+ S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
+ } else {
+ S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
+ }
}
}
}
@@ -545,29 +526,18 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
const double a_ir = alpha_face_id * (nil, Clr);
- double sum_w = 0;
-
for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
S(cell_dof_number[i_id], cell_dof_number[j_id]) -= w_rj(node_id, j_cell) * a_ir;
- sum_w += w_rj(node_id, j_cell);
}
- std::cout << node_id << " ";
if (primal_node_is_on_boundary[node_id]) {
- std::cout << "face: " << sum_w << " -> ";
for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
FaceId l_id = node_to_face_matrix[node_id][l_face];
if (primal_face_is_neumann[l_id]) {
S(cell_dof_number[i_id], face_dof_number[l_id]) -= w_rl(node_id, l_face) * a_ir;
- sum_w += w_rl(node_id, l_face);
- std::cout << "[added value :" << w_rl(node_id, l_face) * a_ir << "] Clr=" << Clr
- << " nil=" << nil << " ";
}
}
- } else {
- std::cout << "inner: ";
}
- std::cout << "sum_w=" << sum_w << '\n';
}
}
}
@@ -575,67 +545,6 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
}
}
- // // test
- // const auto& cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
- // for (const auto& boundary_condition : boundary_condition_list) {
- // std::visit(
- // [&](auto&& bc) {
- // using T = std::decay_t<decltype(bc)>;
- // if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
- // (std::is_same_v<T, FourierBoundaryCondition>)) {
- // const auto& face_list = bc.faceList();
- // // const auto& value_list = bc.valueList();
- // for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- // for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
- // FaceId face_id = face_list[i_face];
- // for (size_t i_face2 = 0; i_face2 < cell_to_face_matrix[cell_id].size(); ++i_face2) {
- // FaceId face_id2 = cell_to_face_matrix[cell_id][i_face2];
- // if (not primal_face_is_on_boundary[face_id2]) {
- // const double alpha_face_id = alpha_l[face_id2];
-
- // const bool is_face_reversed_for_cell_i = primal_face_cell_is_reversed(cell_id, i_face2);
- // for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id2].size(); ++i_node) {
- // NodeId node_id = primal_face_to_node_matrix[face_id2][i_node];
- // const TinyVector<Dimension> nil = [&] {
- // if (is_face_reversed_for_cell_i) {
- // return -primal_nlr(face_id2, i_node);
- // } else {
- // return primal_nlr(face_id2, i_node);
- // }
- // }();
- // if (not is_dirichlet[node_id]) {
- // if (primal_node_is_on_boundary[node_id]) {
- // CellId dual_cell_id = face_dual_cell_id[face_id2];
-
- // for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
- // ++i_dual_node) {
- // const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
- // if (dual_node_primal_node_id[dual_node_id] == node_id) {
- // const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-
- // const double a_ir = alpha_face_id * (nil, Clr);
-
- // for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
- // FaceId face_to_node_id = node_to_face_matrix[node_id][l_face];
- // if (face_to_node_id == face_id) {
- // S(cell_dof_number[cell_id], face_dof_number[face_id]) -=
- // w_rl(node_id, l_face) * a_ir;
- // }
- // }
- // }
- // }
- // }
- // }
- // }
- // }
- // }
- // }
- // }
- // }
- // },
- // boundary_condition);
- // }
- // std::exit(0);
// Termes pour Neumann et Fourier
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
if (primal_face_is_neumann[face_id]) {
@@ -664,8 +573,6 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
FaceId face_id2 = node_to_face_matrix[node_id][l_face];
if (primal_face_is_neumann[face_id2]) {
- std::cout << "TEST***********************"
- << "\n";
S(face_dof_number[face_id], face_dof_number[face_id2]) -=
(1. / nb_node_per_face) * w_rl(node_id, l_face);
}
@@ -797,28 +704,44 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
boundary_condition);
}
- std::cout << S << '\n';
-
- // for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- // if (primal_face_is_neumann[face_id]) {
- // b[face_dof_number[face_id]] = 100;
- // std::cout << "b[" << face_dof_number[face_id] << "] = " << b[face_dof_number[face_id]] << '\n';
- // }
- // }
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (primal_face_is_neumann[face_id]) {
+ CellId cell_id = face_to_cell_matrix[face_id][0];
+ const bool is_face_reversed_for_cell_i =
+ primal_face_cell_is_reversed(cell_id, face_local_numbers_in_their_cells(face_id, 0));
+ for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
+ if (is_dirichlet[node_id]) {
+ const TinyVector<Dimension> nil = [&] {
+ if (is_face_reversed_for_cell_i) {
+ return -primal_nlr(face_id, i_node);
+ } else {
+ return primal_nlr(face_id, i_node);
+ }
+ }();
+ CellId dual_cell_id = face_dual_cell_id[face_id];
+
+ for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
+ const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
+ if (dual_node_primal_node_id[dual_node_id] == node_id) {
+ const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
+ b[cell_dof_number[cell_id]] -= alpha_l[face_id] * (nil, Clr) * dirichlet_value[node_id]; ///
+ }
+ }
+ }
+ }
+ }
+ }
Vector<double> T{number_of_dof};
- T = 1;
- T[6] = 2;
- T[7] = 3;
- T[8] = 4;
- Vector AT = A * T;
- for (size_t i = 0; i < T.size(); ++i) {
- std::cout << " AT[" << i << "] = " << AT[i] << '\t';
- std::cout << " b[" << i << "] = " << b[i] << '\n';
- }
+ T = 0;
BiCGStab{b, A, T, 1000, 1e-15};
+ Vector r = A * T - b;
+
+ std::cout << "real residu = " << std::sqrt((r, r)) << '\n';
+
CellValue<double> Temperature{mesh->connectivity()};
parallel_for(
--
GitLab