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

Try to fix Neumann

parent c612d4c0
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
+ 216
− 3
View file @ 447834c0
...@@ -211,6 +211,22 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( ...@@ -211,6 +211,22 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
} }
} }
FaceValue<bool> primal_face_is_neumann(mesh->connectivity());
if (parallel::size() > 1) {
throw NotImplementedError("Calculation of face_is_neumann is incorrect");
}
primal_face_is_neumann.fill(false);
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
if (primal_face_is_on_boundary[face_id]) {
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 (not is_dirichlet[node_id]) {
primal_face_is_neumann[face_id] = true;
}
}
}
}
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
CellValue<double> Tj = CellValue<double> Tj =
...@@ -270,7 +286,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( ...@@ -270,7 +286,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
} }
} }
LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used}; LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
for (size_t j = 0; j < node_to_cell.size(); j++) { for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) {
A(0, j) = 1; A(0, j) = 1;
} }
for (size_t i = 1; i < Dimension + 1; i++) { for (size_t i = 1; i < Dimension + 1; i++) {
...@@ -302,13 +318,34 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( ...@@ -302,13 +318,34 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
Tr[i_node] += x[j] * Tj[node_to_cell[j]]; Tr[i_node] += x[j] * Tj[node_to_cell[j]];
w_rj(i_node, j) = x[j]; w_rj(i_node, j) = x[j];
} }
int cpt_face = 0; int cpt_face = node_to_cell.size();
for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) { for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
FaceId face_id = node_to_face_matrix[i_node][i_face]; FaceId face_id = node_to_face_matrix[i_node][i_face];
if (primal_face_is_on_boundary[face_id]) { if (primal_face_is_on_boundary[face_id]) {
w_rl(i_node, i_face) = x[cpt_face++]; 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';
}
}
} }
} }
...@@ -393,6 +430,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( ...@@ -393,6 +430,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
return computed_alpha_l; return computed_alpha_l;
}(); }();
SparseMatrixDescriptor S{number_of_dof}; SparseMatrixDescriptor S{number_of_dof};
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
const auto& primal_face_to_cell = face_to_cell_matrix[face_id]; const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
...@@ -478,12 +516,130 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( ...@@ -478,12 +516,130 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
const double a_ir = alpha_face_id * (nil, Clr); 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) { 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]; 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; 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]) {
int nb_node_per_face = primal_face_to_node_matrix[face_id].size();
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 (not is_dirichlet[node_id]) {
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(face_dof_number[face_id], cell_dof_number[j_id]) -= (1. / nb_node_per_face) * w_rj(node_id, j_cell);
}
}
}
}
}
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
if (primal_face_is_neumann[face_id]) {
S(face_dof_number[face_id], face_dof_number[face_id]) += 1;
}
if (primal_face_is_neumann[face_id]) {
int nb_node_per_face = primal_face_to_node_matrix[face_id].size();
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 (not is_dirichlet[node_id]) {
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);
}
} }
} }
} }
...@@ -572,8 +728,65 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( ...@@ -572,8 +728,65 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
} }
} }
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 (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
FaceId face_id = face_list[i_face];
CellId cell_id = face_to_cell_matrix[face_id][0];
Assert(face_to_cell_matrix[face_id].size() == 1);
b[cell_dof_number[cell_id]] += mes_l[face_id] * value_list[i_face];
}
}
},
boundary_condition);
}
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, DirichletBoundaryCondition>) {
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
NodeId node_id = node_list[i_node];
for (size_t i_face = 0; i_face < node_to_face_matrix[node_id].size(); ++i_face) {
FaceId face_id = node_to_face_matrix[node_id][i_face];
if (primal_face_is_neumann[face_id]) {
int nb_node_per_face = primal_face_to_node_matrix[face_id].size();
b[face_dof_number[face_id]] += (1. / nb_node_per_face) * value_list[i_node];
}
}
}
}
},
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';
// }
// }
Vector<double> T{number_of_dof}; Vector<double> T{number_of_dof};
T = 0; 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';
}
BiCGStab{b, A, T, 1000, 1e-15}; BiCGStab{b, A, T, 1000, 1e-15};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment