diff --git a/src/language/modules/KineticSchemeModule.cpp b/src/language/modules/KineticSchemeModule.cpp
index bfe0d148ee68cba722cb8e0c74b66e0069efb682..9c16e61229a9efc533ca37ff26b560d96ecb1f67 100644
--- a/src/language/modules/KineticSchemeModule.cpp
+++ b/src/language/modules/KineticSchemeModule.cpp
@@ -9,6 +9,7 @@
 #include <scheme/EulerKineticSolverAcousticLagrangeFV.hpp>
 #include <scheme/EulerKineticSolverFirstOrderFV.hpp>
 #include <scheme/EulerKineticSolverLagrangeFV.hpp>
+#include <scheme/EulerKineticSolverLagrangeMultiD.hpp>
 #include <scheme/EulerKineticSolverMeanFluxMood.hpp>
 #include <scheme/EulerKineticSolverMoodFD.hpp>
 #include <scheme/EulerKineticSolverMoodFV.hpp>
@@ -799,6 +800,66 @@ KineticSchemeModule::KineticSchemeModule()
 
                                           ));
 
+  this->_addBuiltinFunction("euler_kinetic_lagrange_multiD",
+                            std::function(
+
+                              [](const double& dt, const std::vector<TinyVector<1>>& lambda_vector, const double& gamma,
+                                 const double& eps, const size_t& space_order, const size_t& time_order,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return eulerKineticSolverLagrangeMultiD(dt, lambda_vector, gamma, eps, space_order,
+                                                                        time_order, rho, u, E, p, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_lagrange_multiD",
+                            std::function(
+
+                              [](const double& dt, const std::vector<TinyVector<2>>& lambda_vector, const double& gamma,
+                                 const double& eps, const size_t& space_order, const size_t& time_order,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return eulerKineticSolverLagrangeMultiD(dt, lambda_vector, gamma, eps, space_order,
+                                                                        time_order, rho, u, E, p, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_lagrange_multiD",
+                            std::function(
+
+                              [](const double& dt, const std::vector<TinyVector<3>>& lambda_vector, const double& gamma,
+                                 const double& eps, const size_t& space_order, const size_t& time_order,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return eulerKineticSolverLagrangeMultiD(dt, lambda_vector, gamma, eps, space_order,
+                                                                        time_order, rho, u, E, p, bc_descriptor_list);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("euler_kinetic_acoustic_lagrange_FV",
                             std::function(
 
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 135500f6389fd1e8fcf6b6d982a42b294cc43a65..d563dcd9f95be14792a7c1a2453ee990c9b3f6ea 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -17,6 +17,7 @@ add_library(
   WavesEquationNonUniformKineticSolver.cpp
   WavesEquationNonUniformCellKineticSolver.cpp
   EulerKineticSolver.cpp
+  EulerKineticSolverLagrangeMultiD.cpp
   EulerKineticSolverMeanFluxMood.cpp
   EulerKineticSolverOneFluxMood.cpp
   EulerKineticSolverMoodFD.cpp
diff --git a/src/scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.cpp b/src/scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.cpp
index 88889c40afe4466273c4f34659dcf560875425fa..81ff7dc085add37a8bd76c9419414bfa9bdab2d7 100644
--- a/src/scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.cpp
+++ b/src/scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.cpp
@@ -239,10 +239,11 @@ class EulerKineticSolverAcoustic2LagrangeFVOrder2
       lambda = std::max(lambda, m_lambda_vector[i_wave]);
     }
 
-    for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
-      vec_A[cell_id][0] = lambda * lambda * u[cell_id][0];
-      vec_A[cell_id][1] = p[cell_id];
-    }
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        vec_A[cell_id][0] = lambda * lambda * u[cell_id][0];
+        vec_A[cell_id][1] = p[cell_id];
+      });
 
     return vec_A;
   }
@@ -261,13 +262,14 @@ class EulerKineticSolverAcoustic2LagrangeFVOrder2
     }
     inv_L2norm_lambda_vector = 1. / inv_L2norm_lambda_vector;
 
-    for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
-      for (size_t i = 0; i < nb_waves; ++i) {
-        Mp[cell_id][i] = (1. / nb_waves) * p[cell_id] + inv_L2norm_lambda_vector * A[cell_id][0] * m_lambda_vector[i];
-        Mu[cell_id][i] =
-          (1. / nb_waves) * u[cell_id][0] + inv_L2norm_lambda_vector * A[cell_id][1] * m_lambda_vector[i];
-      }
-    }
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        for (size_t i = 0; i < nb_waves; ++i) {
+          Mp[cell_id][i] = (1. / nb_waves) * p[cell_id] + inv_L2norm_lambda_vector * A[cell_id][0] * m_lambda_vector[i];
+          Mu[cell_id][i] =
+            (1. / nb_waves) * u[cell_id][0] + inv_L2norm_lambda_vector * A[cell_id][1] * m_lambda_vector[i];
+        }
+      });
     return std::make_tuple(DiscreteFunctionP0Vector<double>(Mp), DiscreteFunctionP0Vector<double>(Mu));
   }
 
@@ -280,31 +282,32 @@ class EulerKineticSolverAcoustic2LagrangeFVOrder2
     NodeArray<double> Flux_F{m_mesh.connectivity(), nb_waves};
     Flux_F.fill(0);
 
-    for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
-      for (size_t i = 0; i < nb_waves; ++i) {
-        const size_t N = m_mesh.numberOfCells();
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        for (size_t i = 0; i < nb_waves; ++i) {
+          const size_t N = m_mesh.numberOfCells();
 
-        const auto& cell_to_node   = cell_to_node_matrix[cell_id];
-        const NodeId left_node_id  = cell_to_node[0];
-        const NodeId right_node_id = cell_to_node[1];
+          const auto& cell_to_node   = cell_to_node_matrix[cell_id];
+          const NodeId left_node_id  = cell_to_node[0];
+          const NodeId right_node_id = cell_to_node[1];
 
-        const CellId prev_cell_id = (cell_id + N - 1) % N;
-        const CellId next_cell_id = (cell_id + 1) % N;
+          const CellId prev_cell_id = (cell_id + N - 1) % N;
+          const CellId next_cell_id = (cell_id + 1) % N;
 
-        if (m_lambda_vector[i] < 0) {
-          Flux_F[left_node_id][i]  = F[cell_id][i];
-          Flux_F[right_node_id][i] = F[next_cell_id][i];
+          if (m_lambda_vector[i] < 0) {
+            Flux_F[left_node_id][i]  = F[cell_id][i];
+            Flux_F[right_node_id][i] = F[next_cell_id][i];
 
-        } else if (m_lambda_vector[i] > 0) {
-          Flux_F[left_node_id][i]  = F[prev_cell_id][i];
-          Flux_F[right_node_id][i] = F[cell_id][i];
+          } else if (m_lambda_vector[i] > 0) {
+            Flux_F[left_node_id][i]  = F[prev_cell_id][i];
+            Flux_F[right_node_id][i] = F[cell_id][i];
 
-        } else {
-          Flux_F[right_node_id][i] = 0;
-          Flux_F[left_node_id][i]  = 0;
+          } else {
+            Flux_F[right_node_id][i] = 0;
+            Flux_F[left_node_id][i]  = 0;
+          }
         }
-      }
-    }
+      });
     return Flux_F;
   }
 
diff --git a/src/scheme/EulerKineticSolverLagrangeMultiD.cpp b/src/scheme/EulerKineticSolverLagrangeMultiD.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..70423dfc07916d6f20cec371604b3bd53e1b37f4
--- /dev/null
+++ b/src/scheme/EulerKineticSolverLagrangeMultiD.cpp
@@ -0,0 +1,795 @@
+#include <scheme/EulerKineticSolverLagrangeMultiD.hpp>
+
+#include <language/utils/EvaluateAtPoints.hpp>
+#include <language/utils/InterpolateItemArray.hpp>
+#include <language/utils/InterpolateItemValue.hpp>   //new
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatFaceBoundary.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>   //new
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <scheme/WallBoundaryConditionDescriptor.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <tuple>
+
+template <MeshConcept MeshType>
+class EulerKineticSolverLagrangeMultiD
+{
+ private:
+  constexpr static size_t Dimension = MeshType::Dimension;
+  using Rd                          = TinyVector<Dimension>;
+
+  class SymmetryBoundaryCondition;
+  class WallBoundaryCondition;
+  class InflowListBoundaryCondition;
+  class OutflowBoundaryCondition;
+
+  using BoundaryCondition = std::
+    variant<SymmetryBoundaryCondition, WallBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition>;
+  using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+  const double m_dt;
+  const std::vector<TinyVector<Dimension>> m_lambda_vector;
+  const double m_gamma;
+  const size_t m_spaceOrder;
+  const size_t m_timeOrder;
+  const CellValue<const double> m_rho;
+  const CellValue<const TinyVector<Dimension>> m_u;
+  const CellValue<const double> m_E;
+  const CellValue<const double> m_p;
+  const std::shared_ptr<const MeshType> p_mesh;
+  const MeshType& m_mesh;
+  const CellValue<const double> m_Vj = MeshDataManager::instance().getMeshData(m_mesh).Vj();
+  CellValue<const double> m_inv_Mj;
+  CellValue<const double> m_inv_Vj;
+  CellValue<const double> m_Mj;
+  const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
+  const NodeValue<const TinyVector<Dimension>> m_xr = m_mesh.xr();
+
+  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::symmetry: {
+        bc_list.emplace_back(
+          SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+                                    getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::wall: {
+        bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+                                                   getMeshFaceBoundary(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() != 1 + 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 " << 1 + Dimension;
+          throw NormalError(error_msg.str());
+        }
+        auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+        // CellValue<size_t> is_boundary_cell{p_mesh->connectivity()};
+
+        // is_boundary_cell.fill(0);
+
+        // auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+
+        auto node_list = node_boundary.nodeList();
+        // for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+        //   auto cell_list = node_to_cell_matrix[node_list[i_node]];
+        //   for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+        //     const CellId cell_id      = cell_list[i_cell];
+        //     is_boundary_cell[cell_id] = 1;
+        //   }
+        // }
+
+        // size_t nb_boundary_cells = sum(is_boundary_cell);
+
+        // Array<CellId> cell_list{nb_boundary_cells};
+
+        // size_t index = 0;
+        // for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+        //   if (is_boundary_cell[cell_id]) {
+        //     cell_list[index++] = cell_id;
+        //   }
+        // }
+        //        auto xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
+        auto xr = m_mesh.xr();
+
+        // Table<const double> cell_array_list =
+        //   InterpolateItemArray<double(Rd)>::template interpolate<ItemType::cell>(inflow_list_bc_descriptor
+        //                                                                            .functionSymbolIdList(),
+        //                                                                          xj, cell_list);
+        Table<const double> node_array_list =
+          InterpolateItemArray<double(Rd)>::template interpolate<ItemType::node>(inflow_list_bc_descriptor
+                                                                                   .functionSymbolIdList(),
+                                                                                 xr, node_list);
+        // bc_list.emplace_back(InflowListBoundaryCondition(cell_array_list, cell_list));
+        bc_list.emplace_back(InflowListBoundaryCondition(node_array_list, node_list));
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::outflow: {
+        bc_list.emplace_back(OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+                                                      getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+        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 acoustic solver";
+        throw NormalError(error_msg.str());
+      }
+    }
+
+    return bc_list;
+  }
+
+ public:
+  const CellArray<const double>
+  compute_M_p(const CellValue<const double>& p, const CellValue<const TinyVector<Dimension>>& u, const double& lambda)
+  {
+    const size_t nb_waves = m_lambda_vector.size();
+    CellArray<double> M{p_mesh->connectivity(), nb_waves};
+    TinyVector<Dimension> inv_S = zero;
+    for (size_t d = 0; d < Dimension; ++d) {
+      for (size_t i = 0; i < nb_waves; ++i) {
+        inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
+      }
+      inv_S[d] = 1. / inv_S[d];
+    }
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        for (size_t i = 0; i < nb_waves; ++i) {
+          M[cell_id][i] = (1. / nb_waves) * p[cell_id];
+
+          for (size_t d = 0; d < Dimension; ++d) {
+            M[cell_id][i] += inv_S[d] * m_lambda_vector[i][d] * (lambda * lambda * u[cell_id][d]);
+          }
+        }
+      });
+
+    return M;
+  }
+
+  const CellArray<const TinyVector<Dimension>>
+  compute_M_u(const CellValue<const TinyVector<Dimension>>& u, const CellValue<const double>& p)
+  {
+    const size_t nb_waves = m_lambda_vector.size();
+    CellArray<TinyVector<Dimension>> M{p_mesh->connectivity(), nb_waves};
+    TinyVector<Dimension> inv_S = zero;
+    TinyMatrix<Dimension> I     = identity;
+    for (size_t d = 0; d < Dimension; ++d) {
+      for (size_t i = 0; i < nb_waves; ++i) {
+        inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
+      }
+      inv_S[d] = 1. / inv_S[d];
+    }
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+          M[cell_id][i_wave] = 1. / nb_waves * u[cell_id];
+          for (size_t d1 = 0; d1 < Dimension; ++d1) {
+            for (size_t d2 = 0; d2 < Dimension; ++d2) {
+              M[cell_id][i_wave][d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * p[cell_id] * I(d2, d1);
+            }
+          }
+        }
+      });
+    return M;
+  }
+
+  template <typename T>
+  NodeArray<T>
+  compute_Flux(const CellArray<const T>& Fn, const NodeValue<const bool>& is_boundary_node)
+  {
+    const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+    const size_t nb_waves                                   = m_lambda_vector.size();
+    NodeArray<T> Fr(m_mesh.connectivity(), nb_waves);
+    const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
+    const auto& node_to_cell_matrix               = p_mesh->connectivity().nodeToCellMatrix();
+
+    if constexpr (is_tiny_vector_v<T>) {
+      Fr.fill(zero);
+    } else {
+      Fr.fill(0.);
+    }
+
+    parallel_for(
+      p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+        if (not is_boundary_node[node_id]) {
+          const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+          const auto& node_to_cell                  = node_to_cell_matrix[node_id];
+          for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+            double sum_li_njr = 0;
+
+            for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+              const CellId cell_id      = node_to_cell[i_cell];
+              const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+
+              double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node));
+              if (li_njr > 1e-14) {
+                Fr[node_id][i_wave] += li_njr * Fn[cell_id][i_wave];
+                sum_li_njr += li_njr;
+              }
+            }
+            if (sum_li_njr != 0) {
+              Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
+            }
+          }
+        }
+      });
+
+    return Fr;
+  }
+
+  const NodeValue<TinyVector<Dimension>>
+  compute_velocity_flux(NodeArray<double>& F, const double& lambda)
+  {
+    NodeValue<TinyVector<Dimension>> u{p_mesh->connectivity()};
+    const size_t nb_waves = m_lambda_vector.size();
+    u.fill(zero);
+    const double inv_lambda_squared = (1. / (lambda * lambda));
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+        for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+          u[node_id] += F[node_id][i_wave] * m_lambda_vector[i_wave];
+        }
+        u[node_id] = inv_lambda_squared * u[node_id];
+      });
+    return u;
+  }
+
+  const NodeValue<TinyVector<Dimension>>
+  compute_pressure_flux(NodeArray<TinyVector<Dimension>>& F)
+  {
+    NodeValue<TinyVector<Dimension>> p{p_mesh->connectivity()};
+    const size_t nb_waves = m_lambda_vector.size();
+    p.fill(zero);
+
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+        for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+          for (size_t dim = 0; dim < Dimension; ++dim) {
+            p[node_id][dim] += F[node_id][i_wave][dim] * m_lambda_vector[i_wave][dim];
+          }
+        }
+      });
+    return p;
+  }
+
+  const CellValue<const double>
+  compute_delta_velocity(const NodeValue<const TinyVector<Dimension>>& ur)
+  {
+    const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+    CellValue<double> delta{p_mesh->connectivity()};
+
+    const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+    delta.fill(0.);
+
+    parallel_for(
+      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        const auto& cell_to_node = cell_to_node_matrix[cell_id];
+
+        for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+          const NodeId node_id = cell_to_node[i_node];
+          delta[cell_id] += dot(Cjr(cell_id, i_node), ur[node_id]);
+        }
+      });
+    return delta;
+  }
+
+  const CellValue<const TinyVector<Dimension>>
+  compute_delta_pressure(const NodeValue<const TinyVector<Dimension>>& pr)
+  {
+    const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+    CellValue<TinyVector<Dimension>> delta{p_mesh->connectivity()};
+
+    const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+    delta.fill(zero);
+
+    parallel_for(
+      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        const auto& cell_to_node = cell_to_node_matrix[cell_id];
+
+        for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+          for (size_t dim = 0; dim < Dimension; ++dim) {
+            const NodeId node_id = cell_to_node[i_node];
+            delta[cell_id][dim] += Cjr(cell_id, i_node)[dim] * pr[node_id][dim];
+          }
+        }
+      });
+    return delta;
+  }
+
+  const CellValue<const double>
+  compute_delta_pu(const NodeValue<const TinyVector<Dimension>>& ur, const NodeValue<const TinyVector<Dimension>>& pr)
+  {
+    const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+    CellValue<double> delta{p_mesh->connectivity()};
+
+    const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+    delta.fill(0.);
+
+    parallel_for(
+      p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        const auto& cell_to_node = cell_to_node_matrix[cell_id];
+
+        for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+          const NodeId node_id      = cell_to_node[i_node];
+          TinyVector<Dimension> Fjr = zero;
+          for (size_t dim = 0; dim < Dimension; ++dim) {
+            Fjr[dim] = Cjr(cell_id, i_node)[dim] * pr[node_id][dim];
+          }
+          delta[cell_id] += dot(Fjr, ur[node_id]);
+        }
+      });
+    return delta;
+  }
+
+  const NodeValue<const bool>
+  getBoundaryNode(const BoundaryConditionList& bc_list)
+  {
+    NodeValue<bool> is_boundary_node{p_mesh->connectivity()};
+    is_boundary_node.fill(false);
+    for (auto&& bc_v : bc_list) {
+      std::visit(
+        [=](auto&& bc) {
+          using BCType = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+            auto node_list = bc.nodeList();
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              is_boundary_node[node_list[i_node]] = 1;
+            }
+          } else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
+            auto node_list = bc.nodeList();
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              is_boundary_node[node_list[i_node]] = 1;
+            }
+          }
+        },
+        bc_v);
+    }
+    return is_boundary_node;
+  }
+
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply(const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+  {
+    BoundaryConditionList bc_list = this->_getBCList(m_mesh, bc_descriptor_list);
+    const double lambda           = std::sqrt(dot(m_lambda_vector[0], m_lambda_vector[0]));
+    // const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+    const NodeValue<const bool> is_boundary_node = getBoundaryNode(bc_list);
+
+    const CellValue<double> rho              = copy(m_rho);
+    const CellValue<TinyVector<Dimension>> u = copy(m_u);
+    const CellValue<double> p                = copy(m_p);
+    const CellValue<double> E                = copy(m_E);
+
+    CellValue<TinyVector<Dimension>> u_np1 = copy(u);
+    CellValue<double> p_np1                = copy(p);
+    CellValue<double> E_np1                = copy(E);
+
+    const CellArray<const double> F_p                = compute_M_p(p, u, lambda);
+    const CellArray<const TinyVector<Dimension>> F_u = compute_M_u(u, p);
+
+    NodeArray<double> Fr_p                = compute_Flux<double>(F_p, is_boundary_node);
+    NodeArray<TinyVector<Dimension>> Fr_u = compute_Flux<TinyVector<Dimension>>(F_u, is_boundary_node);
+
+    const NodeValue<TinyVector<Dimension>> ur = compute_velocity_flux(Fr_p, lambda);
+    const NodeValue<TinyVector<Dimension>> pr = compute_pressure_flux(Fr_u);
+
+    for (auto&& bc_v : bc_list) {
+      std::visit(
+        [=](auto&& bc) {
+          using BCType = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<BCType, InflowListBoundaryCondition>) {
+            auto node_list       = bc.nodeList();
+            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];
+              for (size_t dim = 0; dim < Dimension; ++dim) {
+                ur[node_id][dim] = node_array_list[i_node][dim];
+                pr[node_id][dim] = node_array_list[i_node][Dimension];
+              }
+            }
+          }
+        },
+        bc_v);
+    }
+
+    const CellValue<const double> delta_u                = compute_delta_velocity(ur);
+    const CellValue<const TinyVector<Dimension>> delta_p = compute_delta_pressure(pr);
+    const CellValue<const double> delta_pu               = compute_delta_pu(ur, pr);
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        const double dt_over_Mj = m_dt * m_inv_Mj[cell_id];
+        u_np1[cell_id] -= dt_over_Mj * delta_p[cell_id];
+        E_np1[cell_id] -= dt_over_Mj * delta_pu[cell_id];
+      });
+
+    NodeValue<TinyVector<Dimension>> new_xr = copy(m_mesh.xr());
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { new_xr[node_id] += m_dt * ur[node_id]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(m_mesh.shared_connectivity(), new_xr);
+
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+    CellValue<double> rho_np1{p_mesh->connectivity()};
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rho_np1[cell_id] = m_Mj[cell_id] / new_Vj[cell_id]; });
+
+    return std::make_tuple(std::make_shared<MeshVariant>(new_mesh),
+                           std::make_shared<DiscreteFunctionVariant>(
+                             DiscreteFunctionP0<const double>(new_mesh, rho_np1)),
+                           std::make_shared<DiscreteFunctionVariant>(
+                             DiscreteFunctionP0<const TinyVector<Dimension>>(new_mesh, u_np1)),
+                           std::make_shared<DiscreteFunctionVariant>(
+                             DiscreteFunctionP0<const double>(new_mesh, E_np1)));
+  }
+
+  EulerKineticSolverLagrangeMultiD(std::shared_ptr<const MeshType> mesh,
+                                   const double& dt,
+                                   const std::vector<TinyVector<MeshType::Dimension>>& lambda_vector,
+                                   const double& gamma,
+                                   const size_t& space_order,
+                                   const size_t& time_order,
+                                   const CellValue<const double>& rho,
+                                   const CellValue<const TinyVector<MeshType::Dimension>>& u,
+                                   const CellValue<const double>& E,
+                                   const CellValue<const double>& p)
+    : m_dt{dt},
+      m_lambda_vector{lambda_vector},
+      m_gamma{gamma},
+      m_spaceOrder{space_order},
+      m_timeOrder{time_order},
+      m_rho{rho},
+      m_u{u},
+      m_E{E},
+      m_p{p},
+      p_mesh{mesh},
+      m_mesh{*mesh}
+  {
+    m_Mj = [&] {
+      CellValue<double> Mj{m_mesh.connectivity()};
+      parallel_for(
+        m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Mj[cell_id] = m_rho[cell_id] * m_Vj[cell_id]; });
+      return Mj;
+    }();
+
+    m_inv_Mj = [&] {
+      CellValue<double> inv_Mj{m_mesh.connectivity()};
+      parallel_for(
+        m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { inv_Mj[cell_id] = 1. / m_Mj[cell_id]; });
+      return inv_Mj;
+    }();
+
+    m_inv_Vj = [&] {
+      CellValue<double> inv_Vj{m_mesh.connectivity()};
+      parallel_for(
+        m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { inv_Vj[cell_id] = 1. / m_Vj[cell_id]; });
+      return inv_Vj;
+    }();
+  }
+
+  ~EulerKineticSolverLagrangeMultiD() = default;
+};
+
+template <MeshConcept MeshType>
+class EulerKineticSolverLagrangeMultiD<MeshType>::SymmetryBoundaryCondition
+{
+ public:
+  static constexpr const size_t Dimension = MeshType::Dimension;
+
+  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();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_flat_node_boundary.nodeList();
+  }
+
+  size_t
+  numberOfFaces() const
+  {
+    return m_mesh_flat_face_boundary.faceList().size();
+  }
+
+  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 <MeshConcept MeshType>
+class EulerKineticSolverLagrangeMultiD<MeshType>::WallBoundaryCondition
+{
+ public:
+  static constexpr const size_t Dimension = MeshType::Dimension;
+
+  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();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  size_t
+  numberOfFaces() const
+  {
+    return m_mesh_face_boundary.faceList().size();
+  }
+
+  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)
+  {
+    ;
+  }
+
+  ~WallBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class EulerKineticSolverLagrangeMultiD<MeshType>::InflowListBoundaryCondition
+{
+ public:
+  static constexpr const size_t Dimension = MeshType::Dimension;
+
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const Table<const double> m_node_array_list;
+  const Array<const NodeId> m_node_list;
+
+ public:
+  size_t
+  numberOfNodes() const
+  {
+    return m_node_list.size();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_node_list;
+  }
+
+  const Table<const double>&
+  nodeArrayList() const
+  {
+    return m_node_array_list;
+  }
+
+  InflowListBoundaryCondition(const Table<const double>& node_array_list, const Array<const NodeId> node_list)
+    : m_node_array_list(node_array_list), m_node_list(node_list)
+  {
+    ;
+  }
+
+  ~InflowListBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class EulerKineticSolverLagrangeMultiD<MeshType>::OutflowBoundaryCondition
+{
+ public:
+  static constexpr const size_t Dimension = MeshType::Dimension;
+
+  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();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  size_t
+  numberOfFaces() const
+  {
+    return m_mesh_face_boundary.faceList().size();
+  }
+
+  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)
+  {
+    ;
+  }
+
+  ~OutflowBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+std::tuple<std::shared_ptr<const MeshVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverLagrangeMultiD(
+  const double& dt,
+  const std::vector<TinyVector<Dimension>>& lambda_vector,
+  const double& gamma,
+  const double& eps,
+  const size_t& space_order,
+  const size_t& time_order,
+  const std::shared_ptr<const DiscreteFunctionVariant>& rho_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& u_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& E_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& p_n,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+{
+  const std::shared_ptr<const MeshVariant> mesh_v = getCommonMesh({rho_n, u_n, E_n, p_n});
+  if (mesh_v.use_count() == 0) {
+    throw NormalError("rho_n, u_n, E_n and p_n have to be defined on the same mesh");
+  }
+  if (space_order > 1 or time_order > 1) {
+    throw NotImplementedError("Euler kinetic solver in Multi D not implemented at order > 1");
+  }
+  double eps_bis = eps;   // A retirer !!
+  eps_bis += eps_bis;     // A retirer !!
+  return std::visit(
+    [&](auto&& p_mesh)
+      -> std::tuple<std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>> {
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
+      if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == Dimension)) {
+        auto rho = rho_n->get<DiscreteFunctionP0<const double>>();
+        auto u   = u_n->get<DiscreteFunctionP0<const TinyVector<MeshType::Dimension>>>();
+        auto E   = E_n->get<DiscreteFunctionP0<const double>>();
+        auto p   = p_n->get<DiscreteFunctionP0<const double>>();
+
+        EulerKineticSolverLagrangeMultiD solver(p_mesh, dt, lambda_vector, gamma, space_order, time_order,
+                                                rho.cellValues(), u.cellValues(), E.cellValues(), p.cellValues());
+
+        return solver.apply(bc_descriptor_list);
+      } else {
+        throw NotImplementedError("Invalid mesh type for multi-D Lagrangian Kinetic solver");
+      }
+    },
+    mesh_v->variant());
+}
+
+template std::tuple<std::shared_ptr<const MeshVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverLagrangeMultiD(
+  const double& dt,
+  const std::vector<TinyVector<1>>& lambda_vector,
+  const double& gamma,
+  const double& eps,
+  const size_t& space_order,
+  const size_t& time_order,
+  const std::shared_ptr<const DiscreteFunctionVariant>& rho_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& u_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& E_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& p_n,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+
+template std::tuple<std::shared_ptr<const MeshVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverLagrangeMultiD(
+  const double& dt,
+  const std::vector<TinyVector<2>>& lambda_vector,
+  const double& gamma,
+  const double& eps,
+  const size_t& space_order,
+  const size_t& time_order,
+  const std::shared_ptr<const DiscreteFunctionVariant>& rho_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& u_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& E_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& p_n,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+
+template std::tuple<std::shared_ptr<const MeshVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverLagrangeMultiD(
+  const double& dt,
+  const std::vector<TinyVector<3>>& lambda_vector,
+  const double& gamma,
+  const double& eps,
+  const size_t& space_order,
+  const size_t& time_order,
+  const std::shared_ptr<const DiscreteFunctionVariant>& rho_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& u_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& E_n,
+  const std::shared_ptr<const DiscreteFunctionVariant>& p_n,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
diff --git a/src/scheme/EulerKineticSolverLagrangeMultiD.hpp b/src/scheme/EulerKineticSolverLagrangeMultiD.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..acf8e58e1851eecc6a4ba2d58864a217a73e6acb
--- /dev/null
+++ b/src/scheme/EulerKineticSolverLagrangeMultiD.hpp
@@ -0,0 +1,27 @@
+#ifndef EULER_KINETIC_SOLVER_LAGRANGE_MULTID_HPP
+#define EULER_KINETIC_SOLVER_LAGRANGE_MULTID_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+class IBoundaryConditionDescriptor;
+
+template <size_t Dimension>
+std::tuple<std::shared_ptr<const MeshVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverLagrangeMultiD(
+  const double& dt,
+  const std::vector<TinyVector<Dimension>>& lambda_vector,
+  const double& gamma,
+  const double& eps,
+  const size_t& space_order,
+  const size_t& time_order,
+  const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+  const std::shared_ptr<const DiscreteFunctionVariant>& u,
+  const std::shared_ptr<const DiscreteFunctionVariant>& E,
+  const std::shared_ptr<const DiscreteFunctionVariant>& p,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+
+#endif   // EULER_KINETIC_SOLVER_LAGRANGE_MULTID_HPP