Nouveaux fichiers pour solveur nodal

src/scheme/ScalarNodalScheme.cpp

+#include <scheme/ScalarNodalScheme.hpp>
+
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+
+template <size_t Dimension>
+class ScalarNodalSchemeHandler::InterpolationWeightsManager
+{
+ private:
+  std::shared_ptr<const Mesh<Connectivity<Dimension>>> m_mesh;
+  FaceValue<bool> m_primal_face_is_on_boundary;
+  NodeValue<bool> m_primal_node_is_on_boundary;
+  CellValuePerNode<double> m_w_rj;
+  FaceValuePerNode<double> m_w_rl;
+
+ public:
+  CellValuePerNode<double>&
+  wrj()
+  {
+    return m_w_rj;
+  }
+
+  FaceValuePerNode<double>&
+  wrl()
+  {
+    return m_w_rl;
+  }
+
+  void
+  compute()
+  {
+    using MeshDataType      = MeshData<Dimension>;
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+    const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
+
+    const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+    const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
+    const auto& node_to_cell_matrix                  = m_mesh->connectivity().nodeToCellMatrix();
+    const auto& node_to_face_matrix                  = m_mesh->connectivity().nodeToFaceMatrix();
+    CellValuePerNode<double> w_rj{m_mesh->connectivity()};
+    FaceValuePerNode<double> w_rl{m_mesh->connectivity()};
+
+    for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
+      w_rl[i] = std::numeric_limits<double>::signaling_NaN();
+    }
+
+    for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) {
+      SmallVector<double> b{Dimension + 1};
+      b[0] = 1;
+      for (size_t i = 1; i < Dimension + 1; i++) {
+        b[i] = xr[i_node][i - 1];
+      }
+      const auto& node_to_cell = node_to_cell_matrix[i_node];
+
+      if (not m_primal_node_is_on_boundary[i_node]) {
+        SmallMatrix<double> A{Dimension + 1, node_to_cell.size()};
+        for (size_t j = 0; j < node_to_cell.size(); j++) {
+          A(0, j) = 1;
+        }
+        for (size_t i = 1; i < Dimension + 1; i++) {
+          for (size_t j = 0; j < node_to_cell.size(); j++) {
+            const CellId J = node_to_cell[j];
+            A(i, j)        = xj[J][i - 1];
+          }
+        }
+        SmallVector<double> x{node_to_cell.size()};
+        x = zero;
+
+        LeastSquareSolver ls_solver;
+        ls_solver.solveLocalSystem(A, x, b);
+
+        for (size_t j = 0; j < node_to_cell.size(); j++) {
+          w_rj(i_node, j) = x[j];
+        }
+      } else {
+        int nb_face_used = 0;
+        for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+          FaceId face_id = node_to_face_matrix[i_node][i_face];
+          if (m_primal_face_is_on_boundary[face_id]) {
+            nb_face_used++;
+          }
+        }
+        SmallMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
+        for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) {
+          A(0, j) = 1;
+        }
+        for (size_t i = 1; i < Dimension + 1; i++) {
+          for (size_t j = 0; j < node_to_cell.size(); j++) {
+            const CellId J = node_to_cell[j];
+            A(i, j)        = xj[J][i - 1];
+          }
+        }
+        for (size_t i = 1; i < Dimension + 1; i++) {
+          int cpt_face = 0;
+          for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+            FaceId face_id = node_to_face_matrix[i_node][i_face];
+            if (m_primal_face_is_on_boundary[face_id]) {
+              A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
+              cpt_face++;
+            }
+          }
+        }
+
+        SmallVector<double> x{node_to_cell.size() + nb_face_used};
+        x = zero;
+
+        LeastSquareSolver ls_solver;
+        ls_solver.solveLocalSystem(A, x, b);
+
+        for (size_t j = 0; j < node_to_cell.size(); j++) {
+          w_rj(i_node, j) = x[j];
+        }
+        int cpt_face = node_to_cell.size();
+        for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+          FaceId face_id = node_to_face_matrix[i_node][i_face];
+          if (m_primal_face_is_on_boundary[face_id]) {
+            w_rl(i_node, i_face) = x[cpt_face++];
+          }
+        }
+      }
+    }
+    m_w_rj = w_rj;
+    m_w_rl = w_rl;
+  }
+
+  InterpolationWeightsManager(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh,
+                              FaceValue<bool> primal_face_is_on_boundary,
+                              NodeValue<bool> primal_node_is_on_boundary)
+    : m_mesh(mesh),
+      m_primal_face_is_on_boundary(primal_face_is_on_boundary),
+      m_primal_node_is_on_boundary(primal_node_is_on_boundary)
+  {}
+
+  ~InterpolationWeightsManager() = default;
+};
+
+class ScalarNodalSchemeHandler::IScalarNodalScheme
+{
+ public:
+  virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0;
+
+  IScalarNodalScheme()          = default;
+  virtual ~IScalarNodalScheme() = default;
+};
+
+template <size_t Dimension>
+class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHandler::IScalarNodalScheme
+{
+ private:
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<ConnectivityType>;
+  using MeshDataType     = MeshData<Dimension>;
+
+  std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_solution;
+
+  class DirichletBoundaryCondition
+  {
+   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;
+    }
+
+    DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
+      : m_value_list{value_list}, m_face_list{face_list}
+    {
+      Assert(m_value_list.size() == m_face_list.size());
+    }
+
+    ~DirichletBoundaryCondition() = default;
+  };
+
+  class NeumannBoundaryCondition
+  {
+   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;
+    }
+
+    NeumannBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
+      : m_value_list{value_list}, m_face_list{face_list}
+    {
+      Assert(m_value_list.size() == m_face_list.size());
+    }
+
+    ~NeumannBoundaryCondition() = default;
+  };
+
+  class FourierBoundaryCondition
+  {
+   private:
+    const Array<const double> m_coef_list;
+    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;
+    }
+
+    const Array<const double>&
+    coefList() const
+    {
+      return m_coef_list;
+    }
+
+   public:
+    FourierBoundaryCondition(const Array<const FaceId>& face_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}
+    {
+      Assert(m_coef_list.size() == m_face_list.size());
+      Assert(m_value_list.size() == m_face_list.size());
+    }
+
+    ~FourierBoundaryCondition() = default;
+  };
+
+  class SymmetryBoundaryCondition
+  {
+   private:
+    const Array<const FaceId> m_face_list;
+
+   public:
+    const Array<const FaceId>&
+    faceList() const
+    {
+      return m_face_list;
+    }
+
+   public:
+    SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
+
+    ~SymmetryBoundaryCondition() = default;
+  };
+
+ public:
+  std::shared_ptr<const IDiscreteFunction>
+  getSolution() const final
+  {
+    return m_solution;
+  }
+
+  ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
+                      const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
+                      const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mub,
+                      const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mu,
+                      const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
+                      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+  {
+    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mu->mesh(),
+           "diffusion coefficient is not defined on the dual mesh!");
+    Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mub->mesh(),
+           "boundary diffusion coefficient is not defined on the dual mesh!");
+
+    using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition,
+                                           NeumannBoundaryCondition, SymmetryBoundaryCondition>;
+
+    using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+    BoundaryConditionList boundary_condition_list;
+
+    for (const auto& bc_descriptor : bc_descriptor_list) {
+      bool is_valid_boundary_condition = true;
+
+      switch (bc_descriptor->type()) {
+      case IBoundaryConditionDescriptor::Type::symmetry: {
+        throw NotImplementedError("NIY");
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::dirichlet: {
+        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+        if constexpr (Dimension > 1) {
+          MeshFaceBoundary<Dimension> mesh_face_boundary =
+            getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor());
+
+          const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+          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());
+
+          boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), value_list});
+
+        } else {
+          throw NotImplementedError("Dirichlet BC in 1d");
+        }
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::fourier: {
+        throw NotImplementedError("NIY");
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::neumann: {
+        const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor =
+          dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_descriptor);
+
+        if constexpr (Dimension > 1) {
+          MeshFaceBoundary<Dimension> mesh_face_boundary =
+            getMeshFaceBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
+
+          const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
+          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());
+
+          boundary_condition_list.push_back(NeumannBoundaryCondition{mesh_face_boundary.faceList(), value_list});
+
+        } else {
+          throw NotImplementedError("Dirichlet BC in 1d");
+        }
+        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) {
+      const CellValue<const size_t> cell_dof_number = [&] {
+        CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
+        return compute_cell_dof_number;
+      }();
+      size_t number_of_dof = mesh->numberOfCells();
+
+      const FaceValue<const size_t> face_dof_number = [&] {
+        FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
+        compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
+        for (const auto& boundary_condition : boundary_condition_list) {
+          std::visit(
+            [&](auto&& bc) {
+              using T = std::decay_t<decltype(bc)>;
+              if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
+                            (std::is_same_v<T, FourierBoundaryCondition>) or
+                            (std::is_same_v<T, SymmetryBoundaryCondition>) or
+                            (std::is_same_v<T, DirichletBoundaryCondition>)) {
+                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];
+                  if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
+                    std::ostringstream os;
+                    os << "The face " << face_id << " is used at least twice for boundary conditions";
+                    throw NormalError(os.str());
+                  } else {
+                    compute_face_dof_number[face_id] = number_of_dof++;
+                  }
+                }
+              }
+            },
+            boundary_condition);
+        }
+
+        return compute_face_dof_number;
+      }();
+
+      const FaceValue<const bool> primal_face_is_neumann = [&] {
+        FaceValue<bool> face_is_neumann{mesh->connectivity()};
+        face_is_neumann.fill(false);
+        for (const auto& boundary_condition : boundary_condition_list) {
+          std::visit(
+            [&](auto&& bc) {
+              using T = std::decay_t<decltype(bc)>;
+              if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
+                            (std::is_same_v<T, FourierBoundaryCondition>) or
+                            (std::is_same_v<T, SymmetryBoundaryCondition>)) {
+                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];
+                  face_is_neumann[face_id] = true;
+                }
+              }
+            },
+            boundary_condition);
+        }
+
+        return face_is_neumann;
+      }();
+
+      const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+      const auto& face_to_cell_matrix        = mesh->connectivity().faceToCellMatrix();
+      NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
+      if (parallel::size() > 1) {
+        throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
+      }
+
+      primal_node_is_on_boundary.fill(false);
+      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+        if (face_to_cell_matrix[face_id].size() == 1) {
+          for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
+            NodeId node_id                      = primal_face_to_node_matrix[face_id][i_node];
+            primal_node_is_on_boundary[node_id] = true;
+          }
+        }
+      }
+
+      FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
+      if (parallel::size() > 1) {
+        throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
+      }
+
+      primal_face_is_on_boundary.fill(false);
+      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+        if (face_to_cell_matrix[face_id].size() == 1) {
+          primal_face_is_on_boundary[face_id] = true;
+        }
+      }
+
+      FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
+      if (parallel::size() > 1) {
+        throw NotImplementedError("Calculation of face_is_neumann is incorrect");
+      }
+
+      primal_face_is_dirichlet.fill(false);
+      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+        primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]));
+      }
+
+      InterpolationWeightsManager iwm(mesh, primal_face_is_on_boundary, primal_node_is_on_boundary);
+      iwm.compute();
+      CellValuePerNode<double> w_rj = iwm.wrj();
+      FaceValuePerNode<double> w_rl = iwm.wrl();
+
+      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+      const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+      const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
+      const auto& node_to_face_matrix                  = mesh->connectivity().nodeToFaceMatrix();
+
+      {
+        std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh);
+
+        MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
+
+        std::shared_ptr mapper =
+          DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity());
+
+        CellValue<const double> dual_kappaj  = dual_mu->cellValues();
+        CellValue<const double> dual_kappajb = dual_mub->cellValues();
+
+        const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
+
+        const FaceValue<const double> mes_l = [&] {
+          if constexpr (Dimension == 1) {
+            FaceValue<double> compute_mes_l{mesh->connectivity()};
+            compute_mes_l.fill(1);
+            return compute_mes_l;
+          } else {
+            return mesh_data.ll();
+          }
+        }();
+
+        const CellValue<const double> dual_mes_l_j = [=] {
+          CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
+          mapper->toDualCell(mes_l, compute_mes_j);
+
+          return compute_mes_j;
+        }();
+
+        FaceValue<const CellId> face_dual_cell_id = [=]() {
+          FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
+          CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
+          parallel_for(
+            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
+
+          mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id);
+
+          return computed_face_dual_cell_id;
+        }();
+
+        NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
+          CellValue<NodeId> cell_ignored_id{mesh->connectivity()};
+          cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
+
+          NodeValue<NodeId> node_primal_id{mesh->connectivity()};
+
+          parallel_for(
+            mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
+
+          NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
+
+          mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
+
+          return computed_dual_node_primal_node_id;
+        }();
+
+        CellValue<NodeId> primal_cell_dual_node_id = [=]() {
+          CellValue<NodeId> cell_id{mesh->connectivity()};
+          NodeValue<NodeId> node_ignored_id{mesh->connectivity()};
+          node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
+
+          NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()};
+
+          parallel_for(
+            diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; });
+
+          CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()};
+
+          mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id);
+
+          return cell_id;
+        }();
+        const auto& dual_Cjr                     = diamond_mesh_data.Cjr();
+        FaceValue<TinyVector<Dimension>> dualClj = [&] {
+          FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()};
+          const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix();
+          const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
+          parallel_for(
+            mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+              const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+              for (size_t i = 0; i < primal_face_to_cell.size(); i++) {
+                CellId cell_id            = primal_face_to_cell[i];
+                const NodeId dual_node_id = primal_cell_dual_node_id[cell_id];
+                for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size();
+                     i_dual_cell++) {
+                  const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell];
+                  if (face_dual_cell_id[face_id] == dual_cell_id) {
+                    for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
+                         i_dual_node++) {
+                      const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
+                      if (final_dual_node_id == dual_node_id) {
+                        computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node);
+                      }
+                    }
+                  }
+                }
+              }
+            });
+          return computedClj;
+        }();
+
+        FaceValue<TinyVector<Dimension>> nlj = [&] {
+          FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()};
+          parallel_for(
+            mesh->numberOfFaces(),
+            PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; });
+          return computedNlj;
+        }();
+
+        FaceValue<const double> alpha_l = [&] {
+          CellValue<double> alpha_j{diamond_mesh->connectivity()};
+
+          parallel_for(
+            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
+              alpha_j[diamond_cell_id] = dual_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
+            });
+
+          FaceValue<double> computed_alpha_l{mesh->connectivity()};
+          mapper->fromDualCell(alpha_j, computed_alpha_l);
+          return computed_alpha_l;
+        }();
+
+        FaceValue<const double> alphab_l = [&] {
+          CellValue<double> alpha_jb{diamond_mesh->connectivity()};
+
+          parallel_for(
+            diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
+              alpha_jb[diamond_cell_id] = dual_kappajb[diamond_cell_id] / dual_Vj[diamond_cell_id];
+            });
+
+          FaceValue<double> computed_alpha_lb{mesh->connectivity()};
+          mapper->fromDualCell(alpha_jb, computed_alpha_lb);
+          return computed_alpha_lb;
+        }();
+
+        const CellValue<const double> primal_Vj = mesh_data.Vj();
+
+        const Array<int> non_zeros{number_of_dof};
+        non_zeros.fill(Dimension);
+        CRSMatrixDescriptor<double> S(number_of_dof, number_of_dof, non_zeros);
+
+        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+          const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+          const double beta_l             = l2Norm(dualClj[face_id]) * alpha_l[face_id] * mes_l[face_id];
+          const double betab_l            = l2Norm(dualClj[face_id]) * alphab_l[face_id] * mes_l[face_id];
+          for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
+            const CellId cell_id1 = primal_face_to_cell[i_cell];
+            for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
+              const CellId cell_id2 = primal_face_to_cell[j_cell];
+              if (i_cell == j_cell) {
+                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
+                if (primal_face_is_neumann[face_id]) {
+                  S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= betab_l;
+                }
+              } else {
+                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
+              }
+            }
+          }
+        }
+
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          const size_t j = cell_dof_number[cell_id];
+          S(j, j) += (*alpha)[cell_id] * primal_Vj[cell_id];
+        }
+
+        const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
+
+        const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+          const double alpha_face_id  = mes_l[face_id] * alpha_l[face_id];
+          const double alphab_face_id = mes_l[face_id] * alphab_l[face_id];
+
+          for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
+            CellId i_id                            = face_to_cell_matrix[face_id][i_face_cell];
+            const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
+
+            for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
+              NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
+
+              const TinyVector<Dimension> nil = [&] {
+                if (is_face_reversed_for_cell_i) {
+                  return -nlj[face_id];
+                } else {
+                  return nlj[face_id];
+                }
+              }();
+
+              CellId dual_cell_id = face_dual_cell_id[face_id];
+
+              for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
+                const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
+                if (dual_node_primal_node_id[dual_node_id] == node_id) {
+                  const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
+
+                  const double a_ir  = alpha_face_id * dot(nil, Clr);
+                  const double ab_ir = alphab_face_id * dot(nil, Clr);
+
+                  for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
+                    CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
+                    S(cell_dof_number[i_id], cell_dof_number[j_id]) -= w_rj(node_id, j_cell) * a_ir;
+                    if (primal_face_is_neumann[face_id]) {
+                      S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * ab_ir;
+                    }
+                  }
+                  if (primal_node_is_on_boundary[node_id]) {
+                    for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
+                      FaceId l_id = node_to_face_matrix[node_id][l_face];
+                      if (primal_face_is_on_boundary[l_id]) {
+                        S(cell_dof_number[i_id], face_dof_number[l_id]) -= w_rl(node_id, l_face) * a_ir;
+                        if (primal_face_is_neumann[face_id]) {
+                          S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * ab_ir;
+                        }
+                      }
+                    }
+                  }
+                }
+              }
+            }
+          }
+        }
+        for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+          if (primal_face_is_dirichlet[face_id]) {
+            S(face_dof_number[face_id], face_dof_number[face_id]) += 1;
+          }
+        }
+
+        CellValue<const double> fj = f->cellValues();
+
+        CRSMatrix A{S.getCRSMatrix()};
+        Vector<double> b{number_of_dof};
+        b = zero;
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          b[cell_dof_number[cell_id]] = fj[cell_id] * primal_Vj[cell_id];
+        }
+        // Dirichlet on b^L_D
+        {
+          for (const auto& boundary_condition : boundary_condition_list) {
+            std::visit(
+              [&](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& 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];
+                    b[face_dof_number[face_id]] += value_list[i_face];
+                  }
+                }
+              },
+              boundary_condition);
+          }
+        }
+        // EL b^L
+        for (const auto& boundary_condition : boundary_condition_list) {
+          std::visit(
+            [&](auto&& bc) {
+              using T = std::decay_t<decltype(bc)>;
+              if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
+                            (std::is_same_v<T, FourierBoundaryCondition>)) {
+                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) {
+                  FaceId face_id = face_list[i_face];
+                  b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face];
+                }
+              }
+            },
+            boundary_condition);
+        }
+
+        Vector<double> T{number_of_dof};
+        T = zero;
+
+        LinearSolver solver;
+        solver.solveLocalSystem(A, T, b);
+
+        m_solution     = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
+        auto& solution = *m_solution;
+        parallel_for(
+          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_dof_number[cell_id]]; });
+      }
+    }
+  }
+};
+
+std::shared_ptr<const IDiscreteFunction>
+ScalarNodalSchemeHandler::solution() const
+{
+  return m_scheme->getSolution();
+}
+
+ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
+  const std::shared_ptr<const IDiscreteFunction>& alpha,
+  const std::shared_ptr<const IDiscreteFunction>& dual_mub,
+  const std::shared_ptr<const IDiscreteFunction>& dual_mu,
+  const std::shared_ptr<const IDiscreteFunction>& f,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+{
+  const std::shared_ptr i_mesh = getCommonMesh({alpha, f});
+  if (not i_mesh) {
+    throw NormalError("primal discrete functions are not defined on the same mesh");
+  }
+  const std::shared_ptr i_dual_mesh = getCommonMesh({dual_mub, dual_mu});
+  if (not i_dual_mesh) {
+    throw NormalError("dual discrete functions are not defined on the same mesh");
+  }
+  checkDiscretizationType({alpha, dual_mub, dual_mu, f}, DiscreteFunctionType::P0);
+
+  switch (i_mesh->dimension()) {
+  case 1: {
+    using MeshType                   = Mesh<Connectivity<1>>;
+    using DiscreteScalarFunctionType = DiscreteFunctionP0<1, double>;
+
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
+    }
+
+    m_scheme =
+      std::make_unique<ScalarNodalScheme<1>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
+                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
+                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
+                                               bc_descriptor_list);
+    break;
+  }
+  case 2: {
+    using MeshType                   = Mesh<Connectivity<2>>;
+    using DiscreteScalarFunctionType = DiscreteFunctionP0<2, double>;
+
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
+    }
+
+    m_scheme =
+      std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
+                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
+                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
+                                               bc_descriptor_list);
+    break;
+  }
+  case 3: {
+    using MeshType                   = Mesh<Connectivity<3>>;
+    using DiscreteScalarFunctionType = DiscreteFunctionP0<3, double>;
+
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+    if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
+      throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
+    }
+
+    m_scheme =
+      std::make_unique<ScalarNodalScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
+                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
+                                               std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
+                                               bc_descriptor_list);
+    break;
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+ScalarNodalSchemeHandler::~ScalarNodalSchemeHandler() = default;