From c612d4c04e2981067bc2959568c476bd0de0bf0d Mon Sep 17 00:00:00 2001
From: Julie PATELA <julie.patela.ocre@cea.fr>
Date: Wed, 5 Aug 2020 16:52:40 +0200
Subject: [PATCH] Add interpolation of boundary nodes for Neumann BC.
---
.../algorithms/HeatDiamondAlgorithm.cpp | 311 ++++++++++++++----
1 file changed, 255 insertions(+), 56 deletions(-)
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index aef3cd1ae..be6552c27 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -52,6 +52,11 @@ 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;
@@ -81,6 +86,14 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
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});
}
}
@@ -93,9 +106,30 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
break;
}
case IBoundaryConditionDescriptor::Type::neumann: {
- // const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor =
- // dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_descriptor);
- throw NotImplementedError("NIY");
+ const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor =
+ dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*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 == neumann_bc_descriptor.boundaryDescriptor()) {
+ const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
+
+ if constexpr (Dimension > 1) {
+ 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::face>(g_id,
+ 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");
+ }
+ }
+ }
break;
}
default: {
@@ -116,8 +150,66 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
return compute_cell_dof_number;
}();
+ size_t number_of_dof = mesh->numberOfCells();
+
+ const FaceValue<const size_t> face_dof_number = [&] {
+ FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
+ compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
+ 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];
+ if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
+ std::ostringstream os;
+ os << "The face " << face_id << " is used at least twice for boundary conditions";
+ throw NormalError(os.str());
+ } else {
+ compute_face_dof_number[face_id] = number_of_dof++;
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return compute_face_dof_number;
+ }();
+
+ 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());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
+ }
+
+ primal_node_is_on_boundary.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (face_to_cell_matrix[face_id].size() == 1) {
+ 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];
+ primal_node_is_on_boundary[node_id] = true;
+ }
+ }
+ }
- auto dof_number = [&](CellId cell_id) { return cell_dof_number[cell_id]; };
+ FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
+ }
+
+ primal_face_is_on_boundary.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (face_to_cell_matrix[face_id].size() == 1) {
+ primal_face_is_on_boundary[face_id] = true;
+ }
+ }
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
@@ -126,9 +218,17 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
NodeValue<double> Tr(mesh->connectivity());
const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
+
+ const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
CellValuePerNode<double> w_rj{mesh->connectivity()};
+ FaceValuePerNode<double> w_rl{mesh->connectivity()};
+
+ for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
+ w_rl[i] = std::numeric_limits<double>::signaling_NaN();
+ }
for (NodeId i_node = 0; i_node < mesh->numberOfNodes(); ++i_node) {
Vector<double> b{Dimension + 1};
@@ -138,27 +238,77 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
const auto& node_to_cell = node_to_cell_matrix[i_node];
- LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size()};
- for (size_t j = 0; j < node_to_cell.size(); j++) {
- A(0, j) = 1;
- }
- for (size_t i = 1; i < Dimension + 1; i++) {
+ if (not primal_node_is_on_boundary[i_node]) {
+ LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size()};
for (size_t j = 0; j < node_to_cell.size(); j++) {
- const CellId J = node_to_cell[j];
- A(i, j) = xj[J][i - 1];
+ A(0, j) = 1;
}
- }
- Vector<double> x{node_to_cell.size()};
- x = 0;
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ const CellId J = node_to_cell[j];
+ A(i, j) = xj[J][i - 1];
+ }
+ }
+ Vector<double> x{node_to_cell.size()};
+ x = 0;
- LocalRectangularMatrix B = transpose(A) * A;
- Vector y = transpose(A) * b;
- PCG<false>{y, B, B, x, 10, 1e-12};
+ LocalRectangularMatrix B = transpose(A) * A;
+ Vector y = transpose(A) * b;
+ PCG<false>{y, B, B, x, 10, 1e-12};
- Tr[i_node] = 0;
- for (size_t j = 0; j < node_to_cell.size(); j++) {
- Tr[i_node] += x[j] * Tj[node_to_cell[j]];
- w_rj(i_node, j) = x[j];
+ Tr[i_node] = 0;
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ Tr[i_node] += x[j] * Tj[node_to_cell[j]];
+ w_rj(i_node, j) = x[j];
+ }
+ } else {
+ int nb_face_used = 0;
+ 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]) {
+ 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++) {
+ A(0, j) = 1;
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ const CellId J = node_to_cell[j];
+ A(i, j) = xj[J][i - 1];
+ }
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ int cpt_face = 0;
+ 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]) {
+ A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
+ cpt_face++;
+ }
+ }
+ }
+
+ Vector<double> x{node_to_cell.size() + nb_face_used};
+ x = 0;
+
+ LocalRectangularMatrix B = transpose(A) * A;
+ Vector y = transpose(A) * b;
+ PCG<false>{y, B, B, x, 10, 1e-12};
+
+ Tr[i_node] = 0;
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ Tr[i_node] += x[j] * Tj[node_to_cell[j]];
+ w_rj(i_node, j) = x[j];
+ }
+ int cpt_face = 0;
+ 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++];
+ }
+ }
}
}
@@ -243,26 +393,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
return computed_alpha_l;
}();
-
- 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());
- if (parallel::size() > 1) {
- throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
- }
-
- primal_node_is_on_boundary.fill(false);
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- if (face_to_cell_matrix[face_id].size() == 1) {
- 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];
- primal_node_is_on_boundary[node_id] = true;
- }
- }
- }
-
- SparseMatrixDescriptor S{mesh->numberOfCells()};
+ 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];
@@ -272,9 +403,9 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
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(dof_number(cell_id1), dof_number(cell_id2)) += beta_l;
+ S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
} else {
- S(dof_number(cell_id1), dof_number(cell_id2)) -= beta_l;
+ S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
}
}
}
@@ -328,7 +459,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
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 primal_node_is_on_boundary[node_id]) {
+ if (not is_dirichlet[node_id]) {
const TinyVector<Dimension> nil = [&] {
if (is_face_reversed_for_cell_i) {
return -primal_nlr(face_id, i_node);
@@ -349,7 +480,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
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(dof_number(i_id), 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;
}
}
}
@@ -365,11 +496,11 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
const CellValue<const double> primal_Vj = mesh_data.Vj();
CRSMatrix A{S};
- Vector<double> b{mesh->numberOfCells()};
+ 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[dof_number(cell_id)] = fj[cell_id] * primal_Vj[cell_id];
+ b[cell_dof_number[cell_id]] = fj[cell_id] * primal_Vj[cell_id];
}
{ // Dirichlet
@@ -426,7 +557,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
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[dof_number(cell_id)] += alpha_l[face_id] * value_list[i_node] * (nil, Clr);
+ b[cell_dof_number[cell_id]] += alpha_l[face_id] * value_list[i_node] * (nil, Clr);
}
}
}
@@ -439,10 +570,9 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
},
boundary_condition);
}
- // std::exit(0);
}
- Vector<double> T{mesh->numberOfCells()};
+ Vector<double> T{number_of_dof};
T = 0;
BiCGStab{b, A, T, 1000, 1e-15};
@@ -450,7 +580,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
CellValue<double> Temperature{mesh->connectivity()};
parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[dof_number(cell_id)]; });
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_dof_number[cell_id]]; });
Vector<double> error{mesh->numberOfCells()};
parallel_for(
@@ -497,29 +627,98 @@ class HeatDiamondScheme<Dimension>::DirichletBoundaryCondition
DirichletBoundaryCondition(const Array<const NodeId>& node_list, const Array<const double>& value_list)
: m_value_list{value_list}, m_node_list{node_list}
- {}
+ {
+ Assert(m_value_list.size() == m_node_list.size());
+ }
~DirichletBoundaryCondition() = default;
};
template <size_t Dimension>
-class HeatDiamondScheme<Dimension>::FourierBoundaryCondition
+class HeatDiamondScheme<Dimension>::NeumannBoundaryCondition
{
+ private:
+ const Array<const double> m_value_list;
+ const Array<const FaceId> m_face_list;
+
public:
- ~FourierBoundaryCondition() = default;
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ NeumannBoundaryCondition(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_face_list.size());
+ }
+
+ ~NeumannBoundaryCondition() = default;
};
template <size_t Dimension>
-class HeatDiamondScheme<Dimension>::NeumannBoundaryCondition
+class HeatDiamondScheme<Dimension>::FourierBoundaryCondition
{
+ private:
+ const Array<const double> m_coef_list;
+ const Array<const double> m_value_list;
+ const Array<const FaceId> m_face_list;
+
public:
- ~NeumannBoundaryCondition() = default;
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ const Array<const double>&
+ coefList() const
+ {
+ return m_coef_list;
+ }
+
+ public:
+ FourierBoundaryCondition(const Array<const FaceId>& face_list,
+ const Array<const double>& coef_list,
+ const Array<const double>& value_list)
+ : m_coef_list{coef_list}, m_value_list{value_list}, m_face_list{face_list}
+ {
+ Assert(m_coef_list.size() == m_face_list.size());
+ Assert(m_value_list.size() == m_face_list.size());
+ }
+
+ ~FourierBoundaryCondition() = default;
};
template <size_t Dimension>
class HeatDiamondScheme<Dimension>::SymmetryBoundaryCondition
{
+ private:
+ const Array<const FaceId> m_face_list;
+
public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ public:
+ SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
+
~SymmetryBoundaryCondition() = default;
};
--
GitLab