diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index c329feb029f63d3a24c7695c2ecd84556a4fe19b..5cf0b0c5d6cdc73719a38b12a1e235933afbbf0f 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -15,7 +15,12 @@
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
#include <output/VTKWriter.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/FourierBoundaryConditionDescriptor.hpp>
+#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
template <size_t Dimension>
HeatDiamondScheme<Dimension>::HeatDiamondScheme(
@@ -30,9 +35,82 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
using MeshType = Mesh<ConnectivityType>;
using MeshDataType = MeshData<Dimension>;
+ constexpr ItemType FaceType = [] {
+ if constexpr (Dimension > 1) {
+ return ItemType::face;
+ } else {
+ return ItemType::node;
+ }
+ }();
+
+ using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition,
+ SymmetryBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ BoundaryConditionList boundary_condition_list;
+
std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
std::shared_ptr m_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ // const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+ // dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+ throw NotImplementedError("NIY");
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::dirichlet: {
+ const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+ dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*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 == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ MeshNodeBoundary<Dimension> mesh_node_boundary{m_mesh, ref_face_list};
+
+ const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ const auto& node_list = mesh_node_boundary.nodeList();
+
+ Array<const double> value_list =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id,
+ m_mesh->xr(),
+ node_list);
+
+ boundary_condition_list.push_back(DirichletBoundaryCondition{node_list, value_list});
+ }
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::fourier: {
+ // const FourierBoundaryConditionDescriptor& fourier_bc_descriptor =
+ // dynamic_cast<const FourierBoundaryConditionDescriptor&>(*bc_descriptor);
+ throw NotImplementedError("NIY");
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::neumann: {
+ // const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor =
+ // dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_descriptor);
+ throw NotImplementedError("NIY");
+ break;
+ }
+ default: {
+ is_valid_boundary_condition = false;
+ }
+ }
+ if (not is_valid_boundary_condition) {
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for heat equation";
+ throw NormalError(error_msg.str());
+ }
+ }
+
if constexpr (Dimension > 1) {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
@@ -345,8 +423,9 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
Vector<double> error{m_mesh->numberOfCells()};
parallel_for(
- m_mesh->numberOfCells(),
- PUGS_LAMBDA(CellId cell_id) { error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * primal_Vj[cell_id]; });
+ m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]);
+ });
std::cout << "||Error||_2 = " << std::sqrt((error, error)) << "\n";
@@ -365,6 +444,54 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
}
+template <size_t Dimension>
+class HeatDiamondScheme<Dimension>::DirichletBoundaryCondition
+{
+ private:
+ const Array<const double> m_value_list;
+ const Array<const NodeId> m_node_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_node_list;
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ DirichletBoundaryCondition(const Array<const NodeId>& node_list, const Array<const double>& value_list)
+ : m_value_list{value_list}, m_node_list{node_list}
+ {}
+
+ ~DirichletBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class HeatDiamondScheme<Dimension>::FourierBoundaryCondition
+{
+ public:
+ ~FourierBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class HeatDiamondScheme<Dimension>::NeumannBoundaryCondition
+{
+ public:
+ ~NeumannBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class HeatDiamondScheme<Dimension>::SymmetryBoundaryCondition
+{
+ public:
+ ~SymmetryBoundaryCondition() = default;
+};
+
template HeatDiamondScheme<1>::HeatDiamondScheme(
std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.hpp b/src/language/algorithms/HeatDiamondAlgorithm.hpp
index fc8885d618f4aa5071ab74e640cad689e9314d71..c3b3d2e84ff6e13c9fefb32303e7eb2d4b564c01 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.hpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.hpp
@@ -6,11 +6,19 @@
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <memory>
+#include <variant>
#include <vector>
template <size_t Dimension>
-struct HeatDiamondScheme
+class HeatDiamondScheme
{
+ private:
+ class DirichletBoundaryCondition;
+ class FourierBoundaryCondition;
+ class NeumannBoundaryCondition;
+ class SymmetryBoundaryCondition;
+
+ public:
HeatDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& T_id,
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index a30d2a252807c9620e25a93a1a7f0e01dbb8edf5..d2d25382cbb376c6fab9ae67c8044bf719e47b42 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -83,6 +83,44 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("dirichlet",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+ const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
+ const FunctionSymbolId&)>>(
+
+ [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+ const FunctionSymbolId& g_id) -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+ return std::make_shared<DirichletBoundaryConditionDescriptor>("dirichlet", boundary,
+ g_id);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("fourier",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+ const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
+ const FunctionSymbolId&, const FunctionSymbolId&)>>(
+
+ [](std::shared_ptr<const IBoundaryDescriptor> boundary, const FunctionSymbolId& alpha_id,
+ const FunctionSymbolId& g_id) -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+ return std::make_shared<FourierBoundaryConditionDescriptor>("fourier", boundary,
+ alpha_id, g_id);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("neumann",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+ const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
+ const FunctionSymbolId&)>>(
+
+ [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+ const FunctionSymbolId& g_id) -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+ return std::make_shared<NeumannBoundaryConditionDescriptor>("neumann", boundary, g_id);
+ }
+
+ ));
+
this->_addBuiltinFunction("glace",
std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>,