Skip to content
Snippets Groups Projects
Commit 2bbdd33e authored by PATELA Julie's avatar PATELA Julie
Browse files

Fix Neumann boundary conditions in 2d

parent 33ecbe61
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/language/algorithms/HeatDiamondAlgorithm.cpp
+ 46
− 123
View file @ 2bbdd33e
......@@ -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,8 +441,9 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
SparseMatrixDescriptor S{number_of_dof};
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
if (not primal_face_is_neumann[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];
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];
......@@ -477,6 +457,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
}
}
}
FaceValue<const CellId> face_dual_cell_id = [=]() {
FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
......@@ -545,97 +526,25 @@ 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';
}
}
}
}
}
}
}
// // 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]) {
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 (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 (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(
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment