diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
index 447bab9c21d1e513b25a4f7faf8796efa89a3c49..397eb85e7d67a020138516acc5f602652bc063a1 100644
--- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
@@ -540,11 +540,12 @@ class RoeViscousFormEulerianCompositeSolver_v2
               // Assert(face_cell_list.size() == 1);
               // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
+              // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
+              Rd normal(zero);
+              int nbnormal = 0;
 
               for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
                 FaceId node_face_id = node_face_list[i_face];
-                Rd normal(zero);
 
                 for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) {
                   const FaceId facebc_id = face_list[i_facebc];
@@ -556,32 +557,38 @@ class RoeViscousFormEulerianCompositeSolver_v2
                     size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
                     // Normal locale approchée
-                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell);
+                    normalloc *= 1. / l2Norm(normalloc);
+                    normal += normalloc;
+                    ++nbnormal;
+
                     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);
+              if (nbnormal == 0)
+                continue;
+              normal *= 1. / nbnormal;
 
-                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 *= 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);
 
-                  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_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];
 
-                  vectorSym -= dot(vectorSym, normal) * normal;
+                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 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];
-                }
+                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 a2d092067e5910bc88e65d5bb9ffa16ec762c9b2..2c9e1711c951a3fb18f43aa2a1e73ec9094666aa 100644
--- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
@@ -547,11 +547,12 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2
               // Assert(face_cell_list.size() == 1);
               // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
+              // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
+              Rd normal(zero);
+              int nbnormal = 0;
 
               for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
                 FaceId node_face_id = node_face_list[i_face];
-                Rd normal(zero);
 
                 for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) {
                   const FaceId facebc_id = face_list[i_facebc];
@@ -563,31 +564,36 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2
                     size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
                     // Normal locale approchée
-                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell);
+                    normalloc *= 1. / l2Norm(normalloc);
+                    normal += normalloc;
+                    ++nbnormal;
+
                     break;
                   }
                 }
-                if (l2Norm(normal) == 0.)
-                  continue;
+              }
+              if (nbnormal == 0)
+                continue;
+              normal *= 1. / nbnormal;
 
-                normal *= 1. / l2Norm(normal);
-                const auto& node_cell_list = node_to_cell_matrix[node_id];
-                // Assert(face_cell_list.size() == 1);
-                const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
-                for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
-                  CellId node_cell_id              = node_cell_list[i_cell];
-                  size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];
+              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];
+                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;
+                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 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];
               }
             }
 
@@ -963,7 +969,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2
     const double gammaK,
     const double cK,
     const double pK,
-*/
+  */
                               const Rd& normal,
                               const bool check = false) const
   {
diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp
index f8d7311e550baf8235ef3710de131c2beb6b8ae2..2bedc88795d33628eb3c410e22543e1408c8230a 100644
--- a/src/scheme/RusanovEulerianCompositeSolver.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver.cpp
@@ -206,11 +206,12 @@ class RusanovEulerianCompositeSolver
               // Assert(face_cell_list.size() == 1);
               // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
+              // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
+              Rd normal(zero);
+              int nbnormal = 0;
 
               for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
                 FaceId node_face_id = node_face_list[i_face];
-                Rd normal(zero);
 
                 for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) {
                   const FaceId facebc_id = face_list[i_facebc];
@@ -222,32 +223,37 @@ class RusanovEulerianCompositeSolver
                     size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
                     // Normal locale approchée
-                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell);
+                    normalloc *= 1. / l2Norm(normalloc);
+                    normal += normalloc;
+                    ++nbnormal;
                     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);
+              if (nbnormal == 0)
+                continue;
+              normal *= 1. / nbnormal;
 
-                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 *= 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);
 
-                  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_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];
 
-                  vectorSym -= dot(vectorSym, normal) * normal;
+                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 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];
-                }
+                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 4383a00159535e16e211fcbe9218e2fb9566d19c..6fbaafb7a2813fc4f0a0e153d4001b9734c3f397 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp
@@ -214,11 +214,12 @@ class RusanovEulerianCompositeSolver_o2
               // Assert(face_cell_list.size() == 1);
               // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
+              // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
+              Rd normal(zero);
+              int nbnormal = 0;
 
               for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
                 FaceId node_face_id = node_face_list[i_face];
-                Rd normal(zero);
 
                 for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) {
                   const FaceId facebc_id = face_list[i_facebc];
@@ -230,34 +231,40 @@ class RusanovEulerianCompositeSolver_o2
                     size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
                     // Normal locale approchée
-                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell);
+                    normalloc *= 1. / l2Norm(normalloc);
+                    normal += normalloc;
+                    ++nbnormal;
                     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);
+              if (nbnormal == 0)
+                continue;
+              normal *= 1. / nbnormal;
 
-                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 *= 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);
 
-                  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_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];
 
-                  vectorSym -= dot(vectorSym, normal) * normal;
+                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 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];
-                }
+                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];
 
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
index 0fa8469dc38ade3b0d1f596d0634cc3d49f9f93b..af7c70e379caec3ff033ab3d91edb6a4b384b3e1 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
@@ -206,11 +206,11 @@ class RusanovEulerianCompositeSolver_v2
               // Assert(face_cell_list.size() == 1);
               // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
-
+              // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
+              Rd normal(zero);
+              int nbnormal = 0;
               for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
                 FaceId node_face_id = node_face_list[i_face];
-                Rd normal(zero);
 
                 for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) {
                   const FaceId facebc_id = face_list[i_facebc];
@@ -222,32 +222,36 @@ class RusanovEulerianCompositeSolver_v2
                     size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
                     // Normal locale approchée
-                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell);
+                    normalloc *= 1. / l2Norm(normalloc);
+                    normal += normalloc;
+                    ++nbnormal;
                     break;
                   }
                 }
-                if (l2Norm(normal) == 0.)
-                  continue;
+              }
+              if (nbnormal == 0)
+                continue;
+              normal *= 1. / nbnormal;
 
-                normal *= 1. / l2Norm(normal);
-                const auto& node_cell_list = node_to_cell_matrix[node_id];
-                // Assert(face_cell_list.size() == 1);
-                const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+              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];
+              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];
+                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;
+                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;
-                }
+                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;
               }
             }
 
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
index 42cb12e5ea8a0106bc617aaae4a0b524ee4e9d5b..be4c0495efd775e73e684b2effe2c23d1437a130 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
@@ -213,11 +213,12 @@ class RusanovEulerianCompositeSolver_v2_o2
               // Assert(face_cell_list.size() == 1);
               // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id);
 
-              // on va chercher la normale d'une face issue du noeud de CL et contenue dans le faceList
+              // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList
+              Rd normal(zero);
+              int nbnormal = 0;
 
               for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
                 FaceId node_face_id = node_face_list[i_face];
-                Rd normal(zero);
 
                 for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) {
                   const FaceId facebc_id = face_list[i_facebc];
@@ -229,32 +230,37 @@ class RusanovEulerianCompositeSolver_v2_o2
                     size_t face_local_number_in_cell = face_local_numbers_in_their_cells(facebc_id, 0);
 
                     // Normal locale approchée
-                    normal = Cjf(face_cell_id, face_local_number_in_cell);
+                    Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell);
+                    normalloc *= 1. / l2Norm(normalloc);
+                    normal += normalloc;
+                    ++nbnormal;
                     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);
+              if (nbnormal == 0)
+                continue;
+              normal *= 1. / nbnormal;
 
-                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 *= 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);
 
-                  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_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];
 
-                  vectorSym -= dot(vectorSym, normal) * normal;
+                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 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];
-                }
+                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];
               }
             }