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

Add interpolation of boundary nodes for Neumann BC.

parent 2a5eb033
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
+ 255
− 56
View file @ c612d4c0
......@@ -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");
}
auto dof_number = [&](CellId cell_id) { return cell_dof_number[cell_id]; };
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;
}
}
}
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,6 +238,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
const auto& node_to_cell = node_to_cell_matrix[i_node];
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++) {
A(0, j) = 1;
......@@ -160,6 +261,55 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
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;
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment