diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
index d56dc2a5632bfb6f40b9373b6e3ab063f95a30be..447bab9c21d1e513b25a4f7faf8796efa89a3c49 100644
--- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
@@ -519,48 +519,69 @@ class RoeViscousFormEulerianCompositeSolver_v2
             // 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 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_cell_list = node_to_cell_matrix[node_id];
+              const auto& node_face_list = node_to_face_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);
+              // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // Normal locale approchée
-              Rd normal(zero);
-              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];
-                normal += Cjr(node_cell_id, node_local_number_in_cell);
-              }
-              normal *= 1. / node_cell_list.size();
-              normal *= 1. / l2Norm(normal);
+              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
 
-              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 i_face = 0; i_face < node_face_list.size(); ++i_face) {
+                FaceId node_face_id = node_face_list[i_face];
+                Rd normal(zero);
 
-                Rd vectorSym(zero);
-                for (size_t dim = 0; dim < Dimension; ++dim)
-                  vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+                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);
 
-                vectorSym -= dot(vectorSym, normal) * normal;
+                    CellId face_cell_id              = face_cell_list[0];
+                    size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
-                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];
+                    // Normal locale approchée
+                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    break;
+                  }
+                }
+                if (l2Norm(normal) == 0.)
+                  continue;
+
+                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];
+                }
               }
             }
 
diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
index e707e3f8474387385fcb9a5d3a0b9b1f2c97054d..a2d092067e5910bc88e65d5bb9ffa16ec762c9b2 100644
--- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
@@ -527,48 +527,67 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2
             // 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 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_cell_list = node_to_cell_matrix[node_id];
+              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_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+              // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // Normal locale approchée
-              Rd normal(zero);
-              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];
-                normal += Cjr(node_cell_id, node_local_number_in_cell);
-              }
-              normal *= 1. / node_cell_list.size();
-              normal *= 1. / l2Norm(normal);
+              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
 
-              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 i_face = 0; i_face < node_face_list.size(); ++i_face) {
+                FaceId node_face_id = node_face_list[i_face];
+                Rd normal(zero);
 
-                Rd vectorSym(zero);
-                for (size_t dim = 0; dim < Dimension; ++dim)
-                  vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+                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);
 
-                vectorSym -= dot(vectorSym, normal) * normal;
+                    CellId face_cell_id              = face_cell_list[0];
+                    size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
-                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];
+                    // Normal locale approchée
+                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    break;
+                  }
+                }
+                if (l2Norm(normal) == 0.)
+                  continue;
+
+                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];
+                }
               }
             }
 
diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp
index cda7e625cdbc5f83542124d51d1f3fbfd4ad784b..f8d7311e550baf8235ef3710de131c2beb6b8ae2 100644
--- a/src/scheme/RusanovEulerianCompositeSolver.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver.cpp
@@ -185,48 +185,69 @@ class RusanovEulerianCompositeSolver
             // 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 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_cell_list = node_to_cell_matrix[node_id];
+              const auto& node_face_list = node_to_face_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);
+              // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // Normal locale approchée
-              Rd normal(zero);
-              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];
-                normal += Cjr(node_cell_id, node_local_number_in_cell);
-              }
-              normal *= 1. / node_cell_list.size();
-              normal *= 1. / l2Norm(normal);
+              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
 
-              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 i_face = 0; i_face < node_face_list.size(); ++i_face) {
+                FaceId node_face_id = node_face_list[i_face];
+                Rd normal(zero);
 
-                Rd vectorSym(zero);
-                for (size_t dim = 0; dim < Dimension; ++dim)
-                  vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+                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);
 
-                vectorSym -= dot(vectorSym, normal) * normal;
+                    CellId face_cell_id              = face_cell_list[0];
+                    size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
-                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];
+                    // Normal locale approchée
+                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    break;
+                  }
+                }
+                if (l2Norm(normal) == 0.)
+                  continue;
+
+                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];
+                }
               }
             }
 
diff --git a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp
index 7d53c2f884e925b3bba1c442c62ee3f13fd6533d..4383a00159535e16e211fcbe9218e2fb9566d19c 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp
@@ -193,51 +193,71 @@ class RusanovEulerianCompositeSolver_o2
             // 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 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_cell_list = node_to_cell_matrix[node_id];
+              const auto& node_face_list = node_to_face_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);
+              // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // Normal locale approchée
-              Rd normal(zero);
-              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];
-                normal += Cjr(node_cell_id, node_local_number_in_cell);
-              }
-              normal *= 1. / node_cell_list.size();
-              normal *= 1. / l2Norm(normal);
+              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
 
-              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 i_face = 0; i_face < node_face_list.size(); ++i_face) {
+                FaceId node_face_id = node_face_list[i_face];
+                Rd normal(zero);
 
-                Rd vectorSym(zero);
-                for (size_t dim = 0; dim < Dimension; ++dim)
-                  vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+                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);
 
-                vectorSym -= dot(vectorSym, normal) * normal;
+                    CellId face_cell_id              = face_cell_list[0];
+                    size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
-                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];
+                    // Normal locale approchée
+                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    break;
+                  }
+                }
+                if (l2Norm(normal) == 0.)
+                  continue;
+
+                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];
 
@@ -795,7 +815,8 @@ class RusanovEulerianCompositeSolver_o2
     DiscreteFunctionDPk eps_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
 
     // DiscreteFunctionDPk rho_u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
-    // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+    // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const
+    // double>>();
 
     CellValue<double> Limitor_rho(p_mesh->connectivity());
     CellValue<double> Limitor_u(p_mesh->connectivity());
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
index f367cb161b634b9b2a4eff54fa0dcfc2387e2376..0fa8469dc38ade3b0d1f596d0634cc3d49f9f93b 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
@@ -185,49 +185,69 @@ class RusanovEulerianCompositeSolver_v2
             // 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 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_cell_list = node_to_cell_matrix[node_id];
+              const auto& node_face_list = node_to_face_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);
+              // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // Normal locale approchée
-              Rd normal(zero);
-              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];
-                normal += Cjr(node_cell_id, node_local_number_in_cell);
-              }
+              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
 
-              normal *= 1. / node_cell_list.size();
-              normal *= 1. / l2Norm(normal);
+              for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
+                FaceId node_face_id = node_face_list[i_face];
+                Rd normal(zero);
 
-              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 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);
 
-                Rd vectorSym(zero);
-                for (size_t dim = 0; dim < Dimension; ++dim)
-                  vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+                    CellId face_cell_id              = face_cell_list[0];
+                    size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
-                vectorSym -= dot(vectorSym, normal) * normal;
+                    // Normal locale approchée
+                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    break;
+                  }
+                }
+                if (l2Norm(normal) == 0.)
+                  continue;
 
-                for (size_t dim = 0; dim < Dimension; ++dim)
-                  stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
-                // stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = 0;
+                normal *= 1. / l2Norm(normal);
+                const auto& node_cell_list = node_to_cell_matrix[node_id];
+                // Assert(face_cell_list.size() == 1);
+                const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+
+                for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+                  CellId node_cell_id              = node_cell_list[i_cell];
+                  size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];
+
+                  Rd vectorSym(zero);
+                  for (size_t dim = 0; dim < Dimension; ++dim)
+                    vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+
+                  vectorSym -= dot(vectorSym, normal) * normal;
+
+                  for (size_t dim = 0; dim < Dimension; ++dim)
+                    stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
+                  // stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = 0;
+                }
               }
             }
 
@@ -271,7 +291,7 @@ class RusanovEulerianCompositeSolver_v2
 
                 const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
 
-                // Normal locale approchée
+                // Normal locale approchée .. need better def
                 Rd normal(zero);
                 for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
                   CellId edge_cell_id              = edge_cell_list[i_cell];
@@ -302,6 +322,7 @@ class RusanovEulerianCompositeSolver_v2
         boundary_condition);
     }
   }
+
   void
   _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
                                  const MeshType& mesh,
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
index 7a22ea9c73f324fd65fe5598d386a51fce0e5695..42cb12e5ea8a0106bc617aaae4a0b524ee4e9d5b 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
@@ -193,48 +193,68 @@ class RusanovEulerianCompositeSolver_v2_o2
             // 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 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_cell_list = node_to_cell_matrix[node_id];
+              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_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+              // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // Normal locale approchée
-              Rd normal(zero);
-              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];
-                normal += Cjr(node_cell_id, node_local_number_in_cell);
-              }
-              normal *= 1. / node_cell_list.size();
-              normal *= 1. / l2Norm(normal);
+              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
 
-              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 i_face = 0; i_face < node_face_list.size(); ++i_face) {
+                FaceId node_face_id = node_face_list[i_face];
+                Rd normal(zero);
 
-                Rd vectorSym(zero);
-                for (size_t dim = 0; dim < Dimension; ++dim)
-                  vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+                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);
 
-                vectorSym -= dot(vectorSym, normal) * normal;
+                    CellId face_cell_id              = face_cell_list[0];
+                    size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
-                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];
+                    // Normal locale approchée
+                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    break;
+                  }
+                }
+                if (l2Norm(normal) == 0.)
+                  continue;
+
+                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];
+                }
               }
             }