diff --git a/src/language/modules/CompositeSchemeOtherFluxesModule.cpp b/src/language/modules/CompositeSchemeOtherFluxesModule.cpp
index b6e84023d0a988f8ba5c149f09fbef6001b8c7ea..c3612535d13cbf6edc4f21c0772ba2e8710307d7 100644
--- a/src/language/modules/CompositeSchemeOtherFluxesModule.cpp
+++ b/src/language/modules/CompositeSchemeOtherFluxesModule.cpp
@@ -20,6 +20,7 @@
 #include <scheme/HybridRoeRusanovViscousFormEulerianCompositeSolver_v2.hpp>
 #include <scheme/HybridRusanovRoeViscousFormEulerianCompositeSolver_v2.hpp>
 #include <scheme/HybridHLLRusanovEulerianCompositeSolver_v2.hpp>
+#include <scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.hpp>
 
 
 CompositeSchemeOtherFluxesModule::CompositeSchemeOtherFluxesModule() {
@@ -79,6 +80,25 @@ CompositeSchemeOtherFluxesModule::CompositeSchemeOtherFluxesModule() {
                                                   bc_descriptor_list, dt);
           }
     
+          )); 
+
+    this->_addBuiltinFunction("hybrid_hllc_rusanov_eulerian_composite_solver_version2",
+        std::function(
+    
+          [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+             const std::shared_ptr<const DiscreteFunctionVariant>& u,
+             const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
+             const std::shared_ptr<const DiscreteFunctionVariant>& c,
+             const std::shared_ptr<const DiscreteFunctionVariant>& p,   // const size_t& degree,
+             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+               bc_descriptor_list,
+             const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                             std::shared_ptr<const DiscreteFunctionVariant>,
+                                             std::shared_ptr<const DiscreteFunctionVariant>> {
+            return hybridHLLcRusanovEulerianCompositeSolver_v2(rho, u, E, gamma, c, p,   // degree,
+                                                  bc_descriptor_list, dt);
+          }
+    
           )); 
 }
 
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 74f878fccd77fef0a4a97d0e11ec5b3918ce4649..5c645eb9aa3f952b736a0fe01097d8874e01e569 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -35,6 +35,7 @@ add_library(
   HybridRoeRusanovViscousFormEulerianCompositeSolver_v2.cpp
   HybridRusanovRoeViscousFormEulerianCompositeSolver_v2.cpp
   HybridHLLRusanovEulerianCompositeSolver_v2.cpp
+  HybridHLLcRusanovEulerianCompositeSolver_v2.cpp
 )
 
 target_link_libraries(
diff --git a/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp b/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0a7c00c6d095111f14aa46515628a7c947bfe82b
--- /dev/null
+++ b/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp
@@ -0,0 +1,2072 @@
+#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());
+}
diff --git a/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.hpp b/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..124427df443d149e16bc4312e292121876a18c7b
--- /dev/null
+++ b/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.hpp
@@ -0,0 +1,30 @@
+#ifndef HYBRID_HLLC_RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+#define HYBRID_HLLC_RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+
+#include <memory>
+#include <vector>
+
+// double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+//                  const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
+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,
+  const std::shared_ptr<const DiscreteFunctionVariant>& u,
+  const std::shared_ptr<const DiscreteFunctionVariant>& E,
+  const double& gamma,
+  const std::shared_ptr<const DiscreteFunctionVariant>& c,
+  const std::shared_ptr<const DiscreteFunctionVariant>& p,
+  // const size_t& degree,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const double& dt,
+  const bool check = false);
+
+#endif   // HYBRID_HLLC_RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP