diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index 9aeb37209da9d5f340c1f19f4245eaa7aaf085c5..bc2fe33ba8ce85b84cdb418b71c44378c885b5ca 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -506,7 +506,7 @@ class EulerKineticSolverMultiD
 
   template <typename T>
   NodeArray<T>
-  compute_Flux1_glace(const CellArray<const T>& Fn, NodeValue<bool> is_symmetry_boundary_node)
+  compute_Flux1_glace(const CellArray<const T>& Fn, NodeValue<bool> is_boundary_node)
   {
     if constexpr (Dimension > 1) {
       const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
@@ -523,7 +523,7 @@ class EulerKineticSolverMultiD
 
       parallel_for(
         p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
-          if (not is_symmetry_boundary_node[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) {
@@ -553,13 +553,13 @@ class EulerKineticSolverMultiD
   }
 
   template <typename T>
-  FaceArrayPerNode<T>
-  compute_Flux1_eucclhyd(const CellArray<const T>& Fn, NodeValue<bool> is_symmetry_boundary_node)
+  NodeArray<T>
+  compute_Flux1_eucclhyd(const CellArray<const T>& Fn, NodeValue<bool> is_boundary_node)
   {
     if constexpr (Dimension > 1) {
       const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
       const size_t nb_waves                                   = m_lambda_vector.size();
-      FaceArrayPerNode<T> Flr(m_mesh.connectivity(), nb_waves);
+      NodeArray<T> Fr(m_mesh.connectivity(), nb_waves);
       const auto& node_local_numbers_in_their_faces = p_mesh->connectivity().nodeLocalNumbersInTheirFaces();
       const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
       const auto& node_to_face_matrix               = p_mesh->connectivity().nodeToFaceMatrix();
@@ -567,47 +567,49 @@ class EulerKineticSolverMultiD
       const auto& face_cell_is_reversed             = p_mesh->connectivity().cellFaceIsReversed();
 
       if constexpr (is_tiny_vector_v<T>) {
-        Flr.fill(zero);
+        Fr.fill(zero);
       } else {
-        Flr.fill(0.);
+        Fr.fill(0.);
       }
 
       parallel_for(
         p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
-          const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
-          const auto& node_to_face                  = node_to_face_matrix[node_id];
-
-          for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
-            for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
-              const FaceId face_id                      = node_to_face[i_face];
-              const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
-              const auto& face_to_cell                  = face_to_cell_matrix[face_id];
-              const size_t i_node_face                  = node_local_number_in_its_face[i_face];
+          if (not is_boundary_node[node_id]) {
+            const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
+            const auto& node_to_face                  = node_to_face_matrix[node_id];
 
+            for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
               double sum = 0;
 
-              for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
-                const CellId cell_id     = face_to_cell[i_cell];
-                const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
+              for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
+                const FaceId face_id                      = node_to_face[i_face];
+                const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
+                const auto& face_to_cell                  = face_to_cell_matrix[face_id];
+                const size_t i_node_face                  = node_local_number_in_its_face[i_face];
 
-                TinyVector<Dimension> n_face = nlr(face_id, i_node_face);
-                const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+                for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+                  const CellId cell_id     = face_to_cell[i_cell];
+                  const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
 
-                double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
+                  TinyVector<Dimension> n_face = nlr(face_id, i_node_face);
+                  const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+
+                  double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
 
-                if (li_nlr > 0) {
-                  Flr[node_id][i_face][i_wave] += li_nlr * Fn[cell_id][i_wave];
-                  sum += li_nlr;
+                  if (li_nlr > 0) {
+                    Fr[node_id][i_wave] += li_nlr * Fn[cell_id][i_wave];
+                    sum += li_nlr;
+                  }
                 }
               }
-              if ((sum != 0) and (not is_symmetry_boundary_node[node_id])) {
-                Flr[node_id][i_face][i_wave] = 1. / sum * Flr[node_id][i_face][i_wave];
+              if (sum != 0) {
+                Fr[node_id][i_wave] = 1. / sum * Fr[node_id][i_wave];
               }
             }
           }
         });
 
-      return Flr;
+      return Fr;
     } else {
       throw NormalError("Eucclhyd not defined in 1d");
     }
@@ -1030,23 +1032,28 @@ class EulerKineticSolverMultiD
   }
 
   void
-  apply_eucclhyd_bc(FaceArrayPerNode<double> Flr_rho,
-                    FaceArrayPerNode<TinyVector<Dimension>> Flr_rho_u,
-                    FaceArrayPerNode<double> Flr_rho_E,
+  apply_eucclhyd_bc(NodeArray<double> Fr_rho,
+                    NodeArray<TinyVector<Dimension>> Fr_rho_u,
+                    NodeArray<double> Fr_rho_E,
                     const CellValue<const double> rho,
                     const CellValue<const TinyVector<Dimension>> rho_u,
                     const CellValue<const double> rho_E,
+                    const CellArray<const double> Fn_rho,
+                    const CellArray<const TinyVector<Dimension>> Fn_rho_u,
+                    const CellArray<const double> Fn_rho_E,
                     const BoundaryConditionList& bc_list)
   {
     const size_t nb_waves = m_lambda_vector.size();
 
-    // const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
-    const auto& face_local_numbers_in_their_cells           = p_mesh->connectivity().faceLocalNumbersInTheirCells();
-    const auto& face_local_numbers_in_their_nodes           = p_mesh->connectivity().faceLocalNumbersInTheirNodes();
-    const auto& face_to_node_matrix                         = p_mesh->connectivity().faceToNodeMatrix();
-    const auto& face_to_cell_matrix                         = p_mesh->connectivity().faceToCellMatrix();
     const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
     const auto& face_cell_is_reversed                       = p_mesh->connectivity().cellFaceIsReversed();
+    const auto& node_local_numbers_in_their_cells           = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
+    const auto& node_local_numbers_in_their_faces           = p_mesh->connectivity().nodeLocalNumbersInTheirFaces();
+    const auto& face_local_numbers_in_their_cells           = p_mesh->connectivity().faceLocalNumbersInTheirCells();
+    const auto& face_to_cell_matrix                         = p_mesh->connectivity().faceToCellMatrix();
+    const auto& node_to_cell_matrix                         = p_mesh->connectivity().nodeToCellMatrix();
+    const auto& node_to_face_matrix                         = p_mesh->connectivity().nodeToFaceMatrix();
+    const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
 
     TinyVector<MeshType::Dimension> inv_S = zero;
     for (size_t d = 0; d < MeshType::Dimension; ++d) {
@@ -1055,7 +1062,6 @@ class EulerKineticSolverMultiD
       }
       inv_S[d] = 1. / inv_S[d];
     }
-
     FaceArrayPerNode<double> sum_li_nlr(m_mesh.connectivity(), nb_waves);
     sum_li_nlr.fill(0);
 
@@ -1063,42 +1069,64 @@ class EulerKineticSolverMultiD
       std::visit(
         [=, this](auto&& bc) {
           using BCType = std::decay_t<decltype(bc)>;
-          if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
-            auto face_list          = bc.faceList();
-            auto n                  = bc.outgoingNormal();
-            auto nxn                = tensorProduct(n, n);
+          if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
+            auto node_list = bc.nodeList();
+
             TinyMatrix<Dimension> I = identity;
-            auto txt                = I - nxn;
 
-            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-              const FaceId face_id = face_list[i_face];
+            NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()};
+            nr.fill(zero);
 
-              const auto& face_to_node                  = face_to_node_matrix[face_id];
-              const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
-              const auto& face_to_cell                  = face_to_cell_matrix[face_id];
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              const NodeId node_id                      = node_list[i_node];
+              const auto& node_to_cell                  = node_to_cell_matrix[node_id];
+              const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+
+              for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+                const CellId cell_id     = node_to_cell[i_cell];
+                const size_t i_node_cell = node_local_number_in_its_cell[i_cell];
+
+                nr[node_id] += njr(cell_id, i_node_cell);
+              }
+              nr[node_id] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id];
+              auto nxn    = tensorProduct(nr[node_id], nr[node_id]);
+              auto txt    = I - nxn;
 
               for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
-                for (size_t i_node = 0; i_node < face_to_node.size(); ++i_node) {
-                  const NodeId node_id = face_to_node[i_node];
+                const auto& node_to_face                  = node_to_face_matrix[node_id];
+                const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
+
+                double sum                = 0;
+                Fr_rho[node_id][i_wave]   = 0;
+                Fr_rho_u[node_id][i_wave] = zero;
+                Fr_rho_E[node_id][i_wave] = 0;
 
-                  const auto& face_local_number_in_its_node = face_local_numbers_in_their_nodes.itemArray(face_id);
-                  const size_t i_face_node                  = face_local_number_in_its_node[i_node];
+                for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
+                  const FaceId face_id                      = node_to_face[i_face];
+                  const auto& face_to_cell                  = face_to_cell_matrix[face_id];
+                  const size_t i_node_face                  = node_local_number_in_its_face[i_face];
+                  const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
 
                   for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
                     const CellId cell_id     = face_to_cell[i_cell];
                     const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
 
-                    TinyVector<Dimension> n_face = nlr(face_id, i_node);
                     const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+                    TinyVector<Dimension> n_face = sign * nlr(face_id, i_node_face);
 
-                    TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
-                    double li_nlr_prime             = sign * dot(m_lambda_vector[i_wave], nlr_prime);
-                    double li_nlr                   = sign * dot(m_lambda_vector[i_wave], n_face);
+                    double li_nlr = dot(m_lambda_vector[i_wave], n_face);
 
                     if (li_nlr > 0) {
-                      sum_li_nlr[node_id][i_face_node][i_wave] += li_nlr;
+                      Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr;
+                      Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u[cell_id][i_wave];
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr;
+
+                      sum += li_nlr;
                     }
 
+                    TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
+                    double li_nlr_prime             = dot(m_lambda_vector[i_wave], nlr_prime);
+
                     if (li_nlr_prime > 0) {
                       double rhoj_prime                  = rho[cell_id];
                       TinyVector<Dimension> rho_uj_prime = txt * rho_u[cell_id] - nxn * rho_u[cell_id];
@@ -1123,48 +1151,19 @@ class EulerKineticSolverMultiD
                         }
                         Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1];
                       }
-                      Flr_rho[node_id][i_face_node][i_wave] += Fn_rho_prime * li_nlr_prime;
-                      Flr_rho_u[node_id][i_face_node][i_wave] += li_nlr_prime * Fn_rho_u_prime;
-                      Flr_rho_E[node_id][i_face_node][i_wave] += Fn_rho_E_prime * li_nlr_prime;
 
-                      sum_li_nlr[node_id][i_face_node][i_wave] += li_nlr_prime;
+                      Fr_rho[node_id][i_wave] += Fn_rho_prime * li_nlr_prime;
+                      Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u_prime;
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_nlr_prime;
+
+                      sum += li_nlr_prime;
                     }
                   }
                 }
-              }
-            }
-          }
-        },
-        bc_v);
-    }
-
-    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 face_list = bc.faceList();
-
-            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];
-
-              for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
-                for (size_t i_node = 0; i_node < face_to_node.size(); ++i_node) {
-                  const NodeId node_id = face_to_node[i_node];
-
-                  const auto& face_local_number_in_its_node = face_local_numbers_in_their_nodes.itemArray(face_id);
-                  const size_t i_face_node                  = face_local_number_in_its_node[i_node];
-
-                  if (sum_li_nlr[node_id][i_face_node][i_wave] != 0) {
-                    Flr_rho[node_id][i_face_node][i_wave] =
-                      1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho[node_id][i_face_node][i_wave];
-                    Flr_rho_u[node_id][i_face_node][i_wave] =
-                      1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho_u[node_id][i_face_node][i_wave];
-                    Flr_rho_E[node_id][i_face_node][i_wave] =
-                      1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho[node_id][i_face_node][i_wave];
-                  }
+                if (sum != 0) {
+                  Fr_rho[node_id][i_wave] /= sum;
+                  Fr_rho_u[node_id][i_wave] = 1. / sum * Fr_rho_u[node_id][i_wave];
+                  Fr_rho_E[node_id][i_wave] /= sum;
                 }
               }
             }
@@ -1176,7 +1175,7 @@ class EulerKineticSolverMultiD
 
   template <typename T>
   CellArray<T>
-  compute_deltaLFn_glace(NodeArray<T> Fr)
+  compute_deltaLFn_node(NodeArray<T> Fr)
   {
     if constexpr (Dimension > 1) {
       const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
@@ -1206,7 +1205,7 @@ class EulerKineticSolverMultiD
         });
       return deltaLFn;
     } else {
-      throw NormalError("Glace not defined in 1d");
+      throw NormalError("Nodal solver not defined in 1d");
     }
   }
 
@@ -1373,8 +1372,8 @@ class EulerKineticSolverMultiD
         bc_v);
     }
 
-    NodeValue<bool> is_symmetry_boundary_node{p_mesh->connectivity()};
-    is_symmetry_boundary_node.fill(false);
+    NodeValue<bool> is_boundary_node{p_mesh->connectivity()};
+    is_boundary_node.fill(false);
     for (auto&& bc_v : bc_list) {
       std::visit(
         [=](auto&& bc) {
@@ -1382,7 +1381,12 @@ class EulerKineticSolverMultiD
           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_symmetry_boundary_node[node_list[i_node]] = 1;
+              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;
             }
           }
         },
@@ -1411,35 +1415,31 @@ class EulerKineticSolverMultiD
     const CellValue<const TinyVector<Dimension>> rho_u = compute_PFn<TinyVector<Dimension>>(Fn_rho_u);
     const CellValue<const double> rho_E                = compute_PFn<double>(Fn_rho_E);
 
-    NodeArray<double> Fr_rho = compute_Flux1_glace<double>(Fn_rho, is_symmetry_boundary_node);
-    NodeArray<TinyVector<Dimension>> Fr_rho_u =
-      compute_Flux1_glace<TinyVector<Dimension>>(Fn_rho_u, is_symmetry_boundary_node);
-    NodeArray<double> Fr_rho_E = compute_Flux1_glace<double>(Fn_rho_E, is_symmetry_boundary_node);
+    // NodeArray<double> Fr_rho                  = compute_Flux1_glace<double>(Fn_rho, is_boundary_node);
+    // NodeArray<TinyVector<Dimension>> Fr_rho_u = compute_Flux1_glace<TinyVector<Dimension>>(Fn_rho_u,
+    // is_boundary_node); NodeArray<double> Fr_rho_E                = compute_Flux1_glace<double>(Fn_rho_E,
+    // is_boundary_node);
 
-    // FaceArrayPerNode<double> Flr_rho = compute_Flux1_eucclhyd<double>(Fn_rho, is_symmetry_boundary_node);
-    // FaceArrayPerNode<TinyVector<Dimension>> Flr_rho_u =
-    //   compute_Flux1_eucclhyd<TinyVector<Dimension>>(Fn_rho_u, is_symmetry_boundary_node);
-    // FaceArrayPerNode<double> Flr_rho_E = compute_Flux1_eucclhyd<double>(Fn_rho_E, is_symmetry_boundary_node);
+    NodeArray<double> Fr_rho = compute_Flux1_eucclhyd<double>(Fn_rho, is_boundary_node);
+    NodeArray<TinyVector<Dimension>> Fr_rho_u =
+      compute_Flux1_eucclhyd<TinyVector<Dimension>>(Fn_rho_u, is_boundary_node);
+    NodeArray<double> Fr_rho_E = compute_Flux1_eucclhyd<double>(Fn_rho_E, is_boundary_node);
 
     // FaceArray<double> Fl_rho                  = compute_Flux1_face<double>(Fn_rho);
     // FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face<TinyVector<Dimension>>(Fn_rho_u);
     // FaceArray<double> Fl_rho_E                = compute_Flux1_face<double>(Fn_rho_E);
 
     //    apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
-    apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
-    // apply_eucclhyd_bc(Flr_rho, Flr_rho_u, Flr_rho_E, rho, rho_u, rho_E, bc_list);
+    // apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
+    apply_eucclhyd_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
 
     synchronize(Fr_rho);
     synchronize(Fr_rho_u);
     synchronize(Fr_rho_E);
 
-    const CellArray<const double> deltaLFn_rho                  = compute_deltaLFn_glace(Fr_rho);
-    const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_glace(Fr_rho_u);
-    const CellArray<const double> deltaLFn_rho_E                = compute_deltaLFn_glace(Fr_rho_E);
-
-    // const CellArray<const double> deltaLFn_rho                  = compute_deltaLFn_eucclhyd(Flr_rho);
-    // const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_eucclhyd(Flr_rho_u);
-    // const CellArray<const double> deltaLFn_rho_E                = compute_deltaLFn_eucclhyd(Flr_rho_E);
+    const CellArray<const double> deltaLFn_rho                  = compute_deltaLFn_node(Fr_rho);
+    const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u);
+    const CellArray<const double> deltaLFn_rho_E                = compute_deltaLFn_node(Fr_rho_E);
 
     // const CellArray<const double> deltaLFn_rho = compute_deltaLFn_face<double>(Fl_rho);
     // const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u =