diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index be6552c2761c29f822e30265a6bd2ab7280de579..0a6058e73612670524dc9bfece6e2e5ff04432c8 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -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);
CellValue<double> Tj =
@@ -270,7 +286,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
}
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;
}
for (size_t i = 1; i < Dimension + 1; i++) {
@@ -302,13 +318,34 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
Tr[i_node] += x[j] * Tj[node_to_cell[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) {
FaceId face_id = node_to_face_matrix[i_node][i_face];
if (primal_face_is_on_boundary[face_id]) {
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(
return computed_alpha_l;
}();
+
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];
@@ -478,10 +516,29 @@ 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';
}
}
}
@@ -489,6 +546,105 @@ 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]) {
+ 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);
+ }
+ }
+ }
+ }
+ }
+ }
CellValue<double> fj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
@@ -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};
- 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};