New boundary condition treatment for Neumann

src/language/algorithms/HeatDiamondAlgorithm2.cpp

+#include <language/algorithms/HeatDiamondAlgorithm2.hpp>
+
+#include <algebra/BiCGStab.hpp>
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/LocalRectangularMatrix.hpp>
+#include <algebra/PCG.hpp>
+#include <algebra/SparseMatrixDescriptor.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
+#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>
+HeatDiamondScheme2<Dimension>::HeatDiamondScheme2(
+  std::shared_ptr<const IMesh> i_mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const FunctionSymbolId& T_id,
+  const FunctionSymbolId& kappa_id,
+  const FunctionSymbolId& f_id)
+{
+  using ConnectivityType = Connectivity<Dimension>;
+  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 mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+  NodeValue<bool> is_dirichlet{mesh->connectivity()};
+  is_dirichlet.fill(false);
+  NodeValue<double> dirichlet_value{mesh->connectivity()};
+  dirichlet_value.fill(std::numeric_limits<double>::signaling_NaN());
+
+  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 < mesh->connectivity().template numberOfRefItemList<FaceType>();
+           ++i_ref_face_list) {
+        const auto& ref_face_list = 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{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, mesh->xr(),
+                                                                                                      node_list);
+
+          for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+            NodeId node_id = node_list[i_node];
+            if (not is_dirichlet[node_id]) {
+              is_dirichlet[node_id]    = true;
+              dirichlet_value[node_id] = value_list[i_node];
+            }
+          }
+
+          boundary_condition_list.push_back(DirichletBoundaryCondition{node_list, value_list});
+        }
+      }
+
+      if constexpr (Dimension > 1) {
+        for (size_t i_ref_node_list = 0;
+             i_ref_node_list < mesh->connectivity().template numberOfRefItemList<ItemType::node>(); ++i_ref_node_list) {
+          const auto& ref_node_list = mesh->connectivity().template refItemList<ItemType::node>(i_ref_node_list);
+          const RefId& ref          = ref_node_list.refId();
+          if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+            const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+            const auto& node_list = ref_node_list.list();
+
+            Array<const double> value_list =
+              InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id,
+                                                                                                        mesh->xr(),
+                                                                                                        node_list);
+
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              NodeId node_id = node_list[i_node];
+              if (not is_dirichlet[node_id]) {
+                is_dirichlet[node_id]    = true;
+                dirichlet_value[node_id] = value_list[i_node];
+              }
+            }
+
+            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);
+
+      for (size_t i_ref_face_list = 0; i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>();
+           ++i_ref_face_list) {
+        const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+        const RefId& ref          = ref_face_list.refId();
+        if (ref == neumann_bc_descriptor.boundaryDescriptor()) {
+          const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
+
+          if constexpr (Dimension > 1) {
+            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+            Array<const FaceId> face_list = ref_face_list.list();
+            Array<const double> value_list =
+              InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id,
+                                                                                                        mesh_data.xl(),
+                                                                                                        face_list);
+            boundary_condition_list.push_back(NeumannBoundaryCondition{face_list, value_list});
+          } else {
+            throw NotImplementedError("Neumann conditions are not supported 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>)) {
+              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 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_neumann(mesh->connectivity());
+    if (parallel::size() > 1) {
+      throw NotImplementedError("Calculation of face_is_neumann is incorrect");
+    }
+
+    primal_face_is_neumann.fill(false);
+    for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+      if (primal_face_is_on_boundary[face_id]) {
+        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];
+          if (not is_dirichlet[node_id]) {
+            primal_face_is_neumann[face_id] = true;
+          }
+        }
+      }
+    }
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+    CellValue<double> Tj =
+      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
+
+    NodeValue<double> Tr(mesh->connectivity());
+    const NodeValue<const TinyVector<Dimension>>& xr = 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                  = mesh->connectivity().nodeToCellMatrix();
+    const auto& node_to_face_matrix                  = mesh->connectivity().nodeToFaceMatrix();
+    CellValuePerNode<double> w_rj{mesh->connectivity()};
+    FaceValuePerNode<double> w_rl{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 < mesh->numberOfNodes(); ++i_node) {
+      Vector<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 primal_node_is_on_boundary[i_node]) {
+        LocalRectangularMatrix<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];
+          }
+        }
+        Vector<double> x{node_to_cell.size()};
+        x = 0;
+
+        LocalRectangularMatrix B = transpose(A) * A;
+        Vector y                 = transpose(A) * b;
+        PCG<false>{y, B, B, x, 10, 1e-12};
+
+        Tr[i_node] = 0;
+        for (size_t j = 0; j < node_to_cell.size(); j++) {
+          Tr[i_node] += x[j] * Tj[node_to_cell[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 (primal_face_is_on_boundary[face_id]) {
+            nb_face_used++;
+          }
+        }
+        LocalRectangularMatrix<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 (primal_face_is_on_boundary[face_id]) {
+              A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
+              cpt_face++;
+            }
+          }
+        }
+
+        Vector<double> x{node_to_cell.size() + nb_face_used};
+        x = 0;
+
+        LocalRectangularMatrix B = transpose(A) * A;
+        Vector y                 = transpose(A) * b;
+        PCG<false>{y, B, B, x, 10, 1e-12};
+
+        Tr[i_node] = 0;
+        for (size_t j = 0; j < node_to_cell.size(); j++) {
+          Tr[i_node] += x[j] * Tj[node_to_cell[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 (primal_face_is_on_boundary[face_id]) {
+            w_rl(i_node, i_face) = x[cpt_face++];
+          }
+        }
+      }
+    }
+
+    {
+      VTKWriter vtk_writer("T_" + std::to_string(Dimension), 0.01);
+
+      vtk_writer.write(mesh,
+                       {NamedItemValue{"Tr", Tr}, NamedItemValue{"temperature", Tj},
+                        NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+                        NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                        NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
+                       0, true);   // forces last output
+    }
+
+    {
+      std::shared_ptr diamond_mesh = DiamondDualMeshManager::instance().getDiamondDualMesh(mesh);
+
+      MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
+
+      std::shared_ptr mapper =
+        DiamondDualConnectivityManager::instance().getConnectivityToDiamondDualConnectivityDataMapper(
+          mesh->connectivity());
+
+      NodeValue<double> Trd{diamond_mesh->connectivity()};
+
+      mapper->toDualNode(Tr, Tj, Trd);
+
+      CellValue<double> kappaj =
+        InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id,
+                                                                                                  mesh_data.xj());
+
+      CellValue<double> dual_kappaj =
+        InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id,
+                                                                                                  diamond_mesh_data
+                                                                                                    .xj());
+
+      VTKWriter vtk_writer("D_" + std::to_string(Dimension), 0.01);
+
+      vtk_writer.write(diamond_mesh,
+                       {NamedItemValue{"kappaj", dual_kappaj}, NamedItemValue{"coords", diamond_mesh->xr()},
+                        NamedItemValue{"Vj", diamond_mesh_data.Vj()}, NamedItemValue{"xj", diamond_mesh_data.xj()}},
+                       0, true);   // forces last output
+
+      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 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_mes_l_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);
+
+        VTKWriter vtk_writer("DD_" + std::to_string(Dimension), 0.01);
+
+        vtk_writer.write(diamond_mesh,
+                         {NamedItemValue{"alphaj", alpha_j}, NamedItemValue{"xj", diamond_mesh_data.xj()},
+                          NamedItemValue{"Vl", diamond_mesh_data.Vj()}, NamedItemValue{"|l|", dual_mes_l_j}},
+                         0,
+                         true);   // forces last output
+
+        return computed_alpha_l;
+      }();
+
+      SparseMatrixDescriptor S{number_of_dof};
+      // A^{JJ} & A^{LJ}
+      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+        // if (not primal_face_is_neumann[face_id]) { COMMENTE
+        const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+        const double beta_l             = 1. / Dimension * alpha_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]) += beta_l;   // AJOUT EL
+              }
+            } else {
+              S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
+            }
+          }
+        }
+      }
+
+      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;
+      }();
+
+      const auto& dual_Cjr = diamond_mesh_data.Cjr();
+
+      const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
+
+      const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
+      const auto& primal_face_cell_is_reversed                       = mesh->connectivity().cellFaceIsReversed();
+      const auto& face_local_numbers_in_their_cells = mesh->connectivity().faceLocalNumbersInTheirCells();
+      const auto& primal_node_to_cell_matrix        = mesh->connectivity().nodeToCellMatrix();
+      // A^{JR}A^{RJ} & A^{JR}A^{RL} & A^{LR}A^{RL}
+      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+        //        if (face_to_cell_matrix[face_id].size() > 1) { COMMENTE
+        const double alpha_face_id = alpha_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 =
+            primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_face_cell));
+
+          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];
+
+            if (not is_dirichlet[node_id]) {
+              const TinyVector<Dimension> nil = [&] {
+                if (is_face_reversed_for_cell_i) {
+                  return -primal_nlr(face_id, i_node);
+                } else {
+                  return primal_nlr(face_id, i_node);
+                }
+              }();
+
+              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 * (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) * a_ir;   // AJOUT
+                    }
+                  }
+                  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_neumann[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) * a_ir;   // AJOUT
+                        }
+                        // S(face_dof_number[l_id], cell_dof_number[i_id]) += w_rl(node_id, l_face) * a_ir;
+                      }
+                    }
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+      // COMMENTE POUR L INSTANT PAS DE FOURIER
+      // for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+      //   if (primal_face_is_neumann[face_id]) {
+      //     S(face_dof_number[face_id], face_dof_number[face_id]) += mes_l[face_id] * delta[face_id] ;
+      //   }
+      // }
+
+      CellValue<double> fj =
+        InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
+
+      const CellValue<const double> primal_Vj = mesh_data.Vj();
+
+      CRSMatrix A{S};
+      Vector<double> b{number_of_dof};
+
+      const auto& primal_cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
+      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 -A^{JR}b^R & -A^{LR}b^R
+        NodeValue<bool> node_tag{mesh->connectivity()};
+        node_tag.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, DirichletBoundaryCondition>) {
+                const auto& node_list  = bc.nodeList();
+                const auto& value_list = bc.valueList();
+                for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+                  const NodeId node_id = node_list[i_node];
+                  if (not node_tag[node_id]) {
+                    node_tag[node_id] = true;
+
+                    const auto& node_cells = primal_node_to_cell_matrix[node_id];
+                    for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
+                      const CellId cell_id = node_cells[i_cell];
+
+                      const auto& primal_cell_to_face = primal_cell_to_face_matrix[cell_id];
+
+                      for (size_t i_cell_face = 0; i_cell_face < primal_cell_to_face.size(); ++i_cell_face) {
+                        FaceId face_id = primal_cell_to_face_matrix[cell_id][i_cell_face];
+
+                        const bool is_face_reversed_for_cell_i = [=] {
+                          for (size_t i_face = 0; i_face < primal_cell_to_face.size(); ++i_face) {
+                            FaceId primal_face_id = primal_cell_to_face[i_face];
+                            if (primal_face_id == face_id) {
+                              return primal_face_cell_is_reversed(cell_id, i_face);
+                            }
+                          }
+                          Assert(false, "cannot find cell's face");
+                          return false;
+                        }();
+
+                        const auto& primal_face_to_node = primal_face_to_node_matrix[face_id];
+
+                        for (size_t i_face_node = 0; i_face_node < primal_face_to_node.size(); ++i_face_node) {
+                          if (primal_face_to_node[i_face_node] == node_id) {
+                            const TinyVector<Dimension> nil = [=] {
+                              if (is_face_reversed_for_cell_i) {
+                                return -primal_nlr(face_id, i_face_node);
+                              } else {
+                                return primal_nlr(face_id, i_face_node);
+                              }
+                            }();
+
+                            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);
+                                b[cell_dof_number[cell_id]] += alpha_l[face_id] * value_list[i_node] * (nil, Clr);
+                                if (primal_face_is_neumann[face_id]) {
+                                  b[face_dof_number[face_id]] += alpha_l[face_id] * value_list[i_node] * (nil, Clr);
+                                }
+                              }
+                            }
+                          }
+                        }
+                      }
+                    }
+                  }
+                }
+              }
+            },
+            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];
+                Assert(face_to_cell_matrix[face_id].size() == 1);
+                b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face];
+              }
+            }
+          },
+          boundary_condition);
+      }
+      // EL COMMENTE
+      // 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& node_list  = bc.nodeList();
+      //         const auto& value_list = bc.valueList();
+      //         for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+      //           NodeId node_id = node_list[i_node];
+      //           for (size_t i_face = 0; i_face < node_to_face_matrix[node_id].size(); ++i_face) {
+      //             FaceId face_id = node_to_face_matrix[node_id][i_face];
+      //             if (primal_face_is_neumann[face_id]) {   // HARD ajouter contribution dirchlet à face neumann
+      //             A^{LR} b^R
+      //               int nb_node_per_face = primal_face_to_node_matrix[face_id].size();
+      //               b[face_dof_number[face_id]] += (1. / nb_node_per_face) * value_list[i_node];
+      //             }
+      //           }
+      //         }
+      //       }
+      //     },
+      //     boundary_condition);
+      // }
+
+      for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+        if (primal_face_is_neumann[face_id]) {
+          CellId cell_id = face_to_cell_matrix[face_id][0];
+          const bool is_face_reversed_for_cell_i =
+            primal_face_cell_is_reversed(cell_id, face_local_numbers_in_their_cells(face_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];
+            if (is_dirichlet[node_id]) {
+              const TinyVector<Dimension> nil = [&] {
+                if (is_face_reversed_for_cell_i) {
+                  return -primal_nlr(face_id, i_node);
+                } else {
+                  return primal_nlr(face_id, i_node);
+                }
+              }();
+              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);
+                  b[cell_dof_number[cell_id]] -= alpha_l[face_id] * (nil, Clr) * dirichlet_value[node_id];   ///
+                }
+              }
+            }
+          }
+        }
+      }
+
+      Vector<double> T{number_of_dof};
+      T = 0;
+      std::cout << " Matrix " << S << "\n";
+      for (size_t k = 0; k < b.size(); ++k) {
+        std::cout << " b[" << k << "]=" << b[k] << ", ";
+      }
+      std::cout << "\n";
+      BiCGStab{b, A, T, 1000, 1e-15};
+
+      Vector r = A * T - b;
+
+      std::cout << "real residu = " << std::sqrt((r, r)) << '\n';
+
+      CellValue<double> Temperature{mesh->connectivity()};
+
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_dof_number[cell_id]]; });
+
+      Vector<double> error{mesh->numberOfCells()};
+      CellValue<double> cell_error{mesh->connectivity()};
+      double error_max = 0.;
+      size_t cell_max  = 0;
+      parallel_for(
+        mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+          error[cell_id]      = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]);
+          cell_error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]);
+        });
+      for (size_t cell_id = 0; cell_id < mesh->numberOfCells(); cell_id++) {
+        if (error_max < std::abs(error[cell_id])) {
+          error_max = std::abs(error[cell_id]);
+          cell_max  = cell_id;
+        }
+      }
+
+      std::cout << " ||Error||_max = " << error_max << " on cell " << cell_max << "\n";
+      std::cout << "||Error||_2 = " << std::sqrt((error, error)) << "\n";
+
+      {
+        VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01);
+
+        vtk_writer.write(mesh,
+                         {NamedItemValue{"T", Temperature}, NamedItemValue{"Vj", primal_Vj},
+                          NamedItemValue{"Texact", Tj}, NamedItemValue{"error", cell_error}},
+                         0,
+                         true);   // forces last output
+      }
+    }
+  } else {
+    throw NotImplementedError("not done in 1d");
+  }
+}
+
+template <size_t Dimension>
+class HeatDiamondScheme2<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}
+  {
+    Assert(m_value_list.size() == m_node_list.size());
+  }
+
+  ~DirichletBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class HeatDiamondScheme2<Dimension>::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;
+};
+
+template <size_t Dimension>
+class HeatDiamondScheme2<Dimension>::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;
+};
+
+template <size_t Dimension>
+class HeatDiamondScheme2<Dimension>::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;
+};
+
+template HeatDiamondScheme2<1>::HeatDiamondScheme2(
+  std::shared_ptr<const IMesh>,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&);
+
+template HeatDiamondScheme2<2>::HeatDiamondScheme2(
+  std::shared_ptr<const IMesh>,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&);
+
+template HeatDiamondScheme2<3>::HeatDiamondScheme2(
+  std::shared_ptr<const IMesh>,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&);

src/language/algorithms/HeatDiamondAlgorithm2.hpp

+#ifndef HEAT_DIAMOND_ALGORITHM2_HPP
+#define HEAT_DIAMOND_ALGORITHM2_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+#include <variant>
+#include <vector>
+
+template <size_t Dimension>
+class HeatDiamondScheme2
+{
+ private:
+  class DirichletBoundaryCondition;
+  class FourierBoundaryCondition;
+  class NeumannBoundaryCondition;
+  class SymmetryBoundaryCondition;
+
+ public:
+  HeatDiamondScheme2(std::shared_ptr<const IMesh> i_mesh,
+                     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+                     const FunctionSymbolId& T_id,
+                     const FunctionSymbolId& kappa_id,
+                     const FunctionSymbolId& f_id);
+};
+
+#endif   // HEAT_DIAMOND_ALGORITHM2_HPP