Skip to content
Snippets Groups Projects
Commit df9d6468 authored by Axelle Drouard's avatar Axelle Drouard
Browse files

Interpolate BC at the nodes instead of faces

parent 5f95d5bd
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/scheme/ScalarNodalScheme.cpp
+ 56
− 124
View file @ df9d6468
......@@ -53,7 +53,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const Array<const double>& value_list)
: m_value_list{value_list}, m_face_list{face_list}, m_node_list{node_list}
{
Assert(m_value_list.size() == m_face_list.size());
Assert(m_value_list.size() == m_node_list.size());
}
~DirichletBoundaryCondition() = default;
......@@ -90,7 +90,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const Array<const double>& value_list)
: m_value_list{value_list}, m_face_list{face_list}, m_node_list{node_list}
{
Assert(m_value_list.size() == m_face_list.size());
Assert(m_value_list.size() == m_node_list.size());
}
~NeumannBoundaryCondition() = default;
......@@ -101,13 +101,13 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
private:
const Array<const double> m_coef_list;
const Array<const double> m_value_list;
const Array<const FaceId> m_face_list;
const Array<const NodeId> m_node_list;
public:
const Array<const FaceId>&
faceList() const
const Array<const NodeId>&
nodeList() const
{
return m_face_list;
return m_node_list;
}
const Array<const double>&
......@@ -123,13 +123,13 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
public:
FourierBoundaryCondition(const Array<const FaceId>& face_list,
FourierBoundaryCondition(const Array<const NodeId>& node_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}
: m_coef_list{coef_list}, m_value_list{value_list}, m_node_list{node_list}
{
Assert(m_coef_list.size() == m_face_list.size());
Assert(m_value_list.size() == m_face_list.size());
Assert(m_coef_list.size() == m_node_list.size());
Assert(m_value_list.size() == m_node_list.size());
}
~FourierBoundaryCondition() = default;
......@@ -173,13 +173,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
getMeshNodeBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
// MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
Array<const double> value_list =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
mesh_data.xl(),
mesh_face_boundary
.faceList());
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
mesh_node_boundary
.nodeList());
boundary_condition_list.push_back(
DirichletBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
......@@ -203,13 +202,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
getMeshNodeBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
// MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
Array<const double> value_list =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
mesh_data.xl(),
mesh_face_boundary
.faceList());
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
mesh_node_boundary
.nodeList());
boundary_condition_list.push_back(
NeumannBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
......@@ -246,13 +244,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
getMeshNodeBoundary(*mesh, dirichlet_bc_prev_descriptor.boundaryDescriptor());
const FunctionSymbolId g_id = dirichlet_bc_prev_descriptor.rhsSymbolId();
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
// MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
Array<const double> value_list =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
mesh_data.xl(),
mesh_face_boundary
.faceList());
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
mesh_node_boundary
.nodeList());
boundary_prev_condition_list.push_back(
DirichletBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
......@@ -276,13 +273,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
getMeshNodeBoundary(*mesh, neumann_bc_prev_descriptor.boundaryDescriptor());
const FunctionSymbolId g_id = neumann_bc_prev_descriptor.rhsSymbolId();
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
// MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
Array<const double> value_list =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
mesh_data.xl(),
mesh_face_boundary
.faceList());
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(),
mesh_node_boundary
.nodeList());
boundary_prev_condition_list.push_back(
NeumannBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
......@@ -316,7 +312,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const auto is_boundary_face = mesh->connectivity().isBoundaryFace();
const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
// const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
const auto& node_local_numbers_in_their_cells = mesh->connectivity().nodeLocalNumbersInTheirCells();
const CellValue<const double> Vj = mesh_data.Vj();
......@@ -331,18 +327,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
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];
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
const NodeId node_id = face_nodes[i_node];
const auto& node_list = bc.nodeList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
compute_node_is_neumann[node_id] = true;
}
}
}
},
boundary_condition);
}
......@@ -357,19 +347,14 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
const auto& face_list = bc.faceList();
const auto& node_list = bc.nodeList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
const NodeId node_id = face_nodes[i_node];
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
compute_node_is_dirichlet[node_id] = true;
}
}
}
},
boundary_condition);
}
......@@ -486,57 +471,28 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
const auto& face_list = bc.faceList();
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
const NodeId node_id = face_nodes[i_node];
if (not node_is_dirichlet[node_id]) {
if (node_is_angle[node_id]) {
compute_node_boundary_values[node_id] +=
value_list[i_face] * std::abs(dot(nl[face_id], exterior_normal[node_id]));
} else {
compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
sum_mes_l[node_id] += mes_l[face_id];
}
}
}
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
compute_node_boundary_values[node_id] = value_list[i_node];
}
} else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
const auto& face_list = bc.faceList();
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
const NodeId node_id = face_nodes[i_node];
if (not node_is_neumann[node_id]) {
compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
// compute_node_boundary_values[node_id] = xr[node_id][0];
sum_mes_l[node_id] += mes_l[face_id];
} else {
compute_node_boundary_values[node_id] = value_list[i_face];
}
}
compute_node_boundary_values[node_id] = value_list[i_node];
}
}
},
boundary_condition);
}
for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id]) and not node_is_angle[node_id]) {
compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
} else if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
}
// std::cout << node_id << " " << compute_node_boundary_values[node_id] << "\n";
}
return compute_node_boundary_values;
}();
......@@ -550,51 +506,27 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
const auto& face_list = bc.faceList();
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
const NodeId node_id = face_nodes[i_node];
if (not node_is_dirichlet[node_id]) {
compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
sum_mes_l[node_id] += mes_l[face_id];
compute_node_boundary_values[node_id] = value_list[i_node];
}
}
}
} else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
const auto& face_list = bc.faceList();
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
const NodeId node_id = face_nodes[i_node];
if (not node_is_neumann[node_id]) {
compute_node_boundary_values[node_id] += value_list[i_face] * mes_l[face_id];
sum_mes_l[node_id] += mes_l[face_id];
} else {
compute_node_boundary_values[node_id] = value_list[i_face];
}
}
compute_node_boundary_values[node_id] = value_list[i_node];
}
}
},
boundary_condition);
}
for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
if ((not node_is_dirichlet[node_id]) && (node_is_neumann[node_id])) {
compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
} else if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
}
}
return compute_node_boundary_values;
}();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment