From edf6bb03d16763db90a0337848cab1868341d8a5 Mon Sep 17 00:00:00 2001
From: HOCH PHILIPPE <philippe.hoch@gmail.com>
Date: Tue, 6 May 2025 23:27:12 +0200
Subject: [PATCH] Adding Roe composite under Flux Form 2D and 3D

---
 src/scheme/CMakeLists.txt                     |    1 +
 src/scheme/CellbyCellLimitation.hpp           |  264 +-
 .../RoeFluxFormEulerianCompositeSolver_v2.cpp | 2846 +++++++++++++++++
 .../RoeFluxFormEulerianCompositeSolver_v2.hpp |   89 +
 4 files changed, 3175 insertions(+), 25 deletions(-)
 create mode 100644 src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp
 create mode 100644 src/scheme/RoeFluxFormEulerianCompositeSolver_v2.hpp

diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 009bb3bf1..6e79e4324 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -28,6 +28,7 @@ add_library(
   RusanovEulerianCompositeSolver.cpp
   RusanovEulerianCompositeSolver_v2.cpp
   RoeViscousFormEulerianCompositeSolver_v2.cpp
+  RoeFluxFormEulerianCompositeSolver_v2.cpp
   RusanovEulerianCompositeSolver_o2.cpp
   RusanovEulerianCompositeSolver_v2_o2.cpp
   RusanovEulerianCompositeSolver_v2_order_n.cpp
diff --git a/src/scheme/CellbyCellLimitation.hpp b/src/scheme/CellbyCellLimitation.hpp
index a443e71ef..ecd607a0f 100644
--- a/src/scheme/CellbyCellLimitation.hpp
+++ b/src/scheme/CellbyCellLimitation.hpp
@@ -1,9 +1,6 @@
 #ifndef CELL_BY_CELL_LIMITATION_HPP
 #define CELL_BY_CELL_LIMITATION_HPP
 
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-
 #include <language/utils/InterpolateItemArray.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
@@ -24,6 +21,7 @@
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/LimitationTools.hpp>
 #include <scheme/PolynomialReconstruction.hpp>
 #include <scheme/PolynomialReconstructionDescriptor.hpp>
 #include <utils/PugsTraits.hpp>
@@ -44,7 +42,8 @@ class CellByCellLimitation
   density_limiter(const MeshType& mesh,
                   const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
                   const DiscreteFunctionP0<const double>& rho,
-                  DiscreteFunctionDPk<Dimension, double>& rho_L) const
+                  DiscreteFunctionDPk<Dimension, double>& rho_L,
+                  const bool enableWeakBoundPositivityOnly = false) const
   {
     const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
     const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
@@ -95,9 +94,13 @@ class CellByCellLimitation
           for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
             const FaceId face_id = face_list[i_face];
 
-            const Rd& x0 = xr[face_to_node_matrix[face_id][0]];
-            const Rd& x2 = xr[face_to_node_matrix[face_id][1]];
-            const Rd& x1 = .5 * (x0 + x2);
+            const auto& face_to_node = face_to_node_matrix[face_id];
+            const NodeId& nodei      = face_to_node[0];
+            const NodeId& nodef      = face_to_node[face_to_node.size() - 1];
+
+            const Rd x0 = xr[nodei];
+            const Rd x1 = xl[face_id];
+            const Rd x2 = xr[nodef];
 
             const double rho_x0 = rho_L[cell_id](x0);
             const double rho_x1 = rho_L[cell_id](x1);
@@ -106,16 +109,120 @@ class CellByCellLimitation
             rho_bar_min = std::min(rho_bar_min, std::min(rho_x0, std::min(rho_x1, rho_x2)));
             rho_bar_max = std::max(rho_bar_max, std::max(rho_x0, std::max(rho_x1, rho_x2)));
           }
-          // const LineParabolicTransformation<Dimension> t(x0, x1, x2);
-          // for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) {
-          // const double rho_xk = rho_L[cell_id](t(qf.point(iq)));
-
-          // rho_bar_min = std::min(rho_bar_min, rho_xk);
-          // rho_bar_max = std::max(rho_bar_max, rho_xk);
-          //}
         } else {
+          throw NotImplementedError("not implement in 3D");
+        }
+        /*
+        const double eps = 1E-14;
+        double coef1     = 1;
+        if (std::abs(rho_bar_max - rhoj) > eps) {
+          coef1 = (rho_max - rhoj) / ((rho_bar_max - rhoj));
+        }
+
+        double coef2 = 1.;
+        if (std::abs(rho_bar_min - rhoj) > eps) {
+          coef2 = (rho_min - rhoj) / ((rho_bar_min - rhoj));
+        }
+*/
+        const double lambda =   // std::max(0., std::min(1., std::min(coef1, coef2)));
+          toolsLimitation::computeCoefLimitationBarthJespersen(rhoj, rho_bar_min, rho_bar_max, rho_min, rho_max,
+                                                               enableWeakBoundPositivityOnly);
+
+        auto coefficients = rho_L.coefficients(cell_id);
+
+        coefficients[0] = (1 - lambda) * rho[cell_id] + lambda * coefficients[0];
+
+        for (size_t i = 1; i < coefficients.size(); ++i) {
+          coefficients[i] *= lambda;
+        }
+      });
+  }
+
+  void
+  density_limiter(const MeshType& mesh,
+                  const FaceArray<TinyVector<2>>& QuadratureFace,
+                  const EdgeArray<TinyVector<2>>& QuadratureEdge,
+                  const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
+                  const DiscreteFunctionP0<const double>& rho,
+                  DiscreteFunctionDPk<Dimension, double>& rho_L,
+                  const bool enableWeakBoundPositivityOnly = false) const
+  {
+    const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+
+    const auto& xr = mesh.xr();
+    // const auto& xl = mesh.xl();
+
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    //    auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1);
+    auto stencil = StencilManager::instance()
+                     .getCellToCellStencilArray(mesh.connectivity(),
+                                                StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
+                                                symmetry_boundary_descriptor_list);
+    auto xj        = mesh_data.xj();
+    const auto& xl = mesh_data.xl();
+
+    // const QuadratureFormula<1> qf =
+    //   QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        const double rhoj = rho[cell_id];
+
+        double rho_min = rhoj;
+        double rho_max = rhoj;
+
+        const auto cell_stencil = stencil[cell_id];
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          rho_min = std::min(rho_min, rho[cell_stencil[i_cell]]);
+          rho_max = std::max(rho_max, rho[cell_stencil[i_cell]]);
         }
 
+        double rho_bar_min = rhoj;
+        double rho_bar_max = rhoj;
+
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          const CellId cell_k_id = cell_stencil[i_cell];
+          const double rho_xk    = rho_L[cell_id](xj[cell_k_id]);
+
+          rho_bar_min = std::min(rho_bar_min, rho_xk);
+          rho_bar_max = std::max(rho_bar_max, rho_xk);
+        }
+
+        if constexpr (Dimension == 2) {
+          auto face_list = cell_to_face_matrix[cell_id];
+
+          for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+            const FaceId face_id = face_list[i_face];
+
+            const auto& face_to_node = face_to_node_matrix[face_id];
+            const NodeId& nodei      = face_to_node[0];
+            const NodeId& nodef      = face_to_node[face_to_node.size() - 1];
+
+            const Rd& x0 = xr[nodei];
+            const Rd& x1 = xr[nodef];
+
+            const double rho_x0 = rho_L[cell_id](x0);
+            const double rho_x1 = rho_L[cell_id](x1);
+
+            rho_bar_min = std::min(rho_bar_min, std::min(rho_x0, rho_x1));
+            rho_bar_max = std::min(rho_bar_max, std::max(rho_x0, rho_x1));
+
+            for (size_t i = 0; i < QuadratureFace[face_id].size(); ++i) {
+              // const Rd& x_quad = xd + QuadratureFace[face][i][1] * (xf - xd);
+              const Rd& x_quad = xr[nodei] + QuadratureFace[face_id][i][1] * (xr[nodef] - xr[nodei]);
+
+              const double rho_xk = rho_L[cell_id](x_quad);
+
+              rho_bar_min = std::min(rho_bar_min, rho_xk);
+              rho_bar_max = std::max(rho_bar_max, rho_xk);
+            }
+          }
+        } else {
+          throw NotImplementedError("not implement in 3D");
+        }
+        /*
         const double eps = 1E-14;
         double coef1     = 1;
         if (std::abs(rho_bar_max - rhoj) > eps) {
@@ -126,8 +233,11 @@ class CellByCellLimitation
         if (std::abs(rho_bar_min - rhoj) > eps) {
           coef2 = (rho_min - rhoj) / ((rho_bar_min - rhoj));
         }
+*/
+        const double lambda =   // std::max(0., std::min(1., std::min(coef1, coef2)));
 
-        const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
+          toolsLimitation::computeCoefLimitationBarthJespersen(rhoj, rho_bar_min, rho_bar_max, rho_min, rho_max,
+                                                               enableWeakBoundPositivityOnly);
 
         auto coefficients = rho_L.coefficients(cell_id);
 
@@ -139,26 +249,27 @@ class CellByCellLimitation
       });
   }
 
+  //
+  // Idem for epsilon
+  //
+
   void
   specific_internal_nrj_limiter(
     const MeshType& mesh,
     const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
-
-    // const DiscreteFunctionP0<const double>& rho,
-    // const DiscreteFunctionDPk<Dimension, double>& rho_L,
     const DiscreteFunctionP0<const double>& epsilon,
-    // const DiscreteFunctionDPk<Dimension, double>&
-    auto epsilon_R(const CellId cell_id, const Rd& x),
-    CellValue<double>& lambda_epsilon)   // const
-  //                                DiscreteFunctionDPk<Dimension, double>& epsilon_R) const
+    auto epsilon_R,   //(const CellId cell_id, const Rd& x),
+    CellValue<double>& lambda_epsilon,
+    const bool enableWeakBoundPositivityOnly = false) const
   {
     const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
     const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
 
     const auto& xr = mesh.xr();
-    const auto& xl = mesh.xl();
+    // const auto& xl = mesh.xl();
 
     MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    const auto& xl                = mesh_data.xl();
 
     //    auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1);
     auto stencil = StencilManager::instance()
@@ -216,7 +327,7 @@ class CellByCellLimitation
           // epsilon_R_max = std::max(epsilon_R_max, epsilon_xk);
           //}
         }
-
+        /*
         const double eps = 1E-14;
         double coef1     = 1;
         if (std::abs(epsilon_R_max - epsilonj) > eps) {
@@ -227,11 +338,114 @@ class CellByCellLimitation
         if (std::abs(epsilon_R_min - epsilonj) > eps) {
           coef2 = (epsilon_min - epsilonj) / ((epsilon_R_min - epsilonj));
         }
+*/
+        lambda_epsilon[cell_id] =   // std::max(0., std::min(1., std::min(coef1, coef2)));
+          toolsLimitation::computeCoefLimitationBarthJespersen(epsilonj, epsilon_R_min, epsilon_R_min, epsilon_min,
+                                                               epsilon_max, enableWeakBoundPositivityOnly);
+      });
+  }
 
-        lambda_epsilon[cell_id] = std::max(0., std::min(1., std::min(coef1, coef2)));
+  //
+  //
+  //
+  void
+  specific_internal_nrj_limiter(
+    const MeshType& mesh,
+    const FaceArray<TinyVector<2>>& QuadratureFace,
+    const EdgeArray<TinyVector<2>>& QuadratureEdge,
+    const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
+    const DiscreteFunctionP0<const double>& epsilon,
+    auto epsilon_R,   //(const CellId cell_id, const Rd& x),
+    CellValue<double>& lambda_epsilon,
+    const bool enableWeakBoundPositivityOnly = false) const
+  {
+    const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+
+    const auto& xr = mesh.xr();
+    // const auto& xl = mesh.xl();
+
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    //    auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1);
+    auto stencil = StencilManager::instance()
+                     .getCellToCellStencilArray(mesh.connectivity(),
+                                                StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
+                                                symmetry_boundary_descriptor_list);
+    auto xj = mesh_data.xj();
+
+    //    const QuadratureFormula<1> qf =
+    // QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        const double epsilonj = epsilon[cell_id];
+
+        double epsilon_min = epsilonj;
+        double epsilon_max = epsilonj;
+
+        const auto cell_stencil = stencil[cell_id];
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          epsilon_min = std::min(epsilon_min, epsilon[cell_stencil[i_cell]]);
+          epsilon_max = std::max(epsilon_max, epsilon[cell_stencil[i_cell]]);
+        }
+
+        double epsilon_R_min = epsilonj;
+        double epsilon_R_max = epsilonj;
+
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          const CellId cell_k_id  = cell_stencil[i_cell];
+          const double epsilon_xk = epsilon_R(cell_id, xj[cell_k_id]);
+
+          epsilon_R_min = std::min(epsilon_R_min, epsilon_xk);
+          epsilon_R_max = std::max(epsilon_R_max, epsilon_xk);
+        }
+
+        auto face_list = cell_to_face_matrix[cell_id];
+        for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+          const FaceId face_id = face_list[i_face];
+
+          const Rd& x0 = xr[face_to_node_matrix[face_id][0]];
+          // const Rd& x1            = xl[face_id][0];
+          const Rd& x2            = xr[face_to_node_matrix[face_id][1]];
+          const double epsilon_x0 = epsilon_R(cell_id, x0);
+          // const double epsilon_x1 = epsilon_R(cell_id, x1);
+          const double epsilon_x2 = epsilon_R(cell_id, x2);
+          epsilon_R_min           = std::min(epsilon_R_min, std::min(epsilon_x0, epsilon_x2));
+          epsilon_R_max           = std::max(epsilon_R_max, std::max(epsilon_x0, epsilon_x2));
+
+          for (size_t i = 0; i < QuadratureFace[face_id].size(); ++i) {
+            // const Rd& x_quad = xd + QuadratureFace[face][i][1] * (xf - xd);
+            const Rd& x_quad = x0 + QuadratureFace[face_id][i][1] * (x2 - x0);
+
+            const double epsilon_xk = epsilon_R(cell_id, x_quad);
+
+            epsilon_R_min = std::min(epsilon_R_min, epsilon_xk);
+            epsilon_R_max = std::max(epsilon_R_max, epsilon_xk);
+          }
+        }
+        /*
+        const double eps = 1E-14;
+        double coef1     = 1;
+        if (std::abs(epsilon_R_max - epsilonj) > eps) {
+          coef1 = (epsilon_max - epsilonj) / ((epsilon_R_max - epsilonj));
+        }
+
+        double coef2 = 1.;
+        if (std::abs(epsilon_R_min - epsilonj) > eps) {
+          coef2 = (epsilon_min - epsilonj) / ((epsilon_R_min - epsilonj));
+        }
+*/
+        lambda_epsilon[cell_id] =   // std::max(0., std::min(1., std::min(coef1, coef2)));
+          toolsLimitation::computeCoefLimitationBarthJespersen(epsilonj, epsilon_R_min, epsilon_R_max, epsilon_min,
+                                                               epsilon_max, enableWeakBoundPositivityOnly);
       });
   }
 
+  //
+  // Other version (CellValue Coef for limitation )
+  //
+
   void
   computeLimitorVolumicScalarQuantityMinModDukowiczGradient(const MeshType& mesh,
                                                             const CellValue<double>& q,
diff --git a/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp
new file mode 100644
index 000000000..999d343f2
--- /dev/null
+++ b/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp
@@ -0,0 +1,2846 @@
+#include <scheme/RoeFluxFormEulerianCompositeSolver_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 RoeFluxFormEulerianCompositeSolver_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>;
+
+  using Rpxd = TinyMatrix<Dimension + 2, Dimension>;
+
+  class SymmetryBoundaryCondition;
+  class InflowListBoundaryCondition;
+  class OutflowBoundaryCondition;
+  class WallBoundaryCondition;
+  class NeumannflatBoundaryCondition;
+
+  using BoundaryCondition = std::variant<SymmetryBoundaryCondition,
+                                         InflowListBoundaryCondition,
+                                         OutflowBoundaryCondition,
+                                         NeumannflatBoundaryCondition,
+                                         WallBoundaryCondition>;
+
+  using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+  Rpxp
+  SignofASquareDiagonalisableMatrix(const Rpxp& Left, const Rp& Diag_vp, const Rpxp& Right) const
+  {
+    Rpxp MatriceDiag(identity);
+    for (size_t i = 0; i < Left.numberOfRows(); ++i)
+      MatriceDiag(i, i) *= signe(Diag_vp[i]);
+
+    return Left * (MatriceDiag * Right);
+
+    // return (Left * Diag_vp * Right);
+    Rpxp M(zero);
+
+    for (size_t i = 0; i < Left.numberOfRows(); ++i)
+      for (size_t j = 0; j < Left.numberOfColumns(); ++j) {
+        M(i, j) = signe(Diag_vp[i]) * Right(i, j);
+      }
+    return Left * M;
+  };
+
+  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 Roe v2 Eulerian Composite solver";
+        throw NormalError(error_msg.str());
+      }
+    }
+
+    return bc_list;
+  }
+
+ public:
+  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
+  _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] = 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);
+
+              // 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];
+            }
+
+            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];
+                }
+              }
+            }
+          }
+        },
+        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;
+  }
+
+  inline Rpxd
+  Flux(const double& rho, const Rd& U, const double& E, const double gam) const
+  {
+    // const R2 flux_rho   = rhoU;
+    // const R22 flux_rhoU = R22(rhoU.x1() * rhoU.x1() / rho + P, rhoU.x1() * rhoU.x2() / rho, rhoU.x2() * rhoU.x1()
+    // rho,rhoU.x2() * rhoU.x2() / rho + P);
+    // const R2 flux_rhoE  = ((rhoE + P) / rho) * rhoU;
+    /*  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;
+
+    */
+    const Rd& rhoU       = rho * U;
+    const Rdxd rhoUTensU = tensorProduct(rhoU, U);
+    const double p       = pression(rho, E - .5 * dot(U, U), gam);
+
+    Rdxd pid(identity);
+    pid *= p;
+    const Rdxd rhoUTensUPlusPid = rhoUTensU + pid;
+    const double rhoEPlusP      = rho * E + p;
+
+    const Rd& rhoEPlusPtimesU = rhoEPlusP * U;
+
+    Rpxd Fluxx;   // En ligne ci dessous
+
+    Fluxx[0] = rhoU;
+    for (size_t dim = 0; dim < Dimension; ++dim)
+      // for (size_t dim2 = 0; dim2 < Dimension; ++dim2)
+      Fluxx[1 + dim] = rhoUTensUPlusPid[dim];
+    // Fluxx[1 + dim][dim2] = rhoUTensUPlusPid[dim][dim2];
+    Fluxx[1 + Dimension] = rhoEPlusPtimesU;
+    return Fluxx;
+  }
+
+  struct JacobianInformations
+  {
+    Rpxp Jacobian;
+    Rpxp LeftEigenVectors;
+    Rpxp RightEigenVectors;
+    Rp EigenValues;
+    Rpxp AbsJacobian;
+    double maxabsVpNormal;
+    // bool changingSign;
+    double MaxErrorProd;
+
+   public:
+    JacobianInformations(const Rpxp& Jacobiand,
+                         const Rpxp& LeftEigenVectorsd,
+                         const Rpxp& RightEigenVectorsd,
+                         const Rp& EigenValuesd,
+                         const Rpxp& AbsJacobiand,
+                         const double maxabsVpNormald,
+                         const double MaxErrorProdd)
+      : Jacobian(Jacobiand),
+        LeftEigenVectors(LeftEigenVectorsd),
+        RightEigenVectors(RightEigenVectorsd),
+        EigenValues(EigenValuesd),
+        AbsJacobian(AbsJacobiand),
+        maxabsVpNormal(maxabsVpNormald),
+        MaxErrorProd(MaxErrorProdd)
+    {}
+    ~JacobianInformations() {}
+  };
+
+  Rp
+  EvaluateEigenValuesNormal(   // const double rhod,   // rhoJ, uJ, EJ, gammaJ, cJ, pJ, normal
+    const Rd& Ud,
+    // const double Ed,
+    // const double gammad,
+    const double cd,
+    // const double pd,
+    const Rd& normal) const
+  {
+    Rp Eigen;
+    const double uscaln = dot(Ud, normal);
+
+    Eigen[0] = uscaln - cd;
+    for (size_t dim = 0; dim < Dimension; ++dim)
+      Eigen[dim] = uscaln;
+    Eigen[1 + Dimension] = uscaln + cd;
+
+    return Eigen;
+  }
+
+  enum TypeAverageState
+  {
+    RoeAverage,
+    MeanAverage
+  };
+  struct AverageStateStructData
+  {
+    double rho;
+    Rd U;
+    double E;
+    double H;
+    double gamma;
+    double p;
+    double c;
+    AverageStateStructData(const double rhod,
+                           const Rd Ud,
+                           const double Ed,
+                           const double Hd,
+                           const double gammad,
+                           const double pd,
+                           const double cd)
+      : rho(rhod), U(Ud), E(Ed), H(Hd), gamma(gammad), p(pd), c(cd)
+    {}
+  };
+
+  AverageStateStructData
+  RoeAverageState(const double& rhoG,
+                  const Rd& UG,
+                  const double& EG,
+                  const double& gammaG,
+                  const double& pG,
+                  const double& rhoD,
+                  const Rd& UD,
+                  const double& ED,
+                  const double gammaD,
+                  const double pD) const
+
+  {
+    double gamma              = .5 * (gammaG + gammaD);   // ou ponderation racine roG et roD
+    double RacineRoG          = sqrt(rhoG);
+    double RacineRoD          = sqrt(rhoD);
+    double rho_mean           = RacineRoG * RacineRoD;
+    Rd U_mean                 = 1. / (RacineRoG + RacineRoD) * (RacineRoG * UG + RacineRoD * UD);
+    double unorm2             = dot(U_mean, U_mean);
+    double NrjCin             = .5 * unorm2;
+    const double TotEnthalpyG = (EG + pG / rhoG);
+    const double TotEnthalpyD = (ED + pD / rhoD);
+
+    double H_mean = (RacineRoG * TotEnthalpyG + RacineRoD * TotEnthalpyD) / (RacineRoG + RacineRoD);
+
+    double E_mean = H_mean / gamma + ((gamma - 1) / (gamma)) * (NrjCin);
+
+    double P_mean = rho_mean * (H_mean - E_mean);
+
+    double c2 = gamma * P_mean / rho_mean;   // cspeed*cspeed;
+    // assert(fabs((gamma - 1) * rho_mean * (E_mean - .5 * (u_mean, u_mean)) - P_mean) < 1e-13);   // equilibre GP
+    double c_mean = sqrt(c2);   // cspeed_meandof/Area;
+
+    return AverageStateStructData(rho_mean, U_mean, E_mean, H_mean, gamma, P_mean, c_mean);
+  }
+
+  AverageStateStructData
+  MeanAverageState(const double& rhoG,
+                   const Rd& UG,
+                   const double& EG,
+                   const double& gammaG,
+                   const double& pG,
+                   const double& rhoD,
+                   const Rd& UD,
+                   const double& ED,
+                   const double gammaD,
+                   const double pD) const
+
+  {
+    double gamma = .5 * (gammaG + gammaD);   // ou ponderation racine roG et roD
+
+    double rho_mean           = .5 * (rhoG + rhoD);
+    Rd U_mean                 = .5 / (rhoG + rhoD) * (rhoG * UG + rhoD * UD);
+    double unorm2             = dot(U_mean, U_mean);
+    double NrjCin             = .5 * unorm2;
+    const double TotEnthalpyG = (EG + pG / rhoG);
+    const double TotEnthalpyD = (ED + pD / rhoD);
+    // double H_mean = (RacineRoG * TotEnthalpyG + RacineRoD * TotEnthalpyD) / (RacineRoG + RacineRoD);
+
+    double H_mean = .5 / (rhoG + rhoD) * (rhoG * TotEnthalpyG + rhoD * TotEnthalpyD);
+
+    double E_mean = H_mean / gamma + ((gamma - 1) / (gamma)) * (NrjCin);
+
+    double P_mean = rho_mean * (H_mean - E_mean);
+
+    double c2 = gamma * P_mean / rho_mean;   // cspeed*cspeed;
+    // assert(fabs((gamma - 1) * rho_mean * (E_mean - .5 * (u_mean, u_mean)) - P_mean) < 1e-13);   // equilibre GP
+    double c_mean = sqrt(c2);   // cspeed_meandof/Area;
+
+    return AverageStateStructData(rho_mean, U_mean, E_mean, H_mean, gamma, P_mean, c_mean);
+  }
+
+  bool
+  EvaluateChangingSignVpAlongNormal(   // const double rhoJ,
+    const Rd& uJ,
+    // const double EJ,
+    // const double gammaJ,
+    const double cJ,
+    // const double pJ,
+
+    // const double rhoK,
+    const Rd& uK,
+    // const double EK,
+    // const double gammaK,
+    const double cK,
+    // const double pK,
+    const Rd& normal) const
+  {
+    const Rp& EigenJ = EvaluateEigenValuesNormal(   // rhoJ,
+      uJ,                                           // EJ, gammaJ,
+      cJ,                                           // pJ,
+      normal);
+
+    const Rp& EigenK = EvaluateEigenValuesNormal(   // rhoK,
+      uK,                                           // EK, gammaK,
+      cK,                                           // pK,
+      normal);
+
+    if (EigenJ[0] * EigenK[0] < -1e-12)   // <=0)              //< -1e-12)
+      return true;
+    for (size_t dim = 1; dim < Dimension + 1; ++dim)   // vp multiplicite d
+      if (EigenJ[dim] * EigenK[dim] < -1e-12)          // <=0)              //< -1e-12)
+        return true;
+    if (EigenJ[1 + Dimension] * EigenK[1 + Dimension] < -1e-12)   // <=0)              //< -1e-12)
+      return true;
+
+    return false;
+  }
+
+  JacobianInformations
+  JacobianFluxAlongUnitNormal(const AverageStateStructData& RoeState,
+                              /*
+    const double rhoJ,
+    const Rd& uJ,
+    const double EJ,
+    const double gammaJ,
+    const double cJ,
+    const double pJ,
+
+    const double rhoK,
+    const Rd& uK,
+    const double EK,
+    const double gammaK,
+    const double cK,
+    const double pK,
+*/
+                              const Rd& normal,
+                              const bool check = false) const
+  {
+    Assert((l2Norm(normal) - 1) < 1e-12);
+
+    // const double& rho    = RoeState.rho;
+    const Rd& u_mean     = RoeState.U;
+    const double H_mean  = RoeState.H;
+    const double& cspeed = RoeState.c;
+    const double& gamma  = RoeState.gamma;
+    const double& uscaln = dot(u_mean, normal);
+    const double& u2     = dot(u_mean, u_mean);
+    // const R NrjCin=.5*unorm2;
+
+    const double c2  = cspeed * cspeed;
+    const double gm1 = gamma - 1;
+    const double K   = c2 + gm1 * (dot(u_mean, u_mean) - H_mean);
+
+    // Le jacobien est lineaire par rapport a la normale
+    // Ref : PAr ex. le papier de D. Chauvheid sur la tension de surface
+    Rpxp Jacobian;
+
+    Rdxd CentralT = identity;
+    CentralT *= uscaln;
+    CentralT += tensorProduct(u_mean, normal) - gm1 * tensorProduct(normal, u_mean);
+
+    Jacobian(0, 0) = 0.;
+    for (size_t dim = 0; dim < Dimension; ++dim) {
+      Jacobian(0, dim + 1) = normal[dim];
+    }
+    Jacobian(0, Dimension + 1) = 0;
+
+    for (size_t dim = 0; dim < Dimension; ++dim) {
+      Jacobian(dim + 1, 0) = K * normal[dim] - uscaln * u_mean[dim];
+      for (size_t dim2 = 0; dim2 < Dimension; ++dim2)
+        Jacobian(dim + 1, dim2 + 1) = CentralT(dim, dim2);
+      Jacobian(dim + 1, Dimension + 1) = gm1 * normal[dim];
+    }
+
+    Jacobian(Dimension + 1, 0) = (K - H_mean) * uscaln;
+    for (size_t dim = 0; dim < Dimension; ++dim)
+      Jacobian(Dimension + 1, dim + 1) = (H_mean * normal[dim] - gm1 * uscaln * u_mean[dim]);
+    Jacobian(Dimension + 1, Dimension + 1) = gamma * uscaln;
+
+    // Valeur propres..
+    Rp EigenValues;
+
+    EigenValues[0] = uscaln - cspeed;
+    for (size_t dim = 0; dim < Dimension; ++dim)   // vp multiplicite d
+      EigenValues[1 + dim] = uscaln;
+    EigenValues[1 + Dimension] = uscaln + cspeed;
+
+    const double maxabsVpNormal = std::max(fabs(EigenValues[0]), fabs(EigenValues[1 + Dimension]));
+
+    // Vecteur propres a droite et gauche
+
+    // hyper plan ortho à la normale
+    std::vector<Rd> ortho(Dimension - 1);
+    if constexpr (Dimension == 2) {
+      ortho[0] = Rd{normal[1], -normal[0]};   // aussi de norme 1
+    } else {
+      const double a = normal[0];
+      const double b = normal[1];
+      const double c = normal[2];
+
+      if ((a == b) and (b == c)) {
+        static double invsqrt2 = 1. / sqrt(2.);
+        static double invsqrt6 = 1. / sqrt(6.);
+
+        ortho[0] = Rd{invsqrt2, -invsqrt2, 0};
+        ortho[1] = Rd{invsqrt6, invsqrt6, -2 * invsqrt6};   // au signe pres
+
+      } else {
+        ortho[0] = Rd{b - c, -(a - c), a - b};
+        ortho[0] *= 1. / l2Norm(ortho[0]);
+        ortho[1] = Rd{a * (b + c) - b * b - c * c, b * (a + c) - a * a - c * c, c * (a + b) - a * a - b * b};
+        ortho[1] *= 1. / l2Norm(ortho[1]);   // au signe pres
+      }
+    }
+
+    Rpxp RightTligne;
+    RightTligne(0, 0) = 1;
+    for (size_t dim = 1; dim < Dimension + 1; ++dim)
+      RightTligne(0, dim) = u_mean[dim - 1] - cspeed * normal[dim - 1];
+    RightTligne(0, Dimension + 1) = H_mean - uscaln * cspeed;
+
+    RightTligne(1, 0) = 1;
+    for (size_t dim = 1; dim < Dimension + 1; ++dim)
+      RightTligne(1, dim) = u_mean[dim - 1];
+    RightTligne(1, Dimension + 1) = H_mean - c2 / gm1;
+
+    for (size_t dim = 1; dim < Dimension; ++dim) {
+      RightTligne(dim + 1, 0) = 0.;
+      for (size_t dim2 = 1; dim2 < Dimension + 1; ++dim2) {
+        RightTligne(dim + 1, dim2) = ortho[dim - 1][dim2 - 1];
+      }
+      RightTligne(dim + 1, Dimension + 1) = dot(u_mean, ortho[dim - 1]);
+    }
+
+    RightTligne(Dimension + 1, 0) = 1;
+    for (size_t dim = 1; dim < Dimension + 1; ++dim) {
+      RightTligne(Dimension + 1, dim) = u_mean[dim - 1] + cspeed * normal[dim - 1];
+    }
+    RightTligne(Dimension + 1, Dimension + 1) = H_mean + uscaln * cspeed;
+
+    Rpxp Right = transpose(RightTligne);
+
+    const double invc2 = 1. / c2;
+
+    Rpxp Left;   //(zero);
+    Left(0, 0) = .5 * invc2 * (K + uscaln * cspeed);
+    for (size_t dim = 1; dim < Dimension + 1; ++dim)
+      Left(0, dim) = .5 * invc2 * (-gm1 * u_mean[dim - 1] - cspeed * normal[dim - 1]);
+    Left(0, Dimension + 1) = .5 * invc2 * gm1;
+
+    Left(1, 0) = gm1 * invc2 * (H_mean - u2);
+    for (size_t dim = 1; dim < Dimension + 1; ++dim)
+      Left(1, dim) = gm1 * invc2 * u_mean[dim - 1];
+    Left(1, 1 + Dimension) = -gm1 * invc2;
+
+    for (size_t dim = 1; dim < Dimension; ++dim) {
+      Left(1 + dim, 0) = -dot(u_mean, ortho[dim - 1]);
+      for (size_t dim2 = 0; dim2 < Dimension; ++dim2)
+        Left(1 + dim, dim2 + 1) = ortho[dim - 1][dim2];
+      Left(1 + dim, 1 + Dimension) = 0;
+    }
+
+    Left(1 + Dimension, 0) = .5 * invc2 * (K - uscaln * cspeed);
+    for (size_t dim = 1; dim < Dimension + 1; ++dim)
+      Left(1 + Dimension, dim) = .5 * invc2 * (-gm1 * u_mean[dim - 1] + cspeed * normal[dim - 1]);
+    Left(1 + Dimension, Dimension + 1) = .5 * invc2 * gm1;
+
+    // check
+    // Rpxp prod = Right * Left;
+    // std::cout << prod << "\n";
+
+    Rpxp EigenMatAbs(identity);
+    for (size_t dim = 0; dim < Dimension + 2; ++dim) {
+      EigenMatAbs(dim, dim) *= fabs(EigenValues[dim]);
+    }
+    Rpxp ProdLefAbs  = Right * EigenMatAbs;
+    Rpxp AbsJacobian = ProdLefAbs * Left;
+    double maxErr    = 0;
+
+    if (check) {
+      Rpxp EigenAbs(identity);
+      for (size_t dim = 0; dim < Dimension + 2; ++dim) {
+        EigenAbs(dim, dim) *= EigenValues[dim];
+      }
+      Rpxp ProdLeft  = Right * EigenAbs;
+      Rpxp JacobianR = ProdLeft * Left;
+
+      // Rpxp id = identity;
+      maxErr = 0;
+      // double maxErrLeftRightId = 0;
+      for (size_t i = 0; i < Dimension + 2; ++i)
+        for (size_t j = 0; j < Dimension + 2; ++j) {
+          maxErr = std::max(maxErr, fabs(Jacobian(i, j) - JacobianR(i, j)));
+          // maxErrLeftRightId = std::max(maxErrLeftRightId, fabs(id(i, j) - prod(i, j)));
+        }
+      // std::clog << " maxErr " << maxErr << " maxErrUnit " << maxErrLeftRightId;
+      // throw NormalError("STOP");
+    }
+
+    return JacobianInformations(Jacobian, Left, Right, EigenValues, AbsJacobian, maxabsVpNormal, maxErr);
+  }
+
+  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()};
+    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()};
+
+    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();
+
+    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 "signe" 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 Gr, Gf
+
+    NodeValue<Rpxd> Gr{p_mesh->connectivity()};
+    Gr.fill(zero);
+    EdgeValue<Rpxd> Ge{p_mesh->connectivity()};
+    Ge.fill(zero);
+    FaceValue<Rpxd> Gf{p_mesh->connectivity()};
+    Gf.fill(zero);
+
+    NodeValue<double> MaxErrNode{p_mesh->connectivity()};
+    MaxErrNode.fill(0);
+    FaceValue<double> MaxErrFace{p_mesh->connectivity()};
+    MaxErrFace.fill(0);
+    EdgeValue<double> MaxErrEdge{p_mesh->connectivity()};
+    MaxErrEdge.fill(0);
+
+    // NodeValuePerCell<Rpxp> Lambda_rj{p_mesh->connectivity()};
+    CellValuePerNode<Rpxp> Lambda_rj{p_mesh->connectivity()};
+    Lambda_rj.fill(zero);
+    auto is_boundary_node = p_mesh->connectivity().isBoundaryNode();
+
+    parallel_for(
+      p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        const auto& node_to_cell                   = node_to_cell_matrix[r];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+        if (1 == 0) {
+          // if (is_boundary_node[r]) {
+          for (size_t l = 0; l < node_to_cell.size(); ++l) {
+            const CellId& j = node_to_cell[l];
+            const size_t R  = node_local_number_in_its_cells[l];
+
+            // const Rd& Cjr_loc         = Cjr(j, R);
+            // const double norm_Cjr_loc = l2Norm(Cjr_loc);
+            // const Rd& Njr_loc         = 1. / norm_Cjr_loc * Cjr_loc;
+
+            const double rho_bd = StateAtNode[j][R][0];
+            Rd U_bd;
+            for (size_t dim = 0; dim < Dimension; ++dim)
+              U_bd[dim] = StateAtNode[j][R][dim + 1] / rho_bd;
+            const double E_bd = StateAtNode[j][R][Dimension + 1] / rho_bd;
+            const double p_bd = pression(rho_bd, E_bd - .5 * dot(U_bd, U_bd), gamma);
+            AverageStateStructData RoeS(
+              RoeAverageState(rho_bd, U_bd, E_bd, gamma, p_bd, rho_bd, U_bd, E_bd, gamma, p_bd));
+            // MeanAverageState(rho_bd, U_bd, E_bd, gamma, p_bd, rho_bd, U_bd, E_bd, gamma, p_bd));
+
+            const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, njr(j, R));
+            Lambda_rj[r][l] = SignofASquareDiagonalisableMatrix(JacInfoK.RightEigenVectors, JacInfoK.EigenValues,
+                                                                JacInfoK.LeftEigenVectors);
+          }
+        } else {
+          Rpxp somme_Lambda = zero;
+
+          for (size_t l = 0; l < node_to_cell.size(); ++l) {
+            const CellId& j = node_to_cell[l];
+            const size_t R  = node_local_number_in_its_cells[l];
+
+            Rpxp Lambda_risum = zero;
+            //  for (size_t k = 0; k < l; ++k) {
+            for (size_t k = 0; k < node_to_cell.size(); ++k) {
+              if (k == l)
+                continue;
+              const CellId K = node_to_cell[k];
+
+              // Etat de Roe entre j et K
+              // Puis Eval de la Jacobienne dans la dir. njr
+              // Ou bien autre Moyenne
+              AverageStateStructData RoeS(
+                RoeAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K]));
+              // MeanAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K]));
+
+              const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, njr(j, R));
+              Lambda_risum += SignofASquareDiagonalisableMatrix(JacInfoK.RightEigenVectors, JacInfoK.EigenValues,
+                                                                JacInfoK.LeftEigenVectors);
+            }
+
+            // if (node_to_cell.size() > 1)
+            Lambda_risum *= 1. / std::max(1., (node_to_cell.size() - 1.));
+            Lambda_rj[r][l] = Lambda_risum;
+
+            somme_Lambda += Lambda_risum;
+          }
+
+          // somme_Lambda *= (-0.5 / (node_to_cell.size() - 1.));
+          // for (size_t l = 0; l < node_to_cell.size(); ++l)   // boucle sur les mailles voisines
+          // {
+          //   Lambda_rj[r][l] *= (node_to_cell.size() * 0.5 / (node_to_cell.size() - 1.));
+          //   Lambda_rj[r][l] += somme_Lambda;
+          // }   // fin boucle sur i
+
+          somme_Lambda *= -1. / (node_to_cell.size());
+          for (size_t l = 0; l < node_to_cell.size(); ++l)   // boucle sur les mailles voisines
+          {
+            Lambda_rj[r][l] += somme_Lambda;
+          }
+        }
+      });
+
+    synchronize(Lambda_rj);
+
+    if (checkLocalConservation) {
+      NodeValue<double> MaxFrobNormLambda{p_mesh->connectivity()};
+      MaxFrobNormLambda.fill(0.);
+      parallel_for(
+        p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+          if (not is_boundary_node[r]) {
+            const auto& node_to_cell = node_to_cell_matrix[r];
+            Rpxp somme_Lambda        = zero;
+            for (size_t l = 0; l < node_to_cell.size(); l++)   // boucle sur les mailles voisines
+              somme_Lambda += Lambda_rj[r][l];
+
+            const double FrobNorm = frobeniusNorm(somme_Lambda);
+            MaxFrobNormLambda[r]  = FrobNorm;
+            if (FrobNorm > 1e-9) {
+              std::cout << " Somme matrice noeud " << somme_Lambda << "\n";
+              throw NormalError("Problem of TLBL at node");
+            }
+          }
+        });
+      std::cout << " Max Frobenius Norm (node) " << max(MaxFrobNormLambda) << "\n";
+    }
+
+    // Ici on met jf ..
+    FaceValuePerCell<Rpxp> Lambda_jf{p_mesh->connectivity()};
+    Lambda_jf.fill(zero);
+    auto is_boundary_face = p_mesh->connectivity().isBoundaryFace();
+
+    parallel_for(
+      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_to_face = cell_to_face_matrix[j];
+        // const auto& cell_local_number_in_its_faces = cell_local_number_in_its_faces.itemArray(j);
+
+        for (size_t l = 0; l < cell_to_face.size(); ++l) {
+          const FaceId& face = cell_to_face[l];
+          if (is_boundary_face[face]) {
+            const double rho_bd = StateAtFace[j][l][0];
+            Rd U_bd;
+            for (size_t dim = 0; dim < Dimension; ++dim)
+              U_bd[dim] = StateAtFace[j][l][dim + 1] / rho_bd;
+            const double E_bd = StateAtFace[j][l][Dimension + 1] / rho_bd;
+            const double p_bd = pression(rho_bd, E_bd - .5 * dot(U_bd, U_bd), gamma);
+            AverageStateStructData RoeS(
+              RoeAverageState(rho_bd, U_bd, E_bd, gamma, p_bd, rho_bd, U_bd, E_bd, gamma, p_bd));
+            // MeanAverageState(rho_bd, U_bd, E_bd, gamma, p_bd, rho_bd, U_bd, E_bd, gamma, p_bd));
+
+            const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, njf(j, l));
+
+            Lambda_jf[j][l] = SignofASquareDiagonalisableMatrix(JacInfoK.RightEigenVectors, JacInfoK.EigenValues,
+                                                                JacInfoK.LeftEigenVectors);
+            // continue;
+          } else {
+            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);
+            //  const double norm_Cjf_loc = l2Norm(Cjf_loc);
+            // const Rd& Njf_loc = njf(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];
+              }
+            }
+
+            AverageStateStructData RoeS(
+              RoeAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K]));
+            // MeanAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K]));
+
+            const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, njf(j, l));
+
+            Lambda_jf[j][l] = SignofASquareDiagonalisableMatrix(JacInfoK.RightEigenVectors, JacInfoK.EigenValues,
+                                                                JacInfoK.LeftEigenVectors);
+          }
+        }
+      });
+    synchronize(Lambda_jf);
+
+    //
+    // Like node .. node / cells
+    //
+    CellValuePerEdge<Rpxp> Lambda_ej{p_mesh->connectivity()};
+    Lambda_ej.fill(zero);
+    auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge();
+
+    const EdgeValuePerCell<const Rd> Cje = mesh_data.Cje();
+    const EdgeValuePerCell<const Rd> nje = mesh_data.nje();
+
+    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();
+
+      parallel_for(
+        p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId e) {
+          const auto& edge_to_cell                   = edge_to_cell_matrix[e];
+          const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(e);
+
+          Rpxp somme_Lambda = zero;
+
+          for (size_t l = 0; l < edge_to_cell.size(); ++l) {
+            const CellId& j = edge_to_cell[l];
+            const size_t R  = edge_local_number_in_its_cells[l];
+            // const Rd& Cje_loc          = Cje(j, R);
+            // const double& norm_Cje_loc = l2Norm(Cje_loc);
+            // const Rd& Nje_loc          = 1. / (norm_Cje_loc)*Cje_loc;
+
+            Rpxp Lambda_risum = zero;
+
+            for (size_t k = 0; k < edge_to_cell.size(); ++k) {
+              if (k == l)
+                continue;
+              const CellId K = edge_to_cell[k];
+              // const size_t s            = edge_local_number_in_its_cells[k];
+              // const Rd& Cks_loc         = Cje(K, s);
+              // const double norm_Cks_loc = l2Norm(Cks_loc);
+              // const Rd& Nks_loc         = 1. / norm_Cks_loc * Cks_loc;
+
+              // Etat de Roe entre j et K
+              // Puis Eval de la Jacobienne dans la dir. Nks_loc
+              // Ou bien autre Moyenne
+              AverageStateStructData RoeS(
+                RoeAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K]));
+              // MeanAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K]));
+
+              const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, nje(j, R));
+              Lambda_risum += SignofASquareDiagonalisableMatrix(JacInfoK.LeftEigenVectors, JacInfoK.EigenValues,
+                                                                JacInfoK.RightEigenVectors);
+            }
+            // if (edge_to_cell.size() > 1)
+            Lambda_risum *= 1. / std::max(1., (edge_to_cell.size() - 1.));
+            Lambda_ej[e][l] = Lambda_risum;
+
+            somme_Lambda += Lambda_risum;
+          }
+          // somme_Lambda *= (-0.5 / (node_to_cell.size() - 1.));
+          // for (size_t l = 0; l < node_to_cell.size(); ++l)   // boucle sur les mailles voisines
+          // {
+          //   Lambda_rj[r][l] *= (node_to_cell.size() * 0.5 / (node_to_cell.size() - 1.));
+          //   Lambda_rj[r][l] += somme_Lambda;
+          // }   // fin boucle sur i
+
+          somme_Lambda *= -1. / (edge_to_cell.size());
+          for (size_t l = 0; l < edge_to_cell.size(); l++)   // boucle sur les mailles voisines
+          {
+            Lambda_ej[e][l] += somme_Lambda;
+          }
+        });
+      synchronize(Lambda_ej);
+    }
+    if (checkLocalConservation) {
+      const auto& edge_to_cell_matrix = p_mesh->connectivity().edgeToCellMatrix();
+
+      EdgeValue<double> MaxFrobNormLambda{p_mesh->connectivity()};
+      MaxFrobNormLambda.fill(0.);
+      // auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge();
+      parallel_for(
+        p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId e) {
+          if (not is_boundary_edge[e]) {
+            const auto& edge_to_cell = edge_to_cell_matrix[e];
+            Rpxp somme_Lambda        = zero;
+            for (size_t l = 0; l < edge_to_cell.size(); l++)   // boucle sur les mailles voisines
+              somme_Lambda += Lambda_ej[e][l];
+            const double FrobNorm = frobeniusNorm(somme_Lambda);
+            MaxFrobNormLambda[e]  = FrobNorm;
+
+            if (FrobNorm > 1e-9) {
+              std::cout << " Somme matrice edge " << somme_Lambda << "\n";
+              throw NormalError("Problem of TLBL at edge");
+            }
+          }
+        });
+      std::cout << " Max Frobenius Norm (edge) " << max(MaxFrobNormLambda) << "\n";
+    }
+
+    NodeValue<int> nbChangingSignNode{p_mesh->connectivity()};
+    nbChangingSignNode.fill(0);
+
+    parallel_for(
+      p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node) {
+        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);
+        // if (is_boundary_node[node]) {
+        if (1 == 0) {
+          Rpxd FluxMoy;
+          FluxMoy = zero;
+          Rp Umoy = zero;
+          for (size_t l = 0; l < node_to_cell.size(); ++l) {
+            const CellId& j = node_to_cell[l];
+            const size_t R  = node_local_number_in_its_cells[l];
+
+            const Rpxp Lambda_rl        = Lambda_rj[node][l];
+            Rpxp MatriceIDPlusLambda_rj = identity;
+            MatriceIDPlusLambda_rj += Lambda_rl;
+
+            Umoy += StateAtNode[j][R];
+
+            Rpxd Flux;
+            for (size_t d = 0; d < Dimension; ++d) {
+              Flux(0, d) = Flux_rhoAtCellNode[j][R][d];   // dot(Flux_rhoAtCellNode[j][R], Cjr_loc);
+              // dot(Flux_rho[K], Cjr_loc);
+              for (size_t dd = 0; dd < Dimension; ++dd)
+                Flux(dd + 1, d) = Flux_qtmvtAtCellNode[j][R](dd, d);   //  u_Cjr[d];
+
+              Flux(Dimension + 1, d) = Flux_totnrjAtCellNode[j][R][d];   // dot(Flux_totnrjAtCellNode[j][R], Cjr_loc);
+                                                                         // dot(Flux_totnrj[K], Cjr_loc);
+            }
+            FluxMoy += MatriceIDPlusLambda_rj * Flux;
+          }
+          // Gr[node] = (1. / node_to_cell.size()) * FluxMoy;   // MatriceIDPlusLambda_rj * Flux;
+          Umoy *= 1. / node_to_cell.size();
+          const double rhonode = Umoy[0];
+          Rd rhoUnode;
+          for (size_t d = 0; d < Dimension; ++d)
+            rhoUnode[d] = Umoy[d + 1];
+          const double rhoEnode    = Umoy[Dimension + 1];
+          const double Enode       = rhoEnode / rhonode;
+          const Rd Unode           = 1. / rhonode * rhoUnode;
+          const double epsilonnode = Enode - .5 * dot(Unode, 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;
+
+          Rpxd Flux_Umoy;
+          for (size_t d = 0; d < Dimension; ++d) {
+            Flux_Umoy(0, d) = rhoUnode[d];
+            for (size_t dd = 0; dd < Dimension; ++dd)
+              Flux_Umoy(dd + 1, d) = rhoUtensUPlusPidnode(dd, d);
+            Flux_Umoy(Dimension + 1, d) = rhoEnodePlusP * Unode[d];
+          }
+          Gr[node] = Flux_Umoy;   //(1. / node_to_cell.size()) * FluxMoy;   // MatriceIDPlusLambda_rj * Flux;
+        } else {
+          for (size_t l = 0; l < node_to_cell.size(); ++l) {
+            const CellId& j = node_to_cell[l];
+            const size_t R  = node_local_number_in_its_cells[l];
+            // const Rd& Cjr_loc    = Cjr(j, R);
+            const Rpxp Lambda_rl = Lambda_rj[node][l];
+
+            Rpxp MatriceIDPlusLambda_rj = identity;
+            MatriceIDPlusLambda_rj += Lambda_rl;
+
+            Rpxd Flux;
+            for (size_t d = 0; d < Dimension; ++d) {
+              Flux(0, d) = Flux_rhoAtCellNode[j][R][d];
+              for (size_t dd = 0; dd < Dimension; ++dd)
+                Flux(dd + 1, d) = Flux_qtmvtAtCellNode[j][R](dd, d);
+              Flux(Dimension + 1, d) = Flux_totnrjAtCellNode[j][R][d];
+            }
+            Gr[node] += MatriceIDPlusLambda_rj * Flux;
+          }
+          Gr[node] *= 1. / (node_to_cell.size());
+        }
+      });
+    synchronize(Gr);
+
+    if (checkLocalConservation) {
+      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 += Gr[l] * Cjr(K, R);
+              maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gr[l] * Cjr(K, R)));
+            }
+            // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr));
+            MaxErrorConservationNode[l] = l2Norm(SumGjr) / maxGjrAbs;
+          }
+        });
+      std::cout << " Max Error Node " << max(MaxErrorConservationNode) << " Max Error RoeMatrice  " << max(MaxErrNode)
+                << " nb Chg Sign " << sum(nbChangingSignNode) << "\n";
+    }
+
+    FaceValue<int> nbChangingSignFace{p_mesh->connectivity()};
+    nbChangingSignFace.fill(0);
+    const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+
+    parallel_for(
+      p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId f) {
+        const auto& face_to_cell                   = face_to_cell_matrix[f];
+        const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(f);
+
+        // if (1 == 0) {   //}
+        if (is_boundary_face[f]) {
+          const CellId& j = face_to_cell[0];
+          const size_t R  = face_local_number_in_its_cells[0];
+
+          Rpxd FluxJ;
+          for (size_t d = 0; d < Dimension; ++d) {
+            FluxJ(0, d) = Flux_rhoAtCellFace[j][R][d];
+            for (size_t dd = 0; dd < Dimension; ++dd)
+              FluxJ(dd + 1, d) = Flux_qtmvtAtCellFace[j][R](dd, d);
+            FluxJ(Dimension + 1, d) = Flux_totnrjAtCellFace[j][R][d];
+          }
+
+          Gf[f] = FluxJ;   // MatriceIdplusLambda * FluxJ;
+        } else {
+          for (size_t l = 0; l < face_to_cell.size(); ++l) {
+            const CellId& j = face_to_cell[l];
+            const size_t R  = face_local_number_in_its_cells[l];
+
+            const Rpxp& Lambda_fj = Lambda_jf[j][R];
+
+            Rpxp MatriceIDPlusLambda_fj = identity;
+            MatriceIDPlusLambda_fj += Lambda_fj;
+
+            // const Rp& statediff = StateAtNode[j][l] - StateAtNode[K][R];   // State[j] - State[K];
+            // const Rp& diff      = ViscosityMatrixJK * statediff;
+
+            Rpxd Flux;
+            for (size_t d = 0; d < Dimension; ++d) {
+              Flux(0, d) = Flux_rhoAtCellFace[j][R][d];
+              for (size_t dd = 0; dd < Dimension; ++dd)
+                Flux(dd + 1, d) = Flux_qtmvtAtCellFace[j][R](dd, d);
+              Flux(Dimension + 1, d) = Flux_totnrjAtCellFace[j][R][d];
+            }
+
+            Gf[f] += MatriceIDPlusLambda_fj * Flux;
+          }
+          Gf[f] *= 1. / (face_to_cell.size());
+        }
+      });
+    synchronize(Gf);
+
+    if (checkLocalConservation) {
+      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 += Gf[l] * Cjf(K, R);
+              maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gf[l] * Cjf(K, R)));
+            }
+            MaxErrorConservationFace[l] = l2Norm(SumGjf) / maxGjrAbs;
+            // MaxErrorConservationFace   = std::max(MaxErrorConservationFace, l2Norm(SumGjf));
+          }
+        });
+      std::cout << " Max Error Face " << max(MaxErrorConservationFace) << " Max Error RoeMatrice  " << max(MaxErrFace)
+                << " nb Chg Sign " << sum(nbChangingSignFace) << "\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();
+
+      EdgeValue<int> nbChangingSignEdge(p_mesh->connectivity());
+      nbChangingSignEdge.fill(0);
+
+      parallel_for(
+        p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId e) {
+          // Edge
+
+          const auto& edge_to_cell                   = edge_to_cell_matrix[e];
+          const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(e);
+
+          for (size_t l = 0; l < edge_to_cell.size(); ++l) {
+            const CellId& j = edge_to_cell[l];
+            const size_t R  = edge_local_number_in_its_cells[l];
+            // const Rd& Cjr_loc    = Cjr(j, R);
+            const Rpxp Lambda_el = Lambda_ej[e][l];
+
+            Rpxp MatriceIDPlusLambda_ej = identity;
+            MatriceIDPlusLambda_ej += Lambda_el;
+
+            Rpxd Flux;
+            for (size_t d = 0; d < Dimension; ++d) {
+              Flux(0, d) = Flux_rhoAtCellEdge[j][R][d];
+              for (size_t dd = 0; dd < Dimension; ++dd)
+                Flux(dd + 1, d) = Flux_qtmvtAtCellEdge[j][R](dd, d);
+              Flux(Dimension + 1, d) = Flux_totnrjAtCellEdge[j][R][d];
+            }
+            Ge[e] += MatriceIDPlusLambda_ej * Flux;
+          }
+          Ge[e] *= 1. / (edge_to_cell.size());
+          /*
+          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);
+
+            // Moyenne de 2 etats
+            // 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;
+
+              // La moyenne de Roe  entre les etats jk
+              AverageStateStructData RoeS(
+                RoeAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K]));
+              // MeanAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K],
+              gamma, p_n[K]));
+            //  Viscosity j k
+            const double nrmCke                  = l2Norm(Cke_loc);
+            const Rd& normal                     = 1. / nrmCke * Cke_loc;
+            const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, normal);
+
+            const double nrmCje                  = l2Norm(Cje_loc);
+            const Rd& normalJ                    = 1. / nrmCje * Cje_loc;
+            const JacobianInformations& JacInfoJ = JacobianFluxAlongUnitNormal(RoeS, normalJ);
+
+            // Matrice de Viscosite : Valeur Abs de la jacobienne en l'etat de Roe
+            Rpxp ViscosityMatrixJK = .5 * (nrmCje * JacInfoJ.AbsJacobian + nrmCke * JacInfoK.AbsJacobian);
+
+            // Test Rajout Viscosité type Rusanov v2 si chgt signe vp ..
+
+            bool anySignChgJ = EvaluateChangingSignVpAlongNormal(   // rho_n[j],
+              u_n[j],                                               // E_n[j], gamma,
+              c_n[j],                                               // p_n[j], rho_n[K],
+              u_n[K],
+              // E_n[K], gamma,
+              c_n[K],   // p_n[K],
+              normalJ);
+            bool anySignChgK = EvaluateChangingSignVpAlongNormal(   // rho_n[j],
+              u_n[j],                                               // E_n[j], gamma,
+              c_n[j],                                               // p_n[j], rho_n[K],
+              u_n[K],
+              // E_n[K], gamma,
+              c_n[K],   // p_n[K],
+              normal);
+
+            bool anySignDif = false;
+            if (anySignChgJ or anySignChgK)
+              anySignDif = true;   // false;   // true;
+
+            if (anySignDif) {
+              Rpxp AddedViscosity(identity);
+              AddedViscosity *= std::max(JacInfoJ.maxabsVpNormal * nrmCje, JacInfoK.maxabsVpNormal * nrmCke);
+              ViscosityMatrixJK += AddedViscosity;
+              nbChangingSignEdge[edge] = 1;
+            }
+
+            MaxErrEdge[edge] = std::max(MaxErrEdge[edge], JacInfoK.MaxErrorProd);
+
+            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(Ge);
+
+      if (checkLocalConservation) {
+        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 += Ge[l] * Cje(K, R);
+                maxGjrAbs = std::max(maxGjrAbs, l2Norm(Ge[l] * Cje(K, R)));
+              }
+              // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje));
+              MaxErrorConservationEdge[l] = l2Norm(SumGje) / maxGjrAbs;
+            }
+          });
+        std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << " Max Error RoeMatrice  " << max(MaxErrEdge)
+                  << " nb Chg Sign " << sum(nbChangingSignEdge) << "\n";
+      }
+
+    }   // dim 3
+
+    // Pour les assemblages
+    double theta = 2. / 3.;   // std::atan(1);   // pi / 4.;   // 2. / 3.;   // 2. / 3.;   // 2. / 3.;   //.5;sth
+    double eta   = 1. / 6.;   //.2;
+    if constexpr (Dimension == 2) {
+      eta = 0;
+      // theta = 1;   // pour schema aux aretes
+      // theta=0; //pour schema aux noeuds
+    }
+
+    // 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) {
+          const NodeId& node = cell_to_node[l];
+          SumFluxesNode += Gr[node] * Cjr(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) {
+          const EdgeId& edge = cell_to_edge[l];
+
+          SumFluxesEdge += Ge[edge] * Cje(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) {
+          const FaceId& face = cell_to_face[l];
+
+          SumFluxesFace += Gf[face] * Cjf(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));
+  }
+
+  RoeFluxFormEulerianCompositeSolver_v2()  = default;
+  ~RoeFluxFormEulerianCompositeSolver_v2() = default;
+};
+
+template <MeshConcept MeshType>
+class RoeFluxFormEulerianCompositeSolver_v2<MeshType>::WallBoundaryCondition
+{
+};
+
+template <>
+class RoeFluxFormEulerianCompositeSolver_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 RoeFluxFormEulerianCompositeSolver_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 RoeFluxFormEulerianCompositeSolver_v2<MeshType>::NeumannflatBoundaryCondition
+{
+};
+template <>
+class RoeFluxFormEulerianCompositeSolver_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 RoeFluxFormEulerianCompositeSolver_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 RoeFluxFormEulerianCompositeSolver_v2<MeshType>::SymmetryBoundaryCondition
+{
+};
+
+template <>
+class RoeFluxFormEulerianCompositeSolver_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 RoeFluxFormEulerianCompositeSolver_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 RoeFluxFormEulerianCompositeSolver_v2<MeshType>::InflowListBoundaryCondition
+{
+};
+
+template <>
+class RoeFluxFormEulerianCompositeSolver_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 RoeFluxFormEulerianCompositeSolver_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 RoeFluxFormEulerianCompositeSolver_v2<MeshType>::OutflowBoundaryCondition
+{
+};
+
+template <>
+class RoeFluxFormEulerianCompositeSolver_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 RoeFluxFormEulerianCompositeSolver_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>>
+roeFluxFormEulerianCompositeSolver_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("RoeFluxFormFormEulerianCompositeSolver v2 is not available in 1D");
+        } else {
+          if constexpr (is_polygonal_mesh_v<MeshType>) {
+            return RoeFluxFormEulerianCompositeSolver_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>>(),   // degree,
+                     bc_descriptor_list, dt, check);
+          } else {
+            throw NormalError("RoeFluxFormEulerianCompositeSolver v2 is only defined on polygonal meshes");
+          }
+        }
+      },
+    mesh_v->variant());
+}
diff --git a/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.hpp b/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.hpp
new file mode 100644
index 000000000..e7f6519c6
--- /dev/null
+++ b/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.hpp
@@ -0,0 +1,89 @@
+#ifndef ROE_FLUX_FORM_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+#define ROE_FLUX_FORM_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+
+#include <mesh/MeshTraits.hpp>
+
+#include "JacobianAndStructuralInfoForSystemsofEquations.hpp"
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+
+#include <memory>
+#include <tuple>
+#include <vector>
+
+inline double
+signe(const double a)
+{
+  if (a < 0)
+    return -1.;
+  if (a > 0)
+    return 1;
+
+  return 0;
+}
+class Rp;
+class Rpxp;
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+roeFluxFormEulerianCompositeSolver_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);
+/*
+class DiscreteFunctionVariant;
+class IBoundaryConditionDescriptor;
+class MeshVariant;
+template <size_t Dimension>
+class Mesh;
+
+template <size_t Dimension, int number_of_unknowns>
+class RoeFluxFormEulerianCompositeSolver_v2
+{
+ private:
+  static constexpr int p = number_of_unknowns + (Dimension - 1);
+  using Rp               = TinyVector<p>;   // number_of_unknowns + (Dimension - 1)>;
+
+  using Rpxp = TinyMatrix<p, p, double>;
+  // private:
+ public:
+  JacobianStructuralInfoSystemsOfEquations<Dimension, number_of_unknowns> JacobianInfos;
+
+  // KN<std::vector<std::vector<Rn> > > _sign_jacobian_matrix_dof_in_cell;
+  // KN<std::vector<std::vector<Rn> > > _vp_jacobian_matrix_dof_in_cell;
+  // //KN<std::vector<Rnm> > _jac_jacobian_matrix_dof_in_cell;
+  // KN<std::vector<std::vector<Rnm> > > _right_eigenv_matrix_dof_in_cell;
+  // KN<std::vector<std::vector<Rnm> > >  _left_eigenv_matrix_dof_in_cell;
+  // //KN<std::vector<Rnm> >  _rotation_matrix_dof_in_cell;
+  // KN<std::vector<R2nm> >  _MatriceU_dof_in_cell;
+
+  NodeValuePerCell<const Rp> _vp_jacobian_matrix_node_in_cell;
+  EdgeValuePerCell<const Rp> _vp_jacobian_matrix_edge_in_cell;
+  FaceValuePerCell<const Rp> _vp_jacobian_matrix_face_in_cell;
+
+  NodeValuePerCell<const Rpxp> _right_eigenv_matrix_node_in_cell;
+  EdgeValuePerCell<const Rpxp> _right_eigenv_matrix_edge_in_cell;
+  FaceValuePerCell<const Rpxp> _right_eigenv_matrix_face_in_cell;
+
+  NodeValuePerCell<const Rpxp> _left_eigenv_matrix_node_in_cell;
+  EdgeValuePerCell<const Rpxp> _left_eigenv_matrix_edge_in_cell;
+  FaceValuePerCell<const Rpxp> _left_eigenv_matrix_face_in_cell;
+
+  NodeValuePerCell<const Rpxp> _MatriceU_node_in_cell;
+  EdgeValuePerCell<const Rpxp> _MatriceU_edge_in_cell;
+  FaceValuePerCell<const Rpxp> _MatriceU_face_in_cell;
+
+ public:
+  RoeFluxFormEulerianCompositeSolver_v2() {}
+
+  ~RoeFluxFormEulerianCompositeSolver_v2() {}
+};
+*/
+#endif
-- 
GitLab