Skip to content
Snippets Groups Projects
Commit 96e80be0 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Finalize boundary condition rewriting for acoustic solver (Glace)

Now boundary condition treatment is embedded within the scheme which
makes lots more sense
parent cfe40c29
No related branches found
No related tags found
1 merge request!43Feature/glace boundary conditions
...@@ -112,7 +112,7 @@ struct GlaceScheme ...@@ -112,7 +112,7 @@ struct GlaceScheme
{ {
std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n'; std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
std::vector<BoundaryConditionHandler> bc_list; typename AcousticSolver<MeshType>::BoundaryConditionList bc_list;
{ {
constexpr ItemType FaceType = [] { constexpr ItemType FaceType = [] {
if constexpr (Dimension > 1) { if constexpr (Dimension > 1) {
...@@ -125,6 +125,8 @@ struct GlaceScheme ...@@ -125,6 +125,8 @@ struct GlaceScheme
for (const auto& bc_descriptor : bc_descriptor_list) { for (const auto& bc_descriptor : bc_descriptor_list) {
switch (bc_descriptor->type()) { switch (bc_descriptor->type()) {
case IBoundaryConditionDescriptor::Type::symmetry: { case IBoundaryConditionDescriptor::Type::symmetry: {
using SymmetryBoundaryCondition = typename AcousticSolver<MeshType>::SymmetryBoundaryCondition;
const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor = const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor); dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
for (size_t i_ref_face_list = 0; for (size_t i_ref_face_list = 0;
...@@ -132,16 +134,15 @@ struct GlaceScheme ...@@ -132,16 +134,15 @@ struct GlaceScheme
const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list); const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
const RefId& ref = ref_face_list.refId(); const RefId& ref = ref_face_list.refId();
if (ref == sym_bc_descriptor.boundaryDescriptor()) { if (ref == sym_bc_descriptor.boundaryDescriptor()) {
SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc = bc_list.push_back(
new SymmetryBoundaryCondition<MeshType::Dimension>( SymmetryBoundaryCondition{MeshFlatNodeBoundary<MeshType::Dimension>(m_mesh, ref_face_list)});
MeshFlatNodeBoundary<MeshType::Dimension>(m_mesh, ref_face_list));
std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
bc_list.push_back(BoundaryConditionHandler(bc));
} }
} }
break; break;
} }
case IBoundaryConditionDescriptor::Type::velocity: { case IBoundaryConditionDescriptor::Type::velocity: {
using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition;
const VelocityBoundaryConditionDescriptor& velocity_bc_descriptor = const VelocityBoundaryConditionDescriptor& velocity_bc_descriptor =
dynamic_cast<const VelocityBoundaryConditionDescriptor&>(*bc_descriptor); dynamic_cast<const VelocityBoundaryConditionDescriptor&>(*bc_descriptor);
for (size_t i_ref_face_list = 0; for (size_t i_ref_face_list = 0;
...@@ -149,43 +150,51 @@ struct GlaceScheme ...@@ -149,43 +150,51 @@ struct GlaceScheme
const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list); const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
const RefId& ref = ref_face_list.refId(); const RefId& ref = ref_face_list.refId();
if (ref == velocity_bc_descriptor.boundaryDescriptor()) { if (ref == velocity_bc_descriptor.boundaryDescriptor()) {
MeshNodeBoundary<Dimension> mesh_node_boundary{m_mesh, ref_face_list};
const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId(); const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId();
if constexpr (Dimension == 1) { const auto& node_list = mesh_node_boundary.nodeList();
const auto& node_list = ref_face_list.list();
Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>( Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
TinyVector<Dimension>)>::template interpolate<FaceType>(velocity_id, m_mesh->xr(), node_list); TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, m_mesh->xr(), node_list);
std::shared_ptr bc = bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
std::make_shared<VelocityBoundaryCondition<MeshType::Dimension>>(node_list, value_list);
bc_list.push_back(BoundaryConditionHandler(bc));
} else {
const auto& face_list = ref_face_list.list();
const auto& face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
std::set<NodeId> node_set;
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
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) {
node_set.insert(face_nodes[i_node]);
}
}
Array<NodeId> node_list = convert_to_array(node_set); // if constexpr (Dimension == 1) {
// const auto& node_list = ref_face_list.list();
Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>( // Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, m_mesh->xr(), node_list); // TinyVector<Dimension>)>::template interpolate<FaceType>(velocity_id, m_mesh->xr(), node_list);
std::shared_ptr bc = // bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
std::make_shared<VelocityBoundaryCondition<MeshType::Dimension>>(node_list, value_list); // } else {
bc_list.push_back(BoundaryConditionHandler(bc)); // const auto& face_list = ref_face_list.list();
} // const auto& face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
// std::set<NodeId> node_set;
// for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
// 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) {
// node_set.insert(face_nodes[i_node]);
// }
// }
// Array<NodeId> node_list = convert_to_array(node_set);
// Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
// TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, m_mesh->xr(),
// node_list);
// bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
// }
} }
} }
break; break;
} }
case IBoundaryConditionDescriptor::Type::pressure: { case IBoundaryConditionDescriptor::Type::pressure: {
using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition;
const PressureBoundaryConditionDescriptor& pressure_bc_descriptor = const PressureBoundaryConditionDescriptor& pressure_bc_descriptor =
dynamic_cast<const PressureBoundaryConditionDescriptor&>(*bc_descriptor); dynamic_cast<const PressureBoundaryConditionDescriptor&>(*bc_descriptor);
for (size_t i_ref_face_list = 0; for (size_t i_ref_face_list = 0;
...@@ -209,9 +218,7 @@ struct GlaceScheme ...@@ -209,9 +218,7 @@ struct GlaceScheme
} }
}(); }();
std::shared_ptr bc = bc_list.push_back(PressureBoundaryCondition{face_list, face_values});
std::make_shared<PressureBoundaryCondition<MeshType::Dimension>>(face_list, face_values);
bc_list.push_back(BoundaryConditionHandler(bc));
} }
} }
break; break;
...@@ -283,12 +290,12 @@ struct GlaceScheme ...@@ -283,12 +290,12 @@ struct GlaceScheme
while ((t < tmax) and (iteration < itermax)) { while ((t < tmax) and (iteration < itermax)) {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh); MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
vtk_writer.write(m_mesh, // vtk_writer.write(m_mesh,
{NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()}, // {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()}, // NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()}, // NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}}, // NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
t); // t);
double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj); double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
if (t + dt > tmax) { if (t + dt > tmax) {
dt = tmax - t; dt = tmax - t;
......
#ifndef ACOUSTIC_SOLVER_HPP #ifndef ACOUSTIC_SOLVER_HPP
#define ACOUSTIC_SOLVER_HPP #define ACOUSTIC_SOLVER_HPP
#include <rang.hpp>
#include <utils/ArrayUtils.hpp>
#include <scheme/BlockPerfectGas.hpp>
#include <utils/PugsAssert.hpp>
#include <scheme/BoundaryCondition.hpp>
#include <scheme/FiniteVolumesEulerUnknowns.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <algebra/TinyMatrix.hpp> #include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp> #include <algebra/TinyVector.hpp>
#include <mesh/ItemValueUtils.hpp> #include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/SubItemValuePerItem.hpp> #include <mesh/SubItemValuePerItem.hpp>
#include <scheme/BlockPerfectGas.hpp>
#include <scheme/FiniteVolumesEulerUnknowns.hpp>
#include <utils/ArrayUtils.hpp>
#include <utils/Exceptions.hpp> #include <utils/Exceptions.hpp>
#include <utils/Messenger.hpp> #include <utils/Messenger.hpp>
#include <utils/PugsAssert.hpp>
#include <rang.hpp>
#include <iostream> #include <iostream>
template <typename MeshType> template <typename MeshType>
class AcousticSolver class AcousticSolver
{ {
public:
class PressureBoundaryCondition;
class SymmetryBoundaryCondition;
class VelocityBoundaryCondition;
using BoundaryCondition =
std::variant<PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
constexpr static size_t Dimension = MeshType::Dimension; constexpr static size_t Dimension = MeshType::Dimension;
private:
using MeshDataType = MeshData<Dimension>; using MeshDataType = MeshData<Dimension>;
using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>; using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;
std::shared_ptr<const MeshType> m_mesh; std::shared_ptr<const MeshType> m_mesh;
const typename MeshType::Connectivity& m_connectivity; const typename MeshType::Connectivity& m_connectivity;
const std::vector<BoundaryConditionHandler>& m_boundary_condition_list;
using Rd = TinyVector<Dimension>; using Rd = TinyVector<Dimension>;
using Rdd = TinyMatrix<Dimension>; using Rdd = TinyMatrix<Dimension>;
private: const BoundaryConditionList m_boundary_condition_list;
void
_applyPressureBC()
{
for (const auto& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
if constexpr (Dimension == 1) {
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const auto& node_to_cell_matrix = m_connectivity.nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
const auto& node_list = bc.faceList();
const auto& value_list = bc.valueList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_cell_list = node_to_cell_matrix[node_id];
Assert(node_cell_list.size() == 1);
CellId node_cell_id = node_cell_list[0];
size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
m_br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
});
} else {
const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
const auto& face_to_cell_matrix = m_connectivity.faceToCellMatrix();
const auto& face_to_node_matrix = m_connectivity.faceToNodeMatrix();
const auto& face_local_numbers_in_their_cells = m_connectivity.faceLocalNumbersInTheirCells();
const auto& face_cell_is_reversed = m_connectivity.cellFaceIsReversed();
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) {
const FaceId face_id = face_list[i_face];
const auto& face_cell_list = face_to_cell_matrix[face_id];
Assert(face_cell_list.size() == 1);
CellId face_cell_id = face_cell_list[0];
size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
NodeId node_id = face_nodes[i_node];
m_br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
}
}
}
}
},
boundary_condition);
}
}
void
_applySymmetryBC()
{
for (const auto& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
const Rd& n = bc.outgoingNormal();
const Rdd I = identity;
const Rdd nxn = tensorProduct(n, n);
const Rdd P = I - nxn;
const Array<const NodeId>& node_list = bc.nodeList();
parallel_for(
bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
const NodeId r = node_list[r_number];
m_Ar[r] = P * m_Ar[r] * P + nxn;
m_br[r] = P * m_br[r];
});
}
},
boundary_condition);
}
}
void
_applyVelocityBC()
{
for (const auto& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
NodeId node_id = node_list[i_node];
const auto& value = value_list[i_node];
m_Ar[node_id] = identity;
m_br[node_id] = value;
});
}
},
boundary_condition);
}
}
PUGS_INLINE PUGS_INLINE
const CellValue<const double> const CellValue<const double>
computeRhoCj(const CellValue<const double>& rhoj, const CellValue<const double>& cj) computeRhoCj(const CellValue<const double>& rhoj, const CellValue<const double>& cj)
...@@ -121,100 +241,9 @@ class AcousticSolver ...@@ -121,100 +241,9 @@ class AcousticSolver
void void
applyBoundaryConditions() applyBoundaryConditions()
{ {
for (const auto& handler : m_boundary_condition_list) { this->_applyPressureBC();
switch (handler.boundaryCondition().type()) { this->_applySymmetryBC();
case BoundaryCondition::velocity: { this->_applyVelocityBC();
const VelocityBoundaryCondition<Dimension>& velocity_bc =
dynamic_cast<const VelocityBoundaryCondition<Dimension>&>(handler.boundaryCondition());
const auto& node_list = velocity_bc.nodeList();
const auto& value_list = velocity_bc.valueList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
NodeId node_id = node_list[i_node];
const auto& value = value_list[i_node];
m_Ar[node_id] = identity;
m_br[node_id] = value;
});
break;
}
case BoundaryCondition::pressure: {
const PressureBoundaryCondition<Dimension>& pressure_bc =
dynamic_cast<const PressureBoundaryCondition<Dimension>&>(handler.boundaryCondition());
MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
if constexpr (Dimension == 1) {
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const auto& node_to_cell_matrix = m_connectivity.nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
const auto& node_list = pressure_bc.faceList();
const auto& value_list = pressure_bc.valueList();
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_cell_list = node_to_cell_matrix[node_id];
Assert(node_cell_list.size() == 1);
CellId node_cell_id = node_cell_list[0];
size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
m_br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
});
} else {
const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
const auto& face_to_cell_matrix = m_connectivity.faceToCellMatrix();
const auto& face_to_node_matrix = m_connectivity.faceToNodeMatrix();
const auto& face_local_numbers_in_their_cells = m_connectivity.faceLocalNumbersInTheirCells();
const auto& face_cell_is_reversed = m_connectivity.cellFaceIsReversed();
const auto& face_list = pressure_bc.faceList();
const auto& value_list = pressure_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_cell_list = face_to_cell_matrix[face_id];
Assert(face_cell_list.size() == 1);
CellId face_cell_id = face_cell_list[0];
size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
const auto& face_nodes = face_to_node_matrix[face_id];
for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
NodeId node_id = face_nodes[i_node];
m_br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
}
}
}
break;
}
case BoundaryCondition::symmetry: {
const SymmetryBoundaryCondition<Dimension>& symmetry_bc =
dynamic_cast<const SymmetryBoundaryCondition<Dimension>&>(handler.boundaryCondition());
const Rd& n = symmetry_bc.outgoingNormal();
const Rdd I = identity;
const Rdd nxn = tensorProduct(n, n);
const Rdd P = I - nxn;
const Array<const NodeId>& node_list = symmetry_bc.nodeList();
parallel_for(
symmetry_bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
const NodeId r = node_list[r_number];
m_Ar[r] = P * m_Ar[r] * P + nxn;
m_br[r] = P * m_br[r];
});
break;
}
}
}
} }
NodeValue<Rd> NodeValue<Rd>
...@@ -291,7 +320,7 @@ class AcousticSolver ...@@ -291,7 +320,7 @@ class AcousticSolver
CellValue<double> m_Vj_over_cj; CellValue<double> m_Vj_over_cj;
public: public:
AcousticSolver(std::shared_ptr<const MeshType> p_mesh, const std::vector<BoundaryConditionHandler>& bc_list) AcousticSolver(std::shared_ptr<const MeshType> p_mesh, const BoundaryConditionList& bc_list)
: m_mesh(p_mesh), : m_mesh(p_mesh),
m_connectivity(m_mesh->connectivity()), m_connectivity(m_mesh->connectivity()),
m_boundary_condition_list(bc_list), m_boundary_condition_list(bc_list),
...@@ -388,4 +417,124 @@ class AcousticSolver ...@@ -388,4 +417,124 @@ class AcousticSolver
} }
}; };
template <typename MeshType>
class AcousticSolver<MeshType>::PressureBoundaryCondition
{
private:
const Array<const double> m_value_list;
const Array<const FaceId> m_face_list;
public:
const Array<const FaceId>&
faceList() const
{
return m_face_list;
}
const Array<const double>&
valueList() const
{
return m_value_list;
}
PressureBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
: m_value_list{value_list}, m_face_list{face_list}
{}
~PressureBoundaryCondition() = default;
};
template <>
class AcousticSolver<Mesh<Connectivity<1>>>::PressureBoundaryCondition
{
private:
const Array<const double> m_value_list;
const Array<const NodeId> m_face_list;
public:
const Array<const NodeId>&
faceList() const
{
return m_face_list;
}
const Array<const double>&
valueList() const
{
return m_value_list;
}
PressureBoundaryCondition(const Array<const NodeId>& face_list, const Array<const double>& value_list)
: m_value_list{value_list}, m_face_list{face_list}
{}
~PressureBoundaryCondition() = default;
};
template <typename MeshType>
class AcousticSolver<MeshType>::VelocityBoundaryCondition
{
private:
const Array<const TinyVector<Dimension>> m_value_list;
const Array<const NodeId> m_node_list;
public:
const Array<const NodeId>&
nodeList() const
{
return m_node_list;
}
const Array<const TinyVector<Dimension>>&
valueList() const
{
return m_value_list;
}
VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
: m_value_list{value_list}, m_node_list{node_list}
{}
~VelocityBoundaryCondition() = default;
};
template <typename MeshType>
class AcousticSolver<MeshType>::SymmetryBoundaryCondition
{
public:
static constexpr size_t Dimension = MeshType::Dimension;
using Rd = TinyVector<Dimension, double>;
private:
const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
public:
const Rd&
outgoingNormal() const
{
return m_mesh_flat_node_boundary.outgoingNormal();
}
size_t
numberOfNodes() const
{
return m_mesh_flat_node_boundary.nodeList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
: m_mesh_flat_node_boundary(mesh_flat_node_boundary)
{
;
}
~SymmetryBoundaryCondition() = default;
};
#endif // ACOUSTIC_SOLVER_HPP #endif // ACOUSTIC_SOLVER_HPP
#ifndef BOUNDARY_CONDITION_HANDLER_HPP
#define BOUNDARY_CONDITION_HANDLER_HPP
#include <utils/Array.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/RefItemList.hpp>
#include <memory>
#include <vector>
class BoundaryCondition // clazy:exclude=copyable-polymorphic
{
public:
enum Type
{
velocity,
pressure,
symmetry
};
const Type m_type;
protected:
BoundaryCondition(Type type) : m_type(type)
{
;
}
public:
const Type&
type() const
{
return m_type;
}
virtual ~BoundaryCondition() = default;
};
template <size_t Dimension>
class PressureBoundaryCondition final : public BoundaryCondition // clazy:exclude=copyable-polymorphic
{
private:
const Array<const double> m_value_list;
const Array<const FaceId> m_face_list;
public:
const Array<const FaceId>&
faceList() const
{
return m_face_list;
}
const Array<const double>&
valueList() const
{
return m_value_list;
}
PressureBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
: BoundaryCondition{BoundaryCondition::pressure}, m_value_list{value_list}, m_face_list{face_list}
{}
~PressureBoundaryCondition() = default;
};
template <>
class PressureBoundaryCondition<1> final : public BoundaryCondition // clazy:exclude=copyable-polymorphic
{
private:
const Array<const double> m_value_list;
const Array<const NodeId> m_face_list;
public:
const Array<const NodeId>&
faceList() const
{
return m_face_list;
}
const Array<const double>&
valueList() const
{
return m_value_list;
}
PressureBoundaryCondition(const Array<const NodeId>& face_list, const Array<const double>& value_list)
: BoundaryCondition{BoundaryCondition::pressure}, m_value_list{value_list}, m_face_list{face_list}
{}
~PressureBoundaryCondition() = default;
};
template <size_t dimension>
class SymmetryBoundaryCondition : public BoundaryCondition
{
public:
using Rd = TinyVector<dimension, double>;
private:
const MeshFlatNodeBoundary<dimension> m_mesh_flat_node_boundary;
public:
const Rd&
outgoingNormal() const
{
return m_mesh_flat_node_boundary.outgoingNormal();
}
size_t
numberOfNodes() const
{
return m_mesh_flat_node_boundary.nodeList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
SymmetryBoundaryCondition(const MeshFlatNodeBoundary<dimension>& mesh_flat_node_boundary)
: BoundaryCondition(BoundaryCondition::symmetry), m_mesh_flat_node_boundary(mesh_flat_node_boundary)
{
;
}
~SymmetryBoundaryCondition() = default;
};
template <size_t Dimension>
class VelocityBoundaryCondition final : public BoundaryCondition // clazy:exclude=copyable-polymorphic
{
private:
const Array<const TinyVector<Dimension>> m_value_list;
const Array<const NodeId> m_node_list;
public:
const Array<const NodeId>&
nodeList() const
{
return m_node_list;
}
const Array<const TinyVector<Dimension>>&
valueList() const
{
return m_value_list;
}
VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
: BoundaryCondition{BoundaryCondition::velocity}, m_value_list{value_list}, m_node_list{node_list}
{}
~VelocityBoundaryCondition() = default;
};
class BoundaryConditionHandler
{
private:
std::shared_ptr<BoundaryCondition> m_boundary_condition;
public:
const BoundaryCondition&
boundaryCondition() const
{
return *m_boundary_condition;
}
PUGS_INLINE
BoundaryConditionHandler& operator=(BoundaryConditionHandler&) = default;
PUGS_INLINE
BoundaryConditionHandler(const BoundaryConditionHandler&) = default;
PUGS_INLINE
BoundaryConditionHandler(BoundaryConditionHandler&&) = default;
PUGS_INLINE
BoundaryConditionHandler(std::shared_ptr<BoundaryCondition> boundary_condition)
: m_boundary_condition(boundary_condition)
{}
~BoundaryConditionHandler() = default;
};
#endif // BOUNDARY_CONDITION_HANDLER_HPP
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment