Prepare use of boundary conditions comming from BC list

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>>&,

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,

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>,