diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index 010ddd67e32e392a8b4f233db47bc92e1cfe56e0..f299da04d71838d25f34446caadb581e2e7f8ce3 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -46,8 +46,8 @@ getLambdaVector2D(const std::shared_ptr<const MeshVariant>& mesh,
           for (size_t m = 1; m < 4 * M + 1; ++m) {
             size_t i_wave         = (l - 1) * 4 * M + m - 1;
             double ll             = l;
-            lambda_vector[i_wave] = TinyVector<2>{(ll / L) * lambda * std::cos((m * pi) / (2 * M)),
-                                                  (ll / L) * lambda * std::sin((m * pi) / (2 * M))};
+            lambda_vector[i_wave] = TinyVector<2>{(ll / L) * lambda * std::cos((m * pi) / (2 * M) + 0.25 * pi),
+                                                  (ll / L) * lambda * std::sin((m * pi) / (2 * M) + 0.25 * pi)};
           }
         }
         // lambda_vector[4 * L * M] = TinyVector<2>{0, 0};   //!!!!!!!!!!!!!!!!!!!!
diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
index 8651c43999c61cb380933657cb1db271a9055ef1..1df0243085bd8ec08b0c007b7669662ef4a1d694 100644
--- a/src/scheme/EulerKineticSolverMultiDOrder2.cpp
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
@@ -454,7 +454,7 @@ class EulerKineticSolverMultiDOrder2
   FaceArray<T>
   compute_Flux_face(const DiscreteFunctionDPkVector<Dimension, const T>& DPk_Fn,
                     const CellArray<const T>& Fn,
-                    CellValue<bool>& is_wall_boundary_cell)
+                    FaceValue<bool>& is_wall_boundary_face)
   {
     const size_t nb_waves = m_lambda_vector.size();
     FaceArray<T> Fl(m_mesh.connectivity(), nb_waves);
@@ -484,7 +484,7 @@ class EulerKineticSolverMultiDOrder2
             double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
 
             if (li_nl > 0) {
-              if (not is_wall_boundary_cell[cell_id]) {
+              if ((not is_wall_boundary_face[face_id])) {   //} and (not is_tiny_vector_v<T>)) {
                 Fl[face_id][i_wave] += DPk_Fn(cell_id, i_wave)(m_xl[face_id]);
               } else {
                 Fl[face_id][i_wave] += Fn[cell_id][i_wave];
@@ -661,7 +661,7 @@ class EulerKineticSolverMultiDOrder2
 
                   for (size_t k_wave = 0; k_wave < nb_waves; ++k_wave) {
                     rhol += Fn_rho(cell_id, k_wave)(m_xl[face_id]);
-                    rho_ul += Fn_rho_u(cell_id, k_wave)(m_xl[face_id]);
+                    rho_ul += Fn_rho_u(cell_id, k_wave)(m_xl[face_id]);   // Fj_rho_u[cell_id][k_wave];   //
                     rho_El += Fn_rho_E(cell_id, k_wave)(m_xl[face_id]);
                   }
 
@@ -1676,25 +1676,19 @@ class EulerKineticSolverMultiDOrder2
         bc_v);
     }
 
-    CellValue<bool> is_wall_boundary_cell{p_mesh->connectivity()};
-    is_wall_boundary_cell.fill(false);
+    FaceValue<bool> is_wall_boundary_face{p_mesh->connectivity()};
+    is_wall_boundary_face.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, WallBoundaryCondition>) {
-            auto node_list = bc.nodeList();
-
-            is_wall_boundary_cell.fill(0);
+            auto face_list = bc.faceList();
 
-            auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+            is_wall_boundary_face.fill(0);
 
-            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_wall_boundary_cell[cell_id] = 1;
-              }
+            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              is_wall_boundary_face[face_list[i_face]] = 1;
             }
           }
         },
@@ -1729,23 +1723,59 @@ class EulerKineticSolverMultiDOrder2
 
     std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
     // Reconstruction uniquement à l'intérieur donc pas utile
-    // for (auto&& bc_descriptor : bc_descriptor_list) {
-    //   if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
-    //     symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
-    //   }
-    // }
-
-    PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
-                                                                 symmetry_boundary_descriptor_list);
-
-    auto reconstruction =
-      PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho),
-                                                                DiscreteFunctionP0Vector(p_mesh, Fn_rho_u),
-                                                                DiscreteFunctionP0Vector(p_mesh, Fn_rho_E));
-    auto DPk_Fn_rho = reconstruction[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+    for (auto&& bc_descriptor : bc_descriptor_list) {
+      if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
+        symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+      }
+    }
+
+    std::vector<std::shared_ptr<const IBoundaryDescriptor>> boundary_descriptor_list;
+
+    PolynomialReconstructionDescriptor reconstruction_descriptor_scalar(IntegrationMethodType::cell_center, 1,
+                                                                        symmetry_boundary_descriptor_list);
+
+    PolynomialReconstructionDescriptor reconstruction_descriptor_vector(IntegrationMethodType::cell_center, 1,
+                                                                        boundary_descriptor_list);
+
+    auto reconstruction_scalar =
+      PolynomialReconstruction{reconstruction_descriptor_scalar}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho),
+                                                                       // DiscreteFunctionP0Vector(p_mesh, Fn_rho_u),
+                                                                       DiscreteFunctionP0Vector(p_mesh, Fn_rho_E));
+
+    auto reconstruction_vector =
+      PolynomialReconstruction{reconstruction_descriptor_vector}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho_u));
+
+    auto DPk_Fn_rho   = reconstruction_scalar[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+    auto DPk_Fn_rho_E = reconstruction_scalar[1]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+
     auto DPk_Fn_rho_u =
-      reconstruction[1]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
-    auto DPk_Fn_rho_E = reconstruction[2]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+      reconstruction_vector[0]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
+
+    const auto& cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
+    const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+
+    for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+      const auto cell_to_face = cell_to_face_matrix[cell_id];
+      if ((cell_id >= 500) and (cell_id < 510)) {
+        size_t i_face = 2;   // for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+
+        FaceId face_id = cell_to_face[i_face];
+
+        const double sign             = face_cell_is_reversed(cell_id, i_face) ? -1 : 1;
+        const TinyVector<Dimension> n = sign * nl[face_id];
+
+        // for (size_t i_wave = 0; i_wave < m_lambda_vector.size(); ++i_wave) {
+        size_t i_wave = 0;
+        std::cout << "n = " << n << ", x = " << m_xj[cell_id] << ", lambda_1 = " << m_lambda_vector[i_wave] << "\n";
+        std::cout << "F^rho_" << i_wave << "(" << cell_id << "," << face_id
+                  << ") = " << DPk_Fn_rho(cell_id, i_wave)(m_xl[face_id]) << "; F^rho_u_" << i_wave << "(" << cell_id
+                  << "," << face_id << ") = " << DPk_Fn_rho_u(cell_id, i_wave)(m_xl[face_id]) << "; F^rho_E_" << i_wave
+                  << "(" << cell_id << "," << face_id << ") = " << DPk_Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) << "\n";
+        // }
+        //}
+      }
+    }
+    exit(0);
 
     DiscreteFunctionDPkVector<Dimension, double> Fn_rho_lim   = copy(DPk_Fn_rho);
     DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E_lim = copy(DPk_Fn_rho_E);
@@ -1762,10 +1792,10 @@ class EulerKineticSolverMultiDOrder2
     CellArray<double> deltaLFn_rho_E{p_mesh->connectivity(), nb_waves};
 
     if (m_scheme == 1) {
-      FaceArray<double> Fl_rho = compute_Flux_face<double>(Fn_rho_lim, Fn_rho, is_wall_boundary_cell);
+      FaceArray<double> Fl_rho = compute_Flux_face<double>(Fn_rho_lim, Fn_rho, is_wall_boundary_face);
       FaceArray<TinyVector<Dimension>> Fl_rho_u =
-        compute_Flux_face<TinyVector<Dimension>>(Fn_rho_u_lim, Fn_rho_u, is_wall_boundary_cell);
-      FaceArray<double> Fl_rho_E = compute_Flux_face<double>(Fn_rho_E_lim, Fn_rho_E, is_wall_boundary_cell);
+        compute_Flux_face<TinyVector<Dimension>>(Fn_rho_u_lim, Fn_rho_u, is_wall_boundary_face);
+      FaceArray<double> Fl_rho_E = compute_Flux_face<double>(Fn_rho_E_lim, Fn_rho_E, is_wall_boundary_face);
 
       apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, Fn_rho, Fn_rho_u, Fn_rho_E,
                     bc_list);
@@ -1820,18 +1850,40 @@ class EulerKineticSolverMultiDOrder2
     const CellArray<const double> M_rho_E                = compute_scalar_M(rho_E, A_rho_E);
 
     for (size_t i_subtimestep = 0; i_subtimestep < m_TimeOrder; i_subtimestep++) {
-      auto reconstruction_np1 =
-        PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0Vector(p_mesh, Fnp1_rho),
-                                                                  DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_u),
-                                                                  DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_E));
-      auto DPk_Fnp1_rho = reconstruction[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+      auto reconstruction_np1_scalar =
+        PolynomialReconstruction{reconstruction_descriptor_scalar}.build(DiscreteFunctionP0Vector(p_mesh, Fnp1_rho),
+                                                                         DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_E));
+      auto reconstruction_np1_vector =
+        PolynomialReconstruction{reconstruction_descriptor_vector}.build(DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_u));
+
+      auto DPk_Fnp1_rho = reconstruction_scalar[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+      auto DPk_Fnp1_rho_E =
+        reconstruction_scalar[1]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+
       auto DPk_Fnp1_rho_u =
-        reconstruction[1]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
-      auto DPk_Fnp1_rho_E = reconstruction[2]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+        reconstruction_vector[0]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
 
       DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_lim   = copy(DPk_Fnp1_rho);
       DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_E_lim = copy(DPk_Fnp1_rho_E);
 
+      // const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
+      // for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+      //   const auto cell_to_face = cell_to_face_matrix[cell_id];
+      //   if ((cell_id >= 500) and (cell_id < 510)) {
+      //     for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+      //       FaceId face_id = cell_to_face[i_face];
+
+      //       std::cout << "F^rho_1(" << cell_id << "," << face_id << ") = " << DPk_Fnp1_rho(cell_id,
+      //       0)(m_xl[face_id]); std::cout << "F^rho_u_1(" << cell_id << "," << face_id
+      //                 << ") = " << DPk_Fnp1_rho_u(cell_id, 0)(m_xl[face_id]);
+
+      //       std::cout << "F^rho_E_1(" << cell_id << "," << face_id
+      //                 << ") = " << DPk_Fnp1_rho_E(cell_id, 0)(m_xl[face_id]);
+      //     }
+      //   }
+      // }
+      // exit(0);
+
       scalar_limiter(Fnp1_rho, Fnp1_rho_lim);
       scalar_limiter(Fnp1_rho_E, Fnp1_rho_E_lim);
 
@@ -1844,10 +1896,10 @@ class EulerKineticSolverMultiDOrder2
       CellArray<double> deltaLFnp1_rho_E{p_mesh->connectivity(), nb_waves};
 
       if (m_scheme == 1) {
-        FaceArray<double> Fl_rho_np1 = compute_Flux_face<double>(Fnp1_rho_lim, Fn_rho, is_wall_boundary_cell);
+        FaceArray<double> Fl_rho_np1 = compute_Flux_face<double>(Fnp1_rho_lim, Fn_rho, is_wall_boundary_face);
         FaceArray<TinyVector<Dimension>> Fl_rho_u_np1 =
-          compute_Flux_face<TinyVector<Dimension>>(Fnp1_rho_u_lim, Fn_rho_u, is_wall_boundary_cell);
-        FaceArray<double> Fl_rho_E_np1 = compute_Flux_face<double>(Fnp1_rho_E_lim, Fn_rho_E, is_wall_boundary_cell);
+          compute_Flux_face<TinyVector<Dimension>>(Fnp1_rho_u_lim, Fn_rho_u, is_wall_boundary_face);
+        FaceArray<double> Fl_rho_E_np1 = compute_Flux_face<double>(Fnp1_rho_E_lim, Fn_rho_E, is_wall_boundary_face);
 
         apply_face_bc(Fl_rho_np1, Fl_rho_u_np1, Fl_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim, Fn_rho,
                       Fn_rho_u, Fn_rho_E, bc_list);