#include <scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.hpp>

#include <language/utils/InterpolateItemArray.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshEdgeBoundary.hpp>
#include <mesh/MeshFaceBoundary.hpp>
#include <mesh/MeshFlatEdgeBoundary.hpp>
#include <mesh/MeshFlatFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshTraits.hpp>
#include <mesh/MeshVariant.hpp>
#include <mesh/SubItemValuePerItemUtils.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
#include <variant>

template <MeshConcept MeshTypeT>
class HybridHLLcRusanovEulerianCompositeSolver_v2
{
 private:
  using MeshType = MeshTypeT;

  static constexpr size_t Dimension = MeshType::Dimension;

  using Rdxd = TinyMatrix<Dimension>;
  using Rd   = TinyVector<Dimension>;

  using Rpxp = TinyMatrix<Dimension + 2>;
  using Rp   = TinyVector<Dimension + 2>;

  class SymmetryBoundaryCondition;
  class InflowListBoundaryCondition;
  class OutflowBoundaryCondition;
  class WallBoundaryCondition;
  class NeumannflatBoundaryCondition;

  using BoundaryCondition = std::variant<SymmetryBoundaryCondition,
                                         InflowListBoundaryCondition,
                                         OutflowBoundaryCondition,
                                         NeumannflatBoundaryCondition,
                                         WallBoundaryCondition>;

  using BoundaryConditionList = std::vector<BoundaryCondition>;

  BoundaryConditionList
  _getBCList(const MeshType& mesh,
             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
  {
    BoundaryConditionList bc_list;

    for (const auto& bc_descriptor : bc_descriptor_list) {
      bool is_valid_boundary_condition = true;

      switch (bc_descriptor->type()) {
      case IBoundaryConditionDescriptor::Type::wall: {
        if constexpr (Dimension == 2) {
          bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                                     getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
        } else {
          static_assert(Dimension == 3);
          bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                                     getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                                     getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
        }
        break;
      }
      case IBoundaryConditionDescriptor::Type::symmetry: {
        if constexpr (Dimension == 2) {
          bc_list.emplace_back(
            SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                      getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
        } else {
          static_assert(Dimension == 3);
          bc_list.emplace_back(
            SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                      getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                      getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
        }
        break;
      }
      case IBoundaryConditionDescriptor::Type::inflow_list: {
        const InflowListBoundaryConditionDescriptor& inflow_list_bc_descriptor =
          dynamic_cast<const InflowListBoundaryConditionDescriptor&>(*bc_descriptor);
        if (inflow_list_bc_descriptor.functionSymbolIdList().size() != 2 + Dimension) {
          std::ostringstream error_msg;
          error_msg << "invalid number of functions for inflow boundary "
                    << inflow_list_bc_descriptor.boundaryDescriptor() << ", found "
                    << inflow_list_bc_descriptor.functionSymbolIdList().size() << ", expecting " << 2 + Dimension;
          throw NormalError(error_msg.str());
        }

        if constexpr (Dimension == 2) {
          auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
          Table<const double> node_values =
            InterpolateItemArray<double(Rd)>::template interpolate<ItemType::node>(inflow_list_bc_descriptor
                                                                                     .functionSymbolIdList(),
                                                                                   mesh.xr(), node_boundary.nodeList());

          auto xl = MeshDataManager::instance().getMeshData(mesh).xl();

          auto face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
          Table<const double> face_values =
            InterpolateItemArray<double(Rd)>::template interpolate<ItemType::face>(inflow_list_bc_descriptor
                                                                                     .functionSymbolIdList(),
                                                                                   xl, face_boundary.faceList());

          bc_list.emplace_back(InflowListBoundaryCondition(node_boundary, face_boundary, node_values, face_values));
        } else {
          static_assert(Dimension == 3);
          auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
          Table<const double> node_values =
            InterpolateItemArray<double(Rd)>::template interpolate<ItemType::node>(inflow_list_bc_descriptor
                                                                                     .functionSymbolIdList(),
                                                                                   mesh.xr(), node_boundary.nodeList());

          auto xe = MeshDataManager::instance().getMeshData(mesh).xe();

          auto edge_boundary = getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor());
          Table<const double> edge_values =
            InterpolateItemArray<double(Rd)>::template interpolate<ItemType::edge>(inflow_list_bc_descriptor
                                                                                     .functionSymbolIdList(),
                                                                                   xe, edge_boundary.edgeList());

          auto xl = MeshDataManager::instance().getMeshData(mesh).xl();

          auto face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
          Table<const double> face_values =
            InterpolateItemArray<double(Rd)>::template interpolate<ItemType::face>(inflow_list_bc_descriptor
                                                                                     .functionSymbolIdList(),
                                                                                   xl, face_boundary.faceList());

          bc_list.emplace_back(InflowListBoundaryCondition(node_boundary, edge_boundary, face_boundary, node_values,
                                                           edge_values, face_values));
        }
        break;
      }
      case IBoundaryConditionDescriptor::Type::outflow: {
        if constexpr (Dimension == 2) {
          bc_list.emplace_back(
            OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                     getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
        } else {
          static_assert(Dimension == 3);
          bc_list.emplace_back(
            OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                     getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                     getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
        }
        break;
        // std::cout << "outflow not implemented yet\n";
        // 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 Rusanov v2 Eulerian Composite solver";
        throw NormalError(error_msg.str());
      }
    }

    return bc_list;
  }

 public:
  void
  _applyWallBoundaryCondition(const BoundaryConditionList& bc_list,
                              const MeshType& mesh,
                              NodeValuePerCell<Rp>& stateNode,
                              EdgeValuePerCell<Rp>& stateEdge,
                              FaceValuePerCell<Rp>& stateFace) const
  {
    for (const auto& boundary_condition : bc_list) {
      std::visit(
        [&](auto&& bc) {
          using T = std::decay_t<decltype(bc)>;
          if constexpr (std::is_same_v<WallBoundaryCondition, T>) {
            MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
            std::cout << " Traitement WALL local (non flat) \n";
            // const Rd& normal = bc.outgoingNormal();

            const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
            const auto& node_to_face_matrix = mesh.connectivity().nodeToFaceMatrix();
            const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();

            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
            const auto& node_local_numbers_in_their_faces = mesh.connectivity().nodeLocalNumbersInTheirFaces();
            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
            // const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();

            const auto& face_list = bc.faceList();
            const auto& node_list = bc.nodeList();

            // const auto Cjr = mesh_data.Cjr();
            const auto Cjf = mesh_data.Cjf();

            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
              const NodeId node_id = node_list[i_node];

              const auto& node_face_list = node_to_face_matrix[node_id];
              // Assert(face_cell_list.size() == 1);
              // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);

              // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
              Rd normal(zero);
              int nbnormal = 0;
              for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
                FaceId node_face_id = node_face_list[i_face];

                for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) {
                  const FaceId facebc_id = face_list[i_facebc];
                  if (node_face_id == facebc_id) {
                    const auto& face_cell_list = face_to_cell_matrix[facebc_id];
                    Assert(face_cell_list.size() == 1);

                    CellId face_cell_id              = face_cell_list[0];
                    size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);

                    // Normal locale approchée
                    Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell);
                    normalloc *= 1. / l2Norm(normalloc);
                    normal += normalloc;
                    ++nbnormal;
                    break;
                  }
                }
              }
              if (nbnormal == 0)
                continue;
              normal *= 1. / nbnormal;

              normal *= 1. / l2Norm(normal);
              const auto& node_cell_list = node_to_cell_matrix[node_id];
              // Assert(face_cell_list.size() == 1);
              const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);

              for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
                CellId node_cell_id              = node_cell_list[i_cell];
                size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];

                Rd vectorSym(zero);
                for (size_t dim = 0; dim < Dimension; ++dim)
                  vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];

                vectorSym -= dot(vectorSym, normal) * normal;

                for (size_t dim = 0; dim < Dimension; ++dim)
                  stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
                // stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = 0;
              }
            }

            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
              const FaceId face_id = face_list[i_face];

              const auto& face_cell_list = face_to_cell_matrix[face_id];
              Assert(face_cell_list.size() == 1);

              CellId face_cell_id              = face_cell_list[0];
              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);

              // Normal locale approchée
              Rd normal(Cjf(face_cell_id, face_local_number_in_cell));
              normal *= 1. / l2Norm(normal);

              Rd vectorSym(zero);
              for (size_t dim = 0; dim < Dimension; ++dim)
                vectorSym[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];

              vectorSym -= dot(vectorSym, normal) * normal;

              for (size_t dim = 0; dim < Dimension; ++dim)
                stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim];
              // stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = 0;
            }

            if constexpr (Dimension == 3) {
              const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();

              const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();

              const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix();

              const auto& edge_local_numbers_in_their_faces = mesh.connectivity().edgeLocalNumbersInTheirFaces();

              const auto& edge_list = bc.edgeList();

              for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
                const EdgeId edge_id       = edge_list[i_edge];
                const auto& edge_face_list = edge_to_face_matrix[edge_id];

                // on va chercher les normale d'une face issue du edge de CL et contenue dans le faceList
                Rd normal(zero);
                int nbnormal = 0;
                for (size_t i_face = 0; i_face < edge_face_list.size(); ++i_face) {
                  FaceId edge_face_id = edge_face_list[i_face];

                  for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) {
                    const FaceId facebc_id = face_list[i_facebc];
                    if (edge_face_id == facebc_id) {
                      const auto& face_cell_list = face_to_cell_matrix[facebc_id];
                      Assert(face_cell_list.size() == 1);

                      CellId face_cell_id              = face_cell_list[0];
                      size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);

                      // Normal locale approchée
                      Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell);
                      normalloc *= 1. / l2Norm(normalloc);
                      normal += normalloc;
                      ++nbnormal;
                      break;
                    }
                  }
                }

                if (nbnormal == 0)
                  continue;
                normal *= 1. / nbnormal;

                normal *= 1. / l2Norm(normal);

                const auto& edge_cell_list = edge_to_cell_matrix[edge_id];

                const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);

                for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
                  CellId edge_cell_id              = edge_cell_list[i_cell];
                  size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];

                  Rd vectorSym(zero);
                  for (size_t dim = 0; dim < Dimension; ++dim)
                    vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim];

                  vectorSym -= dot(vectorSym, normal) * normal;

                  for (size_t dim = 0; dim < Dimension; ++dim)
                    stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
                  // stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = 0;
                }
              }
            }
          }
        },
        boundary_condition);
    }
  }

  void
  _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
                                 const MeshType& mesh,
                                 NodeValuePerCell<Rp>& stateNode,
                                 EdgeValuePerCell<Rp>& stateEdge,
                                 FaceValuePerCell<Rp>& stateFace) const
  {
    for (const auto& boundary_condition : bc_list) {
      std::visit(
        [&](auto&& bc) {
          using T = std::decay_t<decltype(bc)>;
          if constexpr (std::is_same_v<OutflowBoundaryCondition, T>) {
            std::cout << " Traitement Outflow  \n";
            // const Rd& normal = bc.outgoingNormal();
            /*
            const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
            const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();

            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
            // const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();

            const auto& face_list = bc.faceList();
            const auto& node_list = bc.nodeList();

            const auto xj = mesh.xj();
            const auto xr = mesh.xr();
            const auto xf = mesh.xl();
            const auto xe = mesh.xe();

            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
              const NodeId node_id = node_list[i_node];

              const auto& node_cell_list = node_to_cell_matrix[node_id];
              // Assert(face_cell_list.size() == 1);
              const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);

              for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
                CellId node_cell_id              = node_cell_list[i_cell];
                size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];

                for (size_t dim = 0; dim < Dimension + 2; ++dim)
                  stateNode[node_cell_id][node_local_number_in_cell][dim] += vectorSym[dim];

                Rd vectorSym(zero);
                for (size_t dim = 0; dim < Dimension; ++dim)
                  vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];

                Rdxd MatriceProj(identity);
                MatriceProj -= tensorProduct(normal, normal);
                vectorSym = MatriceProj * vectorSym;

                for (size_t dim = 0; dim < Dimension; ++dim)
                  stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
                //  stateNode[node_cell_id][node_local_number_in_cell][dim] = 0;   // node_array_list[i_node][dim];
              }
            }

            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
              const FaceId face_id = face_list[i_face];

              const auto& face_cell_list = face_to_cell_matrix[face_id];
              Assert(face_cell_list.size() == 1);

              CellId face_cell_id              = face_cell_list[0];
              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);

              Rd vectorSym(zero);
              for (size_t dim = 0; dim < Dimension; ++dim)
                vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];

              Rdxd MatriceProj(identity);
              MatriceProj -= tensorProduct(normal, normal);
              vectorSym = MatriceProj * vectorSym;

              for (size_t dim = 0; dim < Dimension; ++dim)
                stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim];
            }

            if constexpr (Dimension == 3) {
              const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();

              const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();

              const auto& edge_list = bc.edgeList();

              for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
                const EdgeId edge_id = edge_list[i_edge];

                const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
                // Assert(face_cell_list.size() == 1);
                const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);

                for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
                  CellId edge_cell_id              = edge_cell_list[i_cell];
                  size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];

                  Rd vectorSym(zero);
                  for (size_t dim = 0; dim < Dimension; ++dim)
                    vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim];

                  Rdxd MatriceProj(identity);
                  MatriceProj -= tensorProduct(normal, normal);
                  vectorSym = MatriceProj * vectorSym;

                  for (size_t dim = 0; dim < Dimension; ++dim)
                    stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
                }
              }

              //          throw NormalError("Not implemented");
            }
            */
          }
        },
        boundary_condition);
    }
  }

  void
  _applySymmetricBoundaryCondition(const BoundaryConditionList& bc_list,
                                   const MeshType& mesh,
                                   NodeValuePerCell<Rp>& stateNode,
                                   EdgeValuePerCell<Rp>& stateEdge,
                                   FaceValuePerCell<Rp>& stateFace) const
  {
    for (const auto& boundary_condition : bc_list) {
      std::visit(
        [&](auto&& bc) {
          using T = std::decay_t<decltype(bc)>;
          if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
            // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
            std::cout << " Traitement SYMMETRY  \n";
            const Rd& normal = bc.outgoingNormal();

            const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
            const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();

            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
            // const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();

            const auto& face_list = bc.faceList();
            const auto& node_list = bc.nodeList();

            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
              const NodeId node_id = node_list[i_node];

              const auto& node_cell_list = node_to_cell_matrix[node_id];
              // Assert(face_cell_list.size() == 1);
              const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);

              for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
                CellId node_cell_id              = node_cell_list[i_cell];
                size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];

                Rd vectorSym(zero);
                for (size_t dim = 0; dim < Dimension; ++dim)
                  vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];

                Rdxd MatriceProj(identity);
                MatriceProj -= tensorProduct(normal, normal);
                vectorSym = MatriceProj * vectorSym;

                for (size_t dim = 0; dim < Dimension; ++dim)
                  stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
                //  stateNode[node_cell_id][node_local_number_in_cell][dim] = 0;   // node_array_list[i_node][dim];
              }
            }

            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
              const FaceId face_id = face_list[i_face];

              const auto& face_cell_list = face_to_cell_matrix[face_id];
              Assert(face_cell_list.size() == 1);

              CellId face_cell_id              = face_cell_list[0];
              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);

              Rd vectorSym(zero);
              for (size_t dim = 0; dim < Dimension; ++dim)
                vectorSym[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];

              Rdxd MatriceProj(identity);
              MatriceProj -= tensorProduct(normal, normal);
              vectorSym = MatriceProj * vectorSym;

              for (size_t dim = 0; dim < Dimension; ++dim)
                stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim];
            }

            if constexpr (Dimension == 3) {
              const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();

              const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();

              const auto& edge_list = bc.edgeList();

              for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
                const EdgeId edge_id = edge_list[i_edge];

                const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
                // Assert(face_cell_list.size() == 1);
                const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);

                for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
                  CellId edge_cell_id              = edge_cell_list[i_cell];
                  size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];

                  Rd vectorSym(zero);
                  for (size_t dim = 0; dim < Dimension; ++dim)
                    vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim];

                  Rdxd MatriceProj(identity);
                  MatriceProj -= tensorProduct(normal, normal);
                  vectorSym = MatriceProj * vectorSym;

                  for (size_t dim = 0; dim < Dimension; ++dim)
                    stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
                }
              }
            }
          }
        },
        boundary_condition);
    }
  }

  void
  _applyNeumannflatBoundaryCondition(const BoundaryConditionList& bc_list,
                                     const MeshType& mesh,
                                     NodeValuePerCell<Rp>& stateNode,
                                     EdgeValuePerCell<Rp>& stateEdge,
                                     FaceValuePerCell<Rp>& stateFace) const
  {
    for (const auto& boundary_condition : bc_list) {
      std::visit(
        [&](auto&& bc) {
          using T = std::decay_t<decltype(bc)>;
          if constexpr (std::is_same_v<NeumannflatBoundaryCondition, T>) {
            // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
            std::cout << " Traitement WALL  \n";
            const Rd& normal = bc.outgoingNormal();

            const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
            const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();

            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
            // const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();

            const auto& face_list = bc.faceList();
            const auto& node_list = bc.nodeList();

            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
              const NodeId node_id = node_list[i_node];

              const auto& node_cell_list = node_to_cell_matrix[node_id];
              // Assert(face_cell_list.size() == 1);
              const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);

              for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
                CellId node_cell_id              = node_cell_list[i_cell];
                size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];

                Rd vectorSym(zero);
                for (size_t dim = 0; dim < Dimension; ++dim)
                  vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];

                vectorSym -= dot(vectorSym, normal) * normal;

                for (size_t dim = 0; dim < Dimension; ++dim)
                  stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
                //  stateNode[node_cell_id][node_local_number_in_cell][dim] = 0;   // node_array_list[i_node][dim];
              }
            }

            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
              const FaceId face_id = face_list[i_face];

              const auto& face_cell_list = face_to_cell_matrix[face_id];
              Assert(face_cell_list.size() == 1);

              CellId face_cell_id              = face_cell_list[0];
              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);

              Rd vectorSym(zero);
              for (size_t dim = 0; dim < Dimension; ++dim)
                vectorSym[dim] = stateFace[face_cell_id][face_local_number_in_cell][1 + dim];

              vectorSym -= dot(vectorSym, normal) * normal;

              for (size_t dim = 0; dim < Dimension; ++dim)
                stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim];
            }

            if constexpr (Dimension == 3) {
              const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();

              const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();

              const auto& edge_list = bc.edgeList();

              for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
                const EdgeId edge_id = edge_list[i_edge];

                const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
                // Assert(face_cell_list.size() == 1);
                const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);

                for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
                  CellId edge_cell_id              = edge_cell_list[i_cell];
                  size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];

                  Rd vectorSym(zero);
                  for (size_t dim = 0; dim < Dimension; ++dim)
                    vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim];

                  vectorSym -= dot(vectorSym, normal) * normal;

                  for (size_t dim = 0; dim < Dimension; ++dim)
                    stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
                }
              }
            }
          }
        },
        boundary_condition);
    }
  }

  void
  _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list,
                                const MeshType& mesh,
                                NodeValuePerCell<Rp>& stateNode,
                                EdgeValuePerCell<Rp>& stateEdge,
                                FaceValuePerCell<Rp>& stateFace) const
  {
    for (const auto& boundary_condition : bc_list) {
      std::visit(
        [&](auto&& bc) {
          using T = std::decay_t<decltype(bc)>;
          if constexpr (std::is_same_v<InflowListBoundaryCondition, T>) {
            // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
            std::cout << " Traitement INFLOW  \n";

            const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
            const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();

            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
            // const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();

            const auto& face_list = bc.faceList();
            const auto& node_list = bc.nodeList();

            const auto& face_array_list = bc.faceArrayList();
            const auto& node_array_list = bc.nodeArrayList();

            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
              const NodeId node_id = node_list[i_node];

              const auto& node_cell_list = node_to_cell_matrix[node_id];
              // Assert(face_cell_list.size() == 1);
              const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);

              for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
                CellId node_cell_id              = node_cell_list[i_cell];
                size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];

                for (size_t dim = 0; dim < Dimension + 2; ++dim)
                  stateNode[node_cell_id][node_local_number_in_cell][dim] = node_array_list[i_node][dim];
              }
            }

            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
              const FaceId face_id = face_list[i_face];

              const auto& face_cell_list = face_to_cell_matrix[face_id];
              Assert(face_cell_list.size() == 1);

              CellId face_cell_id              = face_cell_list[0];
              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);

              for (size_t dim = 0; dim < Dimension + 2; ++dim)
                stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][dim];
            }

            if constexpr (Dimension == 3) {
              const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();

              const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();
              // const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();

              const auto& edge_list = bc.edgeList();

              const auto& edge_array_list = bc.edgeArrayList();

              for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
                const EdgeId edge_id = edge_list[i_edge];

                const auto& edge_cell_list                 = edge_to_cell_matrix[edge_id];
                const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);

                // Assert(face_cell_list.size() == 1);

                for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
                  CellId edge_cell_id              = edge_cell_list[i_cell];
                  size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];

                  for (size_t dim = 0; dim < Dimension + 2; ++dim)
                    stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim];
                }
              }
            }
          }
        },
        boundary_condition);
    }
  }

 public:
  inline double
  pression(const double rho, const double epsilon, const double gam) const
  {
    return (gam - 1) * rho * epsilon;
  }

  std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
             std::shared_ptr<const DiscreteFunctionVariant>,
             std::shared_ptr<const DiscreteFunctionVariant>>
  solve(const std::shared_ptr<const MeshType>& p_mesh,
        const DiscreteFunctionP0<const double>& rho_n,
        const DiscreteFunctionP0<const Rd>& u_n,
        const DiscreteFunctionP0<const double>& E_n,
        const double& gamma,
        const DiscreteFunctionP0<const double>& c_n,
        const DiscreteFunctionP0<const double>& p_n,
        // const size_t degree,
        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
        const double& dt,
        const bool checkLocalConservation) const
  {
    auto rho = copy(rho_n);
    auto u   = copy(u_n);
    auto E   = copy(E_n);
    // auto c   = copy(c_n);
    // auto p   = copy(p_n);

    auto bc_list = this->_getBCList(*p_mesh, bc_descriptor_list);

    auto rhoE = rho * E;
    auto rhoU = rho * u;

    // Construction du vecteur des variables conservatives
    //
    // Ci dessous juste ordre 1
    //
    CellValue<Rp> State{p_mesh->connectivity()};
    parallel_for(
      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
        State[j][0] = rho[j];
        for (size_t d = 0; d < Dimension; ++d)
          State[j][1 + d] = rhoU[j][d];
        State[j][1 + Dimension] = rhoE[j];
      });

    // CellValue<Rp> State{p_mesh->connectivity()};
    //
    // Remplir ici les reconstructions d'ordre élevé

    //
    const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
    const auto& cell_to_edge_matrix = p_mesh->connectivity().cellToEdgeMatrix();
    const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();

    //    const auto xr = p_mesh->xr();
    // auto xl       = MeshDataManager::instance().getMeshData(*p_mesh).xl();
    // auto xe       = MeshDataManager::instance().getMeshData(*p_mesh).xe();

    NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()};
    StateAtNode.fill(zero);
    parallel_for(
      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtNode[j].fill(State[j]); });

    EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()};
    StateAtEdge.fill(zero);
    parallel_for(
      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtEdge[j].fill(State[j]); });
    FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()};
    StateAtFace.fill(zero);
    parallel_for(
      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtFace[j].fill(State[j]); });

    // Conditions aux limites des etats conservatifs

    _applyInflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
    //_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
    _applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
    _applyNeumannflatBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
    _applyWallBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);

    //
    // Construction du flux .. ok pour l'ordre 1
    //
    CellValue<Rdxd> rhoUtensU{p_mesh->connectivity()};
    CellValue<Rdxd> Pid(p_mesh->connectivity());
    Pid.fill(identity);
    CellValue<Rdxd> rhoUtensUPlusPid{p_mesh->connectivity()};
    rhoUtensUPlusPid.fill(zero);

    parallel_for(
      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
        rhoUtensU[j] = tensorProduct(rhoU[j], u[j]);
        Pid[j] *= p_n[j];
        rhoUtensUPlusPid[j] = rhoUtensU[j] + Pid[j];
      });

    auto Flux_rho    = rhoU;
    auto Flux_qtmvt  = rhoUtensUPlusPid;   // rhoUtensU + Pid;
    auto Flux_totnrj = (rhoE + p_n) * u;

    // Flux avec prise en compte des states at Node/Edge/Face

    NodeValuePerCell<Rd> Flux_rhoAtCellNode{p_mesh->connectivity()};
    EdgeValuePerCell<Rd> Flux_rhoAtCellEdge{p_mesh->connectivity()};
    FaceValuePerCell<Rd> Flux_rhoAtCellFace{p_mesh->connectivity()};
    /*
    parallel_for(
      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
        const auto& cell_to_node = cell_to_node_matrix[j];

        for (size_t l = 0; l < cell_to_node.size(); ++l) {
          for (size_t dim = 0; dim < Dimension; ++dim)
            Flux_rhoAtCellNode[j][l][dim] = StateAtNode[j][l][0] * StateAtNode[j][l][dim];
        }

        const auto& cell_to_face = cell_to_face_matrix[j];

        for (size_t l = 0; l < cell_to_face.size(); ++l) {
          for (size_t dim = 0; dim < Dimension; ++dim)
            Flux_rhoAtCellFace[j][l][dim] = StateAtFace[j][l][0] * StateAtFace[j][l][dim];
        }

        const auto& cell_to_edge = cell_to_edge_matrix[j];

        for (size_t l = 0; l < cell_to_edge.size(); ++l) {
          for (size_t dim = 0; dim < Dimension; ++dim)
            Flux_rhoAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim];
        }
      });
  */
    NodeValuePerCell<Rdxd> Flux_qtmvtAtCellNode{p_mesh->connectivity()};   // = rhoUtensUPlusPid;   // rhoUtensU + Pid;
    EdgeValuePerCell<Rdxd> Flux_qtmvtAtCellEdge{p_mesh->connectivity()};   // = rhoUtensUPlusPid;   // rhoUtensU + Pid;
    FaceValuePerCell<Rdxd> Flux_qtmvtAtCellFace{p_mesh->connectivity()};   // = rhoUtensUPlusPid;   // rhoUtensU + Pid;
                                                                           /*
                                                                           parallel_for(
                                                                             p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
                                                                               const auto& cell_to_node = cell_to_node_matrix[j];
                                                                       for (size_t l = 0; l < cell_to_node.size(); ++l) {
                                                                                   for (size_t dim = 0; dim < Dimension; ++dim)
                                                                                     Flux_qtmvtAtCellNode[j][l][dim] = StateAtNode[j][l][0] * StateAtNode[j][l][dim];
                                                                                 }

                                                                                 const auto& cell_to_face = cell_to_face_matrix[j];

                                                                                 for (size_t l = 0; l < cell_to_face.size(); ++l) {
                                                                                   for (size_t dim = 0; dim < Dimension; ++dim)
                                                                                     Flux_qtmvtAtCellFace[j][l][dim] = StateAtFace[j][l][0] * StateAtFace[j][l][dim];
                                                                                 }

                                                                                 const auto& cell_to_edge = cell_to_edge_matrix[j];

                                                                                 for (size_t l = 0; l < cell_to_edge.size(); ++l) {
                                                                                   for (size_t dim = 0; dim < Dimension; ++dim)
                                                                                     Flux_qtmvtAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim];
                                                                                 }
                                                                               });
                                                                           */
    // parallel_for(
    //   p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
    //     Flux_qtmvtAtCellNode[j] = Flux_qtmvtAtCellEdge[j] = Flux_qtmvtAtCellFace[j] = Flux_qtmvt[j];
    //   });

    NodeValuePerCell<Rd> Flux_totnrjAtCellNode{p_mesh->connectivity()};
    EdgeValuePerCell<Rd> Flux_totnrjAtCellEdge{p_mesh->connectivity()};
    FaceValuePerCell<Rd> Flux_totnrjAtCellFace{p_mesh->connectivity()};

    Flux_rhoAtCellEdge.fill(zero);
    Flux_totnrjAtCellEdge.fill(zero);
    Flux_qtmvtAtCellEdge.fill(zero);

    // Les flux aux nodes/edge/faces
    parallel_for(
      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
        const auto& cell_to_node = cell_to_node_matrix[j];

        for (size_t l = 0; l < cell_to_node.size(); ++l) {
          // Etats conservatifs aux noeuds
          const double rhonode = StateAtNode[j][l][0];
          Rd Unode;
          for (size_t dim = 0; dim < Dimension; ++dim)
            Unode[dim] = StateAtNode[j][l][dim + 1] / rhonode;
          const double rhoEnode = StateAtNode[j][l][Dimension + 1];
          //
          const double Enode       = rhoEnode / rhonode;
          const double epsilonnode = Enode - .5 * dot(Unode, Unode);
          const Rd rhoUnode        = rhonode * Unode;
          const Rdxd rhoUtensUnode = tensorProduct(rhoUnode, Unode);

          const double Pressionnode = pression(rhonode, epsilonnode, gamma);

          const double rhoEnodePlusP = rhoEnode + Pressionnode;

          Rdxd rhoUtensUPlusPidnode(identity);
          rhoUtensUPlusPidnode *= Pressionnode;
          rhoUtensUPlusPidnode += rhoUtensUnode;

          Flux_rhoAtCellNode[j][l]    = rhoUnode;
          Flux_qtmvtAtCellNode[j][l]  = rhoUtensUPlusPidnode;
          Flux_totnrjAtCellNode[j][l] = rhoEnodePlusP * Unode;
        }

        const auto& cell_to_face = cell_to_face_matrix[j];

        for (size_t l = 0; l < cell_to_face.size(); ++l) {
          const double rhoface = StateAtFace[j][l][0];
          Rd Uface;
          for (size_t dim = 0; dim < Dimension; ++dim)
            Uface[dim] = StateAtFace[j][l][dim + 1] / rhoface;
          const double rhoEface = StateAtFace[j][l][Dimension + 1];
          //
          const double Eface       = rhoEface / rhoface;
          const double epsilonface = Eface - .5 * dot(Uface, Uface);
          const Rd rhoUface        = rhoface * Uface;
          const Rdxd rhoUtensUface = tensorProduct(rhoUface, Uface);

          const double Pressionface = pression(rhoface, epsilonface, gamma);

          const double rhoEfacePlusP = rhoEface + Pressionface;

          Rdxd rhoUtensUPlusPidface(identity);
          rhoUtensUPlusPidface *= Pressionface;
          rhoUtensUPlusPidface += rhoUtensUface;

          Flux_rhoAtCellFace[j][l]    = rhoUface;
          Flux_qtmvtAtCellFace[j][l]  = rhoUtensUPlusPidface;
          Flux_totnrjAtCellFace[j][l] = rhoEfacePlusP * Uface;
        }

        if constexpr (Dimension == 3) {
          const auto& cell_to_edge = cell_to_edge_matrix[j];

          for (size_t l = 0; l < cell_to_edge.size(); ++l) {
            const double rhoedge = StateAtEdge[j][l][0];
            Rd Uedge;
            for (size_t dim = 0; dim < Dimension; ++dim)
              Uedge[dim] = StateAtEdge[j][l][dim + 1] / rhoedge;
            const double rhoEedge = StateAtEdge[j][l][Dimension + 1];
            //
            const double Eedge       = rhoEedge / rhoedge;
            const double epsilonedge = Eedge - .5 * dot(Uedge, Uedge);
            const Rd rhoUedge        = rhoedge * Uedge;
            const Rdxd rhoUtensUedge = tensorProduct(rhoUedge, Uedge);

            const double Pressionedge = pression(rhoedge, epsilonedge, gamma);

            const double rhoEedgePlusP = rhoEedge + Pressionedge;

            Rdxd rhoUtensUPlusPidedge(identity);
            rhoUtensUPlusPidedge *= Pressionedge;
            rhoUtensUPlusPidedge += rhoUtensUedge;

            Flux_rhoAtCellEdge[j][l]    = rhoUedge;
            Flux_qtmvtAtCellEdge[j][l]  = rhoUtensUPlusPidedge;
            Flux_totnrjAtCellEdge[j][l] = rhoEedgePlusP * Uedge;
          }
        }
      });

    // parallel_for(
    //   p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
    //     Flux_totnrjAtCellNode[j] = Flux_totnrjAtCellEdge[j] = Flux_totnrjAtCellFace[j] = Flux_totnrj[j];
    //   });

    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);

    auto volumes_cell = mesh_data.Vj();

    // Calcul des matrices de viscosité aux node/edge/face

    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
    const NodeValuePerCell<const Rd> njr = mesh_data.njr();

    const FaceValuePerCell<const Rd> Cjf = mesh_data.Cjf();
    const FaceValuePerCell<const Rd> njf = mesh_data.njf();

    const auto& node_to_cell_matrix               = p_mesh->connectivity().nodeToCellMatrix();
    const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();

    const auto& face_to_cell_matrix               = p_mesh->connectivity().faceToCellMatrix();
    const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();

    // We compute Gjr, Gjf

    NodeValuePerCell<Rp> Gjr{p_mesh->connectivity()};
    Gjr.fill(zero);
    EdgeValuePerCell<Rp> Gje{p_mesh->connectivity()};
    Gje.fill(zero);
    FaceValuePerCell<Rp> Gjf{p_mesh->connectivity()};
    Gjf.fill(zero);

    parallel_for(
      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
        const auto& cell_to_node = cell_to_node_matrix[j];

        for (size_t l = 0; l < cell_to_node.size(); ++l) {
          const NodeId& node                         = cell_to_node[l];
          const auto& node_to_cell                   = node_to_cell_matrix[node];
          const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node);

          const Rd& Cjr_loc = Cjr(j, l);

          for (size_t k = 0; k < node_to_cell.size(); ++k) {
            const CellId K    = node_to_cell[k];
            const size_t R    = node_local_number_in_its_cells[k];
            const Rd& Ckr_loc = Cjr(K, R);

            // Une moyenne entre les etats jk

            Rd uNode     = .5 * (u_n[j] + u_n[K]);
            double cNode = .5 * (c_n[j] + c_n[K]);

            // Viscosity j k
            Rpxp ViscosityMatrixJK(identity);
            const double MaxmaxabsVpNormjk =
              std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uNode, cNode,
                                                                                                    Cjr_loc),
                       toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uNode, cNode,
                                                                                                    Ckr_loc));

            ViscosityMatrixJK *= MaxmaxabsVpNormjk;
            const Rd& u_Cjr = Flux_qtmvtAtCellNode[K][R] * Cjr_loc;   // Flux_qtmvt[K] * Cjr_loc;

            const Rp& statediff = StateAtNode[j][l] - StateAtNode[K][R];   // State[j] - State[K];
            const Rp& diff      = ViscosityMatrixJK * statediff;

            Gjr[j][l][0] += dot(Flux_rhoAtCellNode[K][R], Cjr_loc);   // dot(Flux_rho[K], Cjr_loc);
            for (size_t d = 0; d < Dimension; ++d)
              Gjr[j][l][1 + d] += u_Cjr[d];
            Gjr[j][l][1 + Dimension] += dot(Flux_totnrjAtCellNode[K][R], Cjr_loc);   // dot(Flux_totnrj[K], Cjr_loc);

            Gjr[j][l] += diff;
          }

          Gjr[j][l] *= 1. / node_to_cell.size();
        }
      });

    synchronize(Gjr);
    if (checkLocalConservation) {
      auto is_boundary_node = p_mesh->connectivity().isBoundaryNode();

      NodeValue<double> MaxErrorConservationNode(p_mesh->connectivity());
      MaxErrorConservationNode.fill(0.);
      // double MaxErrorConservationNode = 0;
      parallel_for(
        p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId l) {
          const auto& node_to_cell                   = node_to_cell_matrix[l];
          const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(l);

          if (not is_boundary_node[l]) {
            Rp SumGjr(zero);
            double maxGjrAbs = 0;
            for (size_t k = 0; k < node_to_cell.size(); ++k) {
              const CellId K       = node_to_cell[k];
              const unsigned int R = node_local_number_in_its_cells[k];
              SumGjr += Gjr[K][R];
              maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjr[K][R]));
            }
            // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr));
            MaxErrorConservationNode[l] = l2Norm(SumGjr) / maxGjrAbs;
          }
        });
      std::cout << " Max Error Node " << max(MaxErrorConservationNode) << "\n";
    }
    //
    parallel_for(
      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
        // Edge
        const auto& cell_to_face = cell_to_face_matrix[j];

        for (size_t l = 0; l < cell_to_face.size(); ++l) {
          const FaceId& face                         = cell_to_face[l];
          const auto& face_to_cell                   = face_to_cell_matrix[face];
          const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(face);

          const Rd& Cjf_loc = Cjf(j, l);

          CellId K=face_to_cell[0];
	        unsigned int R =face_local_number_in_its_cells[0];

          if (face_to_cell.size()==1)
          {
            K=j;
            R=l;
          }
          else
          {
            const CellId K1       = face_to_cell[0];
  	        const CellId K2       = face_to_cell[1];

            if (j==K1)
            {
              K = K2;
              R = face_local_number_in_its_cells[1];
            }
          }

          const std::pair<double, double> MinMaxVpNormj = toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(u_n[j], c_n[j], Cjf_loc);
          const std::pair<double, double> MinMaxVpNormk = toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(u_n[K], c_n[K], Cjf_loc); 
          const double MinVpNormjk = std::min(MinMaxVpNormj.first, MinMaxVpNormk.first);   // SL
          const double MaxVpNormjk = std::max(MinMaxVpNormj.second, MinMaxVpNormk.second); //  SR
          Rp SC = zero;

          const Rd& uk_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc;   // Flux_qtmvt[K] * Cjf_loc;
          const Rd& uj_Cjf = Flux_qtmvtAtCellFace[j][l] * Cjf_loc;   // Flux_qtmvt[j] * Cjf_loc;
          const Rp& statediff = StateAtFace[K][R] - StateAtFace[j][l];   // State[K] - State[j];
          const Rp& diff      =  MinVpNormjk * statediff;

          SC[0] = (dot(Flux_rhoAtCellFace[K][R], Cjf_loc) - dot(Flux_rhoAtCellFace[j][l], Cjf_loc) + diff[1 + Dimension]) / diff[0];
          for (size_t d = 0; d < Dimension; ++d)
            SC[1 + d] = (uk_Cjf[d] - uj_Cjf[d] + diff[d]) / diff[d];
          SC[1 + Dimension] = (dot(Flux_totnrjAtCellFace[K][R], Cjf_loc) - dot(Flux_totnrjAtCellFace[j][l], Cjf_loc) + diff[1 + Dimension])  / diff[1 + Dimension];


          if (MinVpNormjk >= 0 ){ // SL >=0
            Gjf[j][l][0] = dot(Flux_rhoAtCellFace[j][l], Cjf_loc);   // dot(Flux_roh[j] , Cjf_loc)
            for (size_t d = 0; d < Dimension; ++d)
              Gjf[j][l][1 + d] = uj_Cjf[d];
            Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[j][l], Cjf_loc);    // dot(Flux_totnrj[K] , Cjf_loc)
          }
          else{
            if (MaxVpNormjk <= 0 ){ // SR <= 0
              Gjf[j][l][0] = dot(Flux_rhoAtCellFace[K][R], Cjf_loc);   // dot(Flux_roh[K] , Cjf_loc)
              for (size_t d = 0; d < Dimension; ++d)
                Gjf[j][l][1 + d] = uk_Cjf[d];
              Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[K][R], Cjf_loc);    // dot(Flux_totnrj[K] , Cjf_loc)
            }
            else{
              Rp diffL = zero;
              Rp diffR = zero;

              if (SC[0] >= 0 ){// SL <= 0 <= SC
                diffL[0] = MinVpNormjk * SC[0] * statediff[0];
                Gjf[j][l][0] = (SC[0] * dot(Flux_rhoAtCellFace[j][l], Cjf_loc) - MinVpNormjk * dot(Flux_rhoAtCellFace[K][R], Cjf_loc) + diffL[0]) / (SC[0] - MinVpNormjk);
              }
              else{ // SC <= 0 <= SR
                diffR[0] = SC[0] * MaxVpNormjk * statediff[0];
                Gjf[j][l][0] = (MaxVpNormjk * dot(Flux_rhoAtCellFace[j][l], Cjf_loc) - SC[0] * dot(Flux_rhoAtCellFace[K][R], Cjf_loc) + diffR[0]) / (MaxVpNormjk - SC[0]);
              }

              for (size_t d = 0; d < Dimension; ++d){
                if (SC[1 + d] >= 0 ){// SL <= 0 <= SC
                  diffL[1 + d] = MinVpNormjk * SC[1 + d] * statediff[1 + d];
                  Gjf[j][l][1 + d] = (SC[1 + d] * uj_Cjf[d] - MinVpNormjk * uk_Cjf[d] + diffL[d])  / (SC[1 + d] - MinVpNormjk);
                }
                else{ // SC <= 0 <= SR
                  diffR[1 + d] = SC[1 + d] * MaxVpNormjk * statediff[1 + d];
                  Gjf[j][l][1 + d] = (MaxVpNormjk * uj_Cjf[d] - SC[1 + d] * uk_Cjf[d] + diffR[d])  / (MaxVpNormjk - SC[1 + d]);
                }
              }

              if (SC[1 + Dimension] >= 0 ){// SL <= 0 <= SC
                diffL[1 + Dimension] = MinVpNormjk * SC[1 + Dimension] * statediff[1 + Dimension];
                Gjf[j][l][1 + Dimension] = (SC[1 + Dimension] * dot(Flux_totnrjAtCellFace[j][l], Cjf_loc) - MinVpNormjk * dot(Flux_totnrjAtCellFace[K][R], Cjf_loc) + diffL[1 + Dimension]) / (SC[1 + Dimension] - MinVpNormjk);
              }
              else{ // SC <= 0 <= SR
                diffR[1 + Dimension] = SC[1 + Dimension] * MaxVpNormjk * statediff[1 + Dimension];
                Gjf[j][l][1 + Dimension] = (MaxVpNormjk * dot(Flux_totnrjAtCellFace[j][l], Cjf_loc) - SC[1 + Dimension] * dot(Flux_totnrjAtCellFace[K][R], Cjf_loc) + diffR[1 + Dimension]) / (MaxVpNormjk - SC[1 + Dimension]); 
              }
            }
          }
        }
      });
    synchronize(Gjf);

    if (checkLocalConservation) {
      auto is_boundary_face = p_mesh->connectivity().isBoundaryFace();

      FaceValue<double> MaxErrorConservationFace(p_mesh->connectivity());
      MaxErrorConservationFace.fill(0.);

      parallel_for(
        p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
          const auto& face_to_cell                   = face_to_cell_matrix[l];
          const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(l);

          if (not is_boundary_face[l]) {
            Rp SumGjf(zero);
            double maxGjrAbs = 0;
            for (size_t k = 0; k < face_to_cell.size(); ++k) {
              const CellId K       = face_to_cell[k];
              const unsigned int R = face_local_number_in_its_cells[k];
              SumGjf += Gjf[K][R];
              maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjf[K][R]));
            }
            MaxErrorConservationFace[l] = l2Norm(SumGjf) / maxGjrAbs;
            // MaxErrorConservationFace   = std::max(MaxErrorConservationFace, l2Norm(SumGjf));
          }
        });
      std::cout << " Max Error Face " << max(MaxErrorConservationFace) << "\n";
    }

    if constexpr (Dimension == 3) {
      const auto& edge_to_cell_matrix               = p_mesh->connectivity().edgeToCellMatrix();
      const auto& edge_local_numbers_in_their_cells = p_mesh->connectivity().edgeLocalNumbersInTheirCells();

      const EdgeValuePerCell<const Rd> Cje = mesh_data.Cje();
      const EdgeValuePerCell<const Rd> nje = mesh_data.nje();

      parallel_for(
        p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
          // Edge
          const auto& cell_to_edge = cell_to_edge_matrix[j];

          for (size_t l = 0; l < cell_to_edge.size(); ++l) {
            const EdgeId& edge                         = cell_to_edge[l];
            const auto& edge_to_cell                   = edge_to_cell_matrix[edge];
            const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge);

            const Rd& Cje_loc = Cje(j, l);

            for (size_t k = 0; k < edge_to_cell.size(); ++k) {
              const CellId K = edge_to_cell[k];
              const size_t R = edge_local_number_in_its_cells[k];

              const Rd& Cke_loc = Cje(K, R);

              // Une moyenne entre les etats jk

              Rd uEdge     = .5 * (u_n[j] + u_n[K]);
              double cEdge = .5 * (c_n[j] + c_n[K]);

              // Viscosity j k
              Rpxp ViscosityMatrixJK(identity);
              const double MaxmaxabsVpNormjk =
                std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uEdge, cEdge,
                                                                                                      Cje_loc),
                         toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uEdge, cEdge,
                                                                                                      Cke_loc));

              ViscosityMatrixJK *= MaxmaxabsVpNormjk;

              const Rd& u_Cje = Flux_qtmvtAtCellEdge[K][R] * Cje_loc;   // Flux_qtmvt[K] * Cje_loc;

              const Rp& statediff = StateAtEdge[j][l] - StateAtEdge[K][R];   // State[j] - State[K];
              const Rp& diff      = ViscosityMatrixJK * statediff;

              Gje[j][l][0] += dot(Flux_rhoAtCellEdge[K][R], Cje_loc);   //  dot(Flux_rho[K], Cje_loc);
              for (size_t d = 0; d < Dimension; ++d)
                Gje[j][l][1 + d] += u_Cje[d];
              Gje[j][l][1 + Dimension] += dot(Flux_totnrjAtCellEdge[K][R], Cje_loc);   // dot(Flux_totnrj[K], Cje_loc);

              Gje[j][l] += diff;
            }

            Gje[j][l] *= 1. / edge_to_cell.size();
          }
        });
      synchronize(Gje);

      if (checkLocalConservation) {
        auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge();

        EdgeValue<double> MaxErrorConservationEdge(p_mesh->connectivity());
        MaxErrorConservationEdge.fill(0.);
        //  double MaxErrorConservationEdge = 0;
        parallel_for(
          p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId l) {
            const auto& edge_to_cell                   = edge_to_cell_matrix[l];
            const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(l);

            if (not is_boundary_edge[l]) {
              Rp SumGje(zero);
              double maxGjrAbs = 0;
              for (size_t k = 0; k < edge_to_cell.size(); ++k) {
                const CellId K       = edge_to_cell[k];
                const unsigned int R = edge_local_number_in_its_cells[k];
                SumGje += Gje[K][R];
                maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gje[K][R]));
              }
              // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje));
              MaxErrorConservationEdge[l] = l2Norm(SumGje) / maxGjrAbs;
            }
          });
        std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << "\n";
      }
    }   // dim 3

    // Pour les assemblages
    double theta = 2. / 3.;   //.5; 2. / 3.
    double eta   = 1. / 6.;   //.2; 1. / 6.
    if constexpr (Dimension == 2) {
      eta = 0;
    }
    // else{
    // theta = 1. / 3.;
    // eta   = 1. / 3.;
    // theta = .5;
    // eta   = 0;
    //}
    //
    parallel_for(
      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
        const auto& cell_to_node = cell_to_node_matrix[j];

        Rp SumFluxesNode(zero);

        for (size_t l = 0; l < cell_to_node.size(); ++l) {
          SumFluxesNode += Gjr[j][l];
        }
        // SumFluxesNode *= (1 - theta);

        const auto& cell_to_edge = cell_to_edge_matrix[j];
        Rp SumFluxesEdge(zero);

        for (size_t l = 0; l < cell_to_edge.size(); ++l) {
          SumFluxesEdge += Gje[j][l];
        }

        const auto& cell_to_face = cell_to_face_matrix[j];
        Rp SumFluxesFace(zero);

        for (size_t l = 0; l < cell_to_face.size(); ++l) {
          SumFluxesFace += Gjf[j][l];
        }
        // SumFluxesEdge *= (theta);

        const Rp SumFluxes = (1 - theta - eta) * SumFluxesNode + eta * SumFluxesEdge + theta * SumFluxesFace;

        State[j] -= dt / volumes_cell[j] * SumFluxes;

        rho[j] = State[j][0];
        for (size_t d = 0; d < Dimension; ++d)
          u[j][d] = State[j][1 + d] / rho[j];
        E[j] = State[j][1 + Dimension] / rho[j];
      });

    return std::make_tuple(std::make_shared<const DiscreteFunctionVariant>(rho),
                           std::make_shared<const DiscreteFunctionVariant>(u),
                           std::make_shared<const DiscreteFunctionVariant>(E));
  }

  HybridHLLcRusanovEulerianCompositeSolver_v2()  = default;
  ~HybridHLLcRusanovEulerianCompositeSolver_v2() = default;
};

template <MeshConcept MeshType>
class HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>::WallBoundaryCondition
{
};

template <>
class HybridHLLcRusanovEulerianCompositeSolver_v2<Mesh<2>>::WallBoundaryCondition
{
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshNodeBoundary m_mesh_node_boundary;
  const MeshFaceBoundary m_mesh_face_boundary;
  // const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
  // const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;

 public:
  size_t
  numberOfNodes() const
  {
    return m_mesh_node_boundary.nodeList().size();
  }

  size_t
  numberOfFaces() const
  {
    return m_mesh_face_boundary.faceList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_face_boundary.faceList();
  }

  WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary)
    : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
  {
    ;
  }
};

template <>
class HybridHLLcRusanovEulerianCompositeSolver_v2<Mesh<3>>::WallBoundaryCondition
{
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshNodeBoundary m_mesh_node_boundary;
  const MeshEdgeBoundary m_mesh_edge_boundary;
  const MeshFaceBoundary m_mesh_face_boundary;

 public:
  size_t
  numberOfNodes() const
  {
    return m_mesh_node_boundary.nodeList().size();
  }
  size_t
  numberOfEdges() const
  {
    return m_mesh_edge_boundary.edgeList().size();
  }

  size_t
  numberOfFaces() const
  {
    return m_mesh_face_boundary.faceList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const EdgeId>&
  edgeList() const
  {
    return m_mesh_edge_boundary.edgeList();
  }

  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_face_boundary.faceList();
  }

  WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
                        const MeshEdgeBoundary& mesh_edge_boundary,
                        const MeshFaceBoundary& mesh_face_boundary)
    : m_mesh_node_boundary(mesh_node_boundary),

      m_mesh_edge_boundary(mesh_edge_boundary),

      m_mesh_face_boundary(mesh_face_boundary)
  {
    ;
  }
};

template <MeshConcept MeshType>
class HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>::NeumannflatBoundaryCondition
{
};
template <>
class HybridHLLcRusanovEulerianCompositeSolver_v2<Mesh<2>>::NeumannflatBoundaryCondition
{
 public:
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
  const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;

 public:
  const Rd&
  outgoingNormal() const
  {
    return m_mesh_flat_node_boundary.outgoingNormal();
  }

  size_t
  numberOfNodes() const
  {
    return m_mesh_flat_node_boundary.nodeList().size();
  }

  size_t
  numberOfFaces() const
  {
    return m_mesh_flat_face_boundary.faceList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_flat_node_boundary.nodeList();
  }

  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_flat_face_boundary.faceList();
  }

  NeumannflatBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
                               const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
    : m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary)
  {
    ;
  }

  ~NeumannflatBoundaryCondition() = default;
};

template <>
class HybridHLLcRusanovEulerianCompositeSolver_v2<Mesh<3>>::NeumannflatBoundaryCondition
{
 public:
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
  const MeshFlatEdgeBoundary<MeshType> m_mesh_flat_edge_boundary;
  const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;

 public:
  const Rd&
  outgoingNormal() const
  {
    return m_mesh_flat_node_boundary.outgoingNormal();
  }

  size_t
  numberOfNodes() const
  {
    return m_mesh_flat_node_boundary.nodeList().size();
  }

  size_t
  numberOfEdges() const
  {
    return m_mesh_flat_edge_boundary.edgeList().size();
  }

  size_t
  numberOfFaces() const
  {
    return m_mesh_flat_face_boundary.faceList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_flat_node_boundary.nodeList();
  }

  const Array<const EdgeId>&
  edgeList() const
  {
    return m_mesh_flat_edge_boundary.edgeList();
  }

  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_flat_face_boundary.faceList();
  }

  NeumannflatBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
                               const MeshFlatEdgeBoundary<MeshType>& mesh_flat_edge_boundary,
                               const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
    : m_mesh_flat_node_boundary(mesh_flat_node_boundary),
      m_mesh_flat_edge_boundary(mesh_flat_edge_boundary),
      m_mesh_flat_face_boundary(mesh_flat_face_boundary)
  {
    ;
  }

  ~NeumannflatBoundaryCondition() = default;
};

template <MeshConcept MeshType>
class HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>::SymmetryBoundaryCondition
{
};

template <>
class HybridHLLcRusanovEulerianCompositeSolver_v2<Mesh<2>>::SymmetryBoundaryCondition
{
 public:
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
  const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;

 public:
  const Rd&
  outgoingNormal() const
  {
    return m_mesh_flat_node_boundary.outgoingNormal();
  }

  size_t
  numberOfNodes() const
  {
    return m_mesh_flat_node_boundary.nodeList().size();
  }

  size_t
  numberOfFaces() const
  {
    return m_mesh_flat_face_boundary.faceList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_flat_node_boundary.nodeList();
  }

  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_flat_face_boundary.faceList();
  }

  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
                            const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
    : m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary)
  {
    ;
  }

  ~SymmetryBoundaryCondition() = default;
};

template <>
class HybridHLLcRusanovEulerianCompositeSolver_v2<Mesh<3>>::SymmetryBoundaryCondition
{
 public:
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
  const MeshFlatEdgeBoundary<MeshType> m_mesh_flat_edge_boundary;
  const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;

 public:
  const Rd&
  outgoingNormal() const
  {
    return m_mesh_flat_node_boundary.outgoingNormal();
  }

  size_t
  numberOfNodes() const
  {
    return m_mesh_flat_node_boundary.nodeList().size();
  }

  size_t
  numberOfEdges() const
  {
    return m_mesh_flat_edge_boundary.edgeList().size();
  }

  size_t
  numberOfFaces() const
  {
    return m_mesh_flat_face_boundary.faceList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_flat_node_boundary.nodeList();
  }

  const Array<const EdgeId>&
  edgeList() const
  {
    return m_mesh_flat_edge_boundary.edgeList();
  }

  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_flat_face_boundary.faceList();
  }

  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
                            const MeshFlatEdgeBoundary<MeshType>& mesh_flat_edge_boundary,
                            const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
    : m_mesh_flat_node_boundary(mesh_flat_node_boundary),
      m_mesh_flat_edge_boundary(mesh_flat_edge_boundary),
      m_mesh_flat_face_boundary(mesh_flat_face_boundary)
  {
    ;
  }

  ~SymmetryBoundaryCondition() = default;
};

template <MeshConcept MeshType>
class HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>::InflowListBoundaryCondition
{
};

template <>
class HybridHLLcRusanovEulerianCompositeSolver_v2<Mesh<2>>::InflowListBoundaryCondition
{
 public:
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshNodeBoundary m_mesh_node_boundary;
  const MeshFaceBoundary m_mesh_face_boundary;
  const Table<const double> m_node_array_list;
  const Table<const double> m_face_array_list;

 public:
  size_t
  numberOfNodes() const
  {
    return m_mesh_node_boundary.nodeList().size();
  }

  size_t
  numberOfFaces() const
  {
    return m_mesh_face_boundary.faceList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_face_boundary.faceList();
  }

  const Table<const double>&
  nodeArrayList() const
  {
    return m_node_array_list;
  }

  const Table<const double>&
  faceArrayList() const
  {
    return m_face_array_list;
  }

  InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
                              const MeshFaceBoundary& mesh_face_boundary,
                              const Table<const double>& node_array_list,
                              const Table<const double>& face_array_list)
    : m_mesh_node_boundary(mesh_node_boundary),
      m_mesh_face_boundary(mesh_face_boundary),
      m_node_array_list(node_array_list),
      m_face_array_list(face_array_list)
  {
    ;
  }

  ~InflowListBoundaryCondition() = default;
};

template <>
class HybridHLLcRusanovEulerianCompositeSolver_v2<Mesh<3>>::InflowListBoundaryCondition
{
 public:
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshNodeBoundary m_mesh_node_boundary;
  const MeshEdgeBoundary m_mesh_edge_boundary;
  const MeshFaceBoundary m_mesh_face_boundary;
  const Table<const double> m_node_array_list;
  const Table<const double> m_edge_array_list;
  const Table<const double> m_face_array_list;

 public:
  size_t
  numberOfNodes() const
  {
    return m_mesh_node_boundary.nodeList().size();
  }

  size_t
  numberOfEdges() const
  {
    return m_mesh_edge_boundary.edgeList().size();
  }

  size_t
  numberOfFaces() const
  {
    return m_mesh_face_boundary.faceList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const EdgeId>&
  edgeList() const
  {
    return m_mesh_edge_boundary.edgeList();
  }

  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_face_boundary.faceList();
  }

  const Table<const double>&
  nodeArrayList() const
  {
    return m_node_array_list;
  }

  const Table<const double>&
  edgeArrayList() const
  {
    return m_edge_array_list;
  }

  const Table<const double>&
  faceArrayList() const
  {
    return m_face_array_list;
  }

  InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
                              const MeshEdgeBoundary& mesh_edge_boundary,
                              const MeshFaceBoundary& mesh_face_boundary,
                              const Table<const double>& node_array_list,
                              const Table<const double>& edge_array_list,
                              const Table<const double>& face_array_list)
    : m_mesh_node_boundary(mesh_node_boundary),
      m_mesh_edge_boundary(mesh_edge_boundary),
      m_mesh_face_boundary(mesh_face_boundary),
      m_node_array_list(node_array_list),
      m_edge_array_list(edge_array_list),
      m_face_array_list(face_array_list)
  {
    ;
  }

  ~InflowListBoundaryCondition() = default;
};

template <MeshConcept MeshType>
class HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>::OutflowBoundaryCondition
{
};

template <>
class HybridHLLcRusanovEulerianCompositeSolver_v2<Mesh<2>>::OutflowBoundaryCondition
{
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshNodeBoundary m_mesh_node_boundary;
  const MeshFaceBoundary m_mesh_face_boundary;

 public:
  size_t
  numberOfNodes() const
  {
    return m_mesh_node_boundary.nodeList().size();
  }

  size_t
  numberOfFaces() const
  {
    return m_mesh_face_boundary.faceList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_face_boundary.faceList();
  }

  OutflowBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary)
    : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
  {
    ;
  }
};

template <>
class HybridHLLcRusanovEulerianCompositeSolver_v2<Mesh<3>>::OutflowBoundaryCondition
{
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshNodeBoundary m_mesh_node_boundary;
  const MeshEdgeBoundary m_mesh_edge_boundary;
  const MeshFaceBoundary m_mesh_face_boundary;

 public:
  size_t
  numberOfNodes() const
  {
    return m_mesh_node_boundary.nodeList().size();
  }
  size_t
  numberOfEdges() const
  {
    return m_mesh_edge_boundary.edgeList().size();
  }

  size_t
  numberOfFaces() const
  {
    return m_mesh_face_boundary.faceList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const EdgeId>&
  edgeList() const
  {
    return m_mesh_edge_boundary.edgeList();
  }

  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_face_boundary.faceList();
  }

  OutflowBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
                           const MeshEdgeBoundary& mesh_edge_boundary,
                           const MeshFaceBoundary& mesh_face_boundary)
    : m_mesh_node_boundary(mesh_node_boundary),

      m_mesh_edge_boundary(mesh_edge_boundary),

      m_mesh_face_boundary(mesh_face_boundary)
  {
    ;
  }
};

std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
           std::shared_ptr<const DiscreteFunctionVariant>,
           std::shared_ptr<const DiscreteFunctionVariant>>
hybridHLLcRusanovEulerianCompositeSolver_v2(
  const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
  const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
  const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
  const double& gamma,
  const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
  const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
  //  const size_t& degree,
  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
  const double& dt,
  const bool check)
{
  std::shared_ptr mesh_v = getCommonMesh({rho_v, u_v, E_v, c_v, p_v});
  if (not mesh_v) {
    throw NormalError("discrete functions are not defined on the same mesh");
  }

  if (not checkDiscretizationType({rho_v, u_v, E_v}, DiscreteFunctionType::P0)) {
    throw NormalError("acoustic solver expects P0 functions");
  }

  return std::visit(
    PUGS_LAMBDA(auto&& p_mesh)
      ->std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
                   std::shared_ptr<const DiscreteFunctionVariant>> {
        using MeshType                    = mesh_type_t<decltype(p_mesh)>;
        static constexpr size_t Dimension = MeshType::Dimension;
        using Rd                          = TinyVector<Dimension>;

        if constexpr (Dimension == 1) {
          throw NormalError("Hybrid HLLc/Rusanov EulerianCompositeSolver v2 is not available in 1D");
        } else {
          if constexpr (is_polygonal_mesh_v<MeshType>) {
            return HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>{}
              .solve(p_mesh, rho_v->get<DiscreteFunctionP0<const double>>(), u_v->get<DiscreteFunctionP0<const Rd>>(),
                     E_v->get<DiscreteFunctionP0<const double>>(), gamma, c_v->get<DiscreteFunctionP0<const double>>(),
                     p_v->get<DiscreteFunctionP0<const double>>(), bc_descriptor_list, dt, check);
          } else {
            throw NormalError("Hybrid HLLc/Rusanov EulerianCompositeSolver v2 is only defined on polygonal meshes");
          }
        }
      },
    mesh_v->variant());
}