From 1fb103d45986c9c5e5cb0b7579bbfdedd7f73045 Mon Sep 17 00:00:00 2001
From: LABOURASSE Emmanuel <labourassee@gmail.com>
Date: Tue, 6 Oct 2020 11:24:46 +0200
Subject: [PATCH] Correction of the Neumann boundary conditions for heat
---
.../algorithms/HeatDiamondAlgorithm.cpp | 449 ++++++------------
1 file changed, 156 insertions(+), 293 deletions(-)
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index c729f97b2..33d6783e5 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -52,85 +52,40 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
- NodeValue<bool> is_dirichlet{mesh->connectivity()};
- is_dirichlet.fill(false);
- NodeValue<double> dirichlet_value{mesh->connectivity()};
- dirichlet_value.fill(std::numeric_limits<double>::signaling_NaN());
-
for (const auto& bc_descriptor : bc_descriptor_list) {
bool is_valid_boundary_condition = true;
switch (bc_descriptor->type()) {
case IBoundaryConditionDescriptor::Type::symmetry: {
- // const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
- // dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
throw NotImplementedError("NIY");
break;
}
case IBoundaryConditionDescriptor::Type::dirichlet: {
const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
-
- for (size_t i_ref_face_list = 0; i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>();
- ++i_ref_face_list) {
- const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
- const RefId& ref = ref_face_list.refId();
- if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
- MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
-
- const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
-
- const auto& node_list = mesh_node_boundary.nodeList();
-
- Array<const double> value_list =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
- node_list);
-
- for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- NodeId node_id = node_list[i_node];
- if (not is_dirichlet[node_id]) {
- is_dirichlet[node_id] = true;
- dirichlet_value[node_id] = value_list[i_node];
- }
- }
-
- boundary_condition_list.push_back(DirichletBoundaryCondition{node_list, value_list});
- }
- }
-
if constexpr (Dimension > 1) {
- for (size_t i_ref_node_list = 0;
- i_ref_node_list < mesh->connectivity().template numberOfRefItemList<ItemType::node>(); ++i_ref_node_list) {
- const auto& ref_node_list = mesh->connectivity().template refItemList<ItemType::node>(i_ref_node_list);
- const RefId& ref = ref_node_list.refId();
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
- const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+ MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
- const auto& node_list = ref_node_list.list();
+ const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+ Array<const FaceId> face_list = ref_face_list.list();
Array<const double> value_list =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id,
- mesh->xr(),
- node_list);
-
- for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- NodeId node_id = node_list[i_node];
- if (not is_dirichlet[node_id]) {
- is_dirichlet[node_id] = true;
- dirichlet_value[node_id] = value_list[i_node];
- }
- }
-
- boundary_condition_list.push_back(DirichletBoundaryCondition{node_list, value_list});
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
+ mesh_data.xl(),
+ face_list);
+ boundary_condition_list.push_back(DirichletBoundaryCondition{face_list, value_list});
}
}
}
-
break;
}
case IBoundaryConditionDescriptor::Type::fourier: {
- // const FourierBoundaryConditionDescriptor& fourier_bc_descriptor =
- // dynamic_cast<const FourierBoundaryConditionDescriptor&>(*bc_descriptor);
throw NotImplementedError("NIY");
break;
}
@@ -154,8 +109,6 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
mesh_data.xl(),
face_list);
boundary_condition_list.push_back(NeumannBoundaryCondition{face_list, value_list});
- } else {
- throw NotImplementedError("Neumann conditions are not supported in 1d");
}
}
}
@@ -190,7 +143,8 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
using T = std::decay_t<decltype(bc)>;
if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
(std::is_same_v<T, FourierBoundaryCondition>) or
- (std::is_same_v<T, SymmetryBoundaryCondition>)) {
+ (std::is_same_v<T, SymmetryBoundaryCondition>) or
+ (std::is_same_v<T, DirichletBoundaryCondition>)) {
const auto& face_list = bc.faceList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
@@ -211,6 +165,30 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
return compute_face_dof_number;
}();
+ const FaceValue<const bool> primal_face_is_neumann = [&] {
+ FaceValue<bool> face_is_neumann{mesh->connectivity()};
+ face_is_neumann.fill(false);
+ 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>) or
+ (std::is_same_v<T, SymmetryBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ face_is_neumann[face_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return face_is_neumann;
+ }();
+
const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
@@ -240,27 +218,22 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
}
- FaceValue<bool> primal_face_is_neumann(mesh->connectivity());
+ FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
if (parallel::size() > 1) {
throw NotImplementedError("Calculation of face_is_neumann is incorrect");
}
- primal_face_is_neumann.fill(false);
+ primal_face_is_dirichlet.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;
- }
- }
- }
+ primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]));
}
+
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
CellValue<double> Tj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
-
+ FaceValue<double> Tl =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
NodeValue<double> Tr(mesh->connectivity());
const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
@@ -352,6 +325,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
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++];
+ Tr[i_node] += w_rl(i_node, i_face) * Tl[face_id];
}
}
}
@@ -441,32 +415,21 @@ 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 = 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;
+ 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;
+ if (primal_face_is_neumann[face_id]) {
+ S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= beta_l;
}
+ } else {
+ S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
}
}
- } else {
- 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];
- Assert(primal_face_to_cell.size() == 1, "Wrong number of cells for a boundary face");
- 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];
- S(cell_dof_number[cell_id1], cell_dof_number[cell_id1]) += beta_l;
- S(cell_dof_number[cell_id1], face_dof_number[face_id]) -= beta_l;
- S(face_dof_number[face_id], face_dof_number[face_id]) += beta_l;
- S(face_dof_number[face_id], cell_dof_number[cell_id1]) -= beta_l;
- }
}
}
@@ -505,48 +468,48 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
const auto& primal_face_cell_is_reversed = mesh->connectivity().cellFaceIsReversed();
const auto& face_local_numbers_in_their_cells = mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
-
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- if (face_to_cell_matrix[face_id].size() > 1) {
- const double alpha_face_id = alpha_l[face_id];
-
- for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
- CellId i_id = face_to_cell_matrix[face_id][i_face_cell];
- const bool is_face_reversed_for_cell_i =
- primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_face_cell));
-
- 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]) {
- 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);
- }
- }();
+ const double alpha_face_id = alpha_l[face_id];
- CellId dual_cell_id = face_dual_cell_id[face_id];
+ for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
+ CellId i_id = face_to_cell_matrix[face_id][i_face_cell];
+ const bool is_face_reversed_for_cell_i =
+ primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_face_cell));
- 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);
+ 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];
- const double a_ir = alpha_face_id * (nil, Clr);
+ 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);
+ }
+ }();
- 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;
- }
- if (primal_node_is_on_boundary[node_id]) {
- 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;
- }
+ 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);
+
+ const double a_ir = alpha_face_id * (nil, Clr);
+
+ 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;
+ if (primal_face_is_neumann[face_id]) {
+ S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * a_ir;
+ }
+ }
+ if (primal_node_is_on_boundary[node_id]) {
+ 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_on_boundary[l_id]) {
+ S(cell_dof_number[i_id], face_dof_number[l_id]) -= w_rl(node_id, l_face) * a_ir;
+ if (primal_face_is_neumann[face_id]) {
+ S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * a_ir;
}
}
}
@@ -556,125 +519,41 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
}
}
- // 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]) {
+ if (primal_face_is_dirichlet[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]) {
- 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());
const CellValue<const double> primal_Vj = mesh_data.Vj();
-
CRSMatrix A{S};
Vector<double> b{number_of_dof};
- const auto& primal_cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
b[cell_dof_number[cell_id]] = fj[cell_id] * primal_Vj[cell_id];
}
-
- { // Dirichlet
- NodeValue<bool> node_tag{mesh->connectivity()};
- node_tag.fill(false);
+ // Dirichlet on b^L_D
+ {
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& face_list = bc.faceList();
const auto& value_list = bc.valueList();
- for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- const NodeId node_id = node_list[i_node];
- if (not node_tag[node_id]) {
- node_tag[node_id] = true;
-
- const auto& node_cells = primal_node_to_cell_matrix[node_id];
- for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
- const CellId cell_id = node_cells[i_cell];
-
- const auto& primal_cell_to_face = primal_cell_to_face_matrix[cell_id];
-
- for (size_t i_cell_face = 0; i_cell_face < primal_cell_to_face.size(); ++i_cell_face) {
- FaceId face_id = primal_cell_to_face_matrix[cell_id][i_cell_face];
-
- const bool is_face_reversed_for_cell_i = [=] {
- for (size_t i_face = 0; i_face < primal_cell_to_face.size(); ++i_face) {
- FaceId primal_face_id = primal_cell_to_face[i_face];
- if (primal_face_id == face_id) {
- return primal_face_cell_is_reversed(cell_id, i_face);
- }
- }
- Assert(false, "cannot find cell's face");
- return false;
- }();
-
- const auto& primal_face_to_node = primal_face_to_node_matrix[face_id];
-
- for (size_t i_face_node = 0; i_face_node < primal_face_to_node.size(); ++i_face_node) {
- if (primal_face_to_node[i_face_node] == node_id) {
- const TinyVector<Dimension> nil = [=] {
- if (is_face_reversed_for_cell_i) {
- return -primal_nlr(face_id, i_face_node);
- } else {
- return primal_nlr(face_id, i_face_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] * value_list[i_node] * (nil, Clr);
- }
- }
- }
- }
- }
- }
- }
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ b[face_dof_number[face_id]] += value_list[i_face]; // sign
}
}
},
boundary_condition);
}
}
-
+ // EL b^L
for (const auto& boundary_condition : boundary_condition_list) {
std::visit(
[&](auto&& bc) {
@@ -685,93 +564,77 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
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];
- }
- }
+ b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face]; // sign
}
}
},
boundary_condition);
}
- 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 = 0;
-
- BiCGStab{b, A, T, 1000, 1e-15};
-
- Vector r = A * T - b;
-
- std::cout << "real residu = " << std::sqrt((r, r)) << '\n';
+ {
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { T[cell_dof_number[cell_id]] = Tj[cell_id]; });
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ if (primal_face_is_on_boundary[face_id]) {
+ T[face_dof_number[face_id]] = Tl[face_id];
+ }
+ });
+ T = 1;
+ }
+ BiCGStab{b, A, T, 1000, 1e-9};
CellValue<double> Temperature{mesh->connectivity()};
-
+ FaceValue<double> Temperature_face{mesh->connectivity()};
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_dof_number[cell_id]]; });
-
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ if (primal_face_is_neumann[face_id]) {
+ Temperature_face[face_id] = T[face_dof_number[face_id]];
+ } else {
+ Temperature_face[face_id] = Tl[face_id];
+ }
+ });
Vector<double> error{mesh->numberOfCells()};
+ CellValue<double> cell_error{mesh->connectivity()};
+ Vector<double> face_error{mesh->numberOfFaces()};
+ double error_max = 0.;
+ size_t cell_max = 0;
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
- error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]);
+ error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]);
+ cell_error[cell_id] = (Temperature[cell_id] - Tj[cell_id]);
+ });
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ if (primal_face_is_on_boundary[face_id]) {
+ face_error[face_id] = (Temperature_face[face_id] - Tl[face_id]) * sqrt(mes_l[face_id]);
+ } else {
+ face_error[face_id] = 0;
+ }
});
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); cell_id++) {
+ if (error_max < std::abs(cell_error[cell_id])) {
+ error_max = std::abs(cell_error[cell_id]);
+ cell_max = cell_id;
+ }
+ }
- std::cout << "||Error||_2 = " << std::sqrt((error, error)) << "\n";
+ std::cout << " ||Error||_max (cell)= " << error_max << " on cell " << cell_max << "\n";
+ std::cout << "||Error||_2 (cell)= " << std::sqrt((error, error)) << "\n";
+ std::cout << "||Error||_2 (face)= " << std::sqrt((face_error, face_error)) << "\n";
+ std::cout << "||Error||_2 (total)= " << std::sqrt((error, error)) + std::sqrt((face_error, face_error)) << "\n";
{
VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01);
vtk_writer.write(mesh,
{NamedItemValue{"T", Temperature}, NamedItemValue{"Vj", primal_Vj},
- NamedItemValue{"Texact", Tj}},
+ NamedItemValue{"Texact", Tj}, NamedItemValue{"error", cell_error}},
0,
true); // forces last output
}
@@ -786,13 +649,13 @@ class HeatDiamondScheme<Dimension>::DirichletBoundaryCondition
{
private:
const Array<const double> m_value_list;
- const Array<const NodeId> m_node_list;
+ const Array<const FaceId> m_face_list;
public:
- const Array<const NodeId>&
- nodeList() const
+ const Array<const FaceId>&
+ faceList() const
{
- return m_node_list;
+ return m_face_list;
}
const Array<const double>&
@@ -801,10 +664,10 @@ class HeatDiamondScheme<Dimension>::DirichletBoundaryCondition
return m_value_list;
}
- DirichletBoundaryCondition(const Array<const NodeId>& node_list, const Array<const double>& value_list)
- : m_value_list{value_list}, m_node_list{node_list}
+ DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
{
- Assert(m_value_list.size() == m_node_list.size());
+ Assert(m_value_list.size() == m_face_list.size());
}
~DirichletBoundaryCondition() = default;
--
GitLab