From cfe40c29c4d3f5f3e4db81048f0dead9be29f266 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Mon, 27 Jul 2020 23:50:37 +0200
Subject: [PATCH] Add velocity boundary conditions for acoustic solver (Glace)
Still requires design improvements and to check data validity
---
src/language/modules/SchemeModule.cpp | 62 ++++++++++++++++++-
src/scheme/AcousticSolver.hpp | 36 +++++++----
src/scheme/BoundaryCondition.hpp | 32 +++++++++-
src/scheme/IBoundaryConditionDescriptor.hpp | 1 +
.../VelocityBoundaryConditionDescriptor.hpp | 55 ++++++++++++++++
5 files changed, 170 insertions(+), 16 deletions(-)
create mode 100644 src/scheme/VelocityBoundaryConditionDescriptor.hpp
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 7903e6ce4..4e19eb473 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -10,6 +10,7 @@
#include <scheme/NumberedBoundaryDescriptor.hpp>
#include <scheme/PressureBoundaryConditionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <scheme/VelocityBoundaryConditionDescriptor.hpp>
#include <memory>
@@ -140,6 +141,50 @@ struct GlaceScheme
}
break;
}
+ case IBoundaryConditionDescriptor::Type::velocity: {
+ const VelocityBoundaryConditionDescriptor& velocity_bc_descriptor =
+ dynamic_cast<const VelocityBoundaryConditionDescriptor&>(*bc_descriptor);
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<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();
+ if (ref == velocity_bc_descriptor.boundaryDescriptor()) {
+ const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId();
+
+ if constexpr (Dimension == 1) {
+ const auto& node_list = ref_face_list.list();
+
+ Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<FaceType>(velocity_id, m_mesh->xr(), node_list);
+
+ std::shared_ptr bc =
+ 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);
+
+ Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, m_mesh->xr(), node_list);
+
+ std::shared_ptr bc =
+ std::make_shared<VelocityBoundaryCondition<MeshType::Dimension>>(node_list, value_list);
+ bc_list.push_back(BoundaryConditionHandler(bc));
+ }
+ }
+ }
+ break;
+ }
case IBoundaryConditionDescriptor::Type::pressure: {
const PressureBoundaryConditionDescriptor& pressure_bc_descriptor =
dynamic_cast<const PressureBoundaryConditionDescriptor&>(*bc_descriptor);
@@ -172,7 +217,9 @@ struct GlaceScheme
break;
}
default: {
- throw UnexpectedError("Unknown BCDescription\n");
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver";
+ throw NormalError(error_msg.str());
}
}
}
@@ -325,6 +372,19 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("velocity",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+ const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
+ const FunctionSymbolId&)>>(
+
+ [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+ const FunctionSymbolId& velocity_id)
+ -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+ return std::make_shared<VelocityBoundaryConditionDescriptor>(boundary, velocity_id);
+ }
+
+ ));
+
this->_addBuiltinFunction("glace",
std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>,
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 75b8eec54..1fcfa105e 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -123,11 +123,22 @@ class AcousticSolver
{
for (const auto& handler : m_boundary_condition_list) {
switch (handler.boundaryCondition().type()) {
- case BoundaryCondition::normal_velocity: {
- throw NotImplementedError("normal_velocity BC");
- }
case BoundaryCondition::velocity: {
- throw NotImplementedError("velocity BC");
+ 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 =
@@ -142,16 +153,17 @@ class AcousticSolver
const auto& node_list = pressure_bc.faceList();
const auto& value_list = pressure_bc.valueList();
- for (size_t i_node = 0; i_node < node_list.size(); ++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);
+ 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);
+ 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);
- }
+ 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();
diff --git a/src/scheme/BoundaryCondition.hpp b/src/scheme/BoundaryCondition.hpp
index cce2f0636..2bd924974 100644
--- a/src/scheme/BoundaryCondition.hpp
+++ b/src/scheme/BoundaryCondition.hpp
@@ -15,7 +15,6 @@ class BoundaryCondition // clazy:exclude=copyable-polymorphic
enum Type
{
velocity,
- normal_velocity,
pressure,
symmetry
};
@@ -58,7 +57,7 @@ class PressureBoundaryCondition final : public BoundaryCondition // clazy:excl
return m_value_list;
}
- PressureBoundaryCondition(const Array<const FaceId>& face_list, Array<const double> 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}
{}
@@ -85,7 +84,7 @@ class PressureBoundaryCondition<1> final : public BoundaryCondition // clazy:e
return m_value_list;
}
- PressureBoundaryCondition(const Array<const NodeId>& face_list, Array<const double> 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}
{}
@@ -129,6 +128,33 @@ class SymmetryBoundaryCondition : public BoundaryCondition
~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:
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index 8abb25371..95f68bd25 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -8,6 +8,7 @@ class IBoundaryConditionDescriptor
public:
enum class Type
{
+ velocity,
pressure,
symmetry
};
diff --git a/src/scheme/VelocityBoundaryConditionDescriptor.hpp b/src/scheme/VelocityBoundaryConditionDescriptor.hpp
new file mode 100644
index 000000000..5ce10e6a5
--- /dev/null
+++ b/src/scheme/VelocityBoundaryConditionDescriptor.hpp
@@ -0,0 +1,55 @@
+#ifndef VELOCITY_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define VELOCITY_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class VelocityBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+ std::ostream&
+ _write(std::ostream& os) const final
+ {
+ os << "velocity(" << *m_boundary_descriptor << ")";
+ return os;
+ }
+
+ std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+ const FunctionSymbolId m_function_symbol_id;
+
+ public:
+ FunctionSymbolId
+ functionSymbolId() const
+ {
+ return m_function_symbol_id;
+ }
+
+ const IBoundaryDescriptor&
+ boundaryDescriptor() const
+ {
+ return *m_boundary_descriptor;
+ }
+
+ Type
+ type() const final
+ {
+ return Type::velocity;
+ }
+
+ VelocityBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+ const FunctionSymbolId& function_symbol_id)
+ : m_boundary_descriptor(boundary_descriptor), m_function_symbol_id{function_symbol_id}
+ {
+ ;
+ }
+
+ VelocityBoundaryConditionDescriptor(const VelocityBoundaryConditionDescriptor&) = delete;
+ VelocityBoundaryConditionDescriptor(VelocityBoundaryConditionDescriptor&&) = delete;
+
+ ~VelocityBoundaryConditionDescriptor() = default;
+};
+
+#endif // VELOCITY!_BOUNDARY_CONDITION_DESCRIPTOR_HPP
--
GitLab