Fix Eucchlyd scheme for kinetic solver

src/scheme/EulerKineticSolver2D.cpp

@@ -352,9 +352,10 @@ class EulerKineticSolver2D
     const size_t nb_waves = m_lambda_vector.size();
     FaceArrayPerNode<double> 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();
     const auto& face_to_cell_matrix               = p_mesh->connectivity().faceToCellMatrix();
-    const auto& cell_to_node_matrix               = p_mesh->connectivity().cellToNodeMatrix();
+    const auto& face_cell_is_reversed             = p_mesh->connectivity().cellFaceIsReversed();
 
     parallel_for(
       p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
@@ -364,38 +365,27 @@ class EulerKineticSolver2D
         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 unsigned int i_node_face = node_local_number_in_its_face[i_face];
+            const size_t i_node_face                  = node_local_number_in_its_face[i_face];
             Fr[node_id][i_face][i_wave]               = 0;
             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 auto& cell_to_node = cell_to_node_matrix[cell_id];
-
-              auto local_node_number_in_cell = [&](NodeId node_number) {
-                for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
-                  if (node_number == cell_to_node[i_node]) {
-                    return i_node;
-                  }
-                }
-                return std::numeric_limits<size_t>::max();
-              };
-
-              const size_t i_node_cell = local_node_number_in_cell(node_id);
+              const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
 
               TinyVector<Dimension> n_face = nlr(face_id, i_node_face);
-              TinyVector<Dimension> n_node = njr(cell_id, i_node_cell);
+              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 (dot(n_face, n_node) < 0) {
-                for (size_t d = 0; d < Dimension; ++d) {
-                  n_face[d] = -n_face[d];
-                }
-              }
-              double li_nlr = dot(m_lambda_vector[i_wave], n_face);
               if (li_nlr > 0) {
                 Fr[node_id][i_face][i_wave] += Fn[cell_id][i_wave] * li_nlr;
                 sum += li_nlr;
+                // std::cout << "i_wave = " << i_wave << ", node_id = " << node_id << ", face_id = " << face_id << ",
+                // cell_id = " << cell_id << ", lambda = " << m_lambda_vector[i_wave] << ", nlr = " << n_face << ", njr
+                // = " << n_node << "\n";
               }
             }
             if (sum != 0) {
@@ -440,14 +430,13 @@ class EulerKineticSolver2D
     DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
 
     const auto& cell_to_face_matrix   = p_mesh->connectivity().cellToFaceMatrix();
-    const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
     const auto& face_to_node_matrix   = p_mesh->connectivity().faceToNodeMatrix();
     const auto& node_to_face_matrix   = p_mesh->connectivity().nodeToFaceMatrix();
+    const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
 
     parallel_for(
       p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
         const auto& cell_to_face = cell_to_face_matrix[cell_id];
-        const auto& cell_to_node = cell_to_node_matrix[cell_id];
 
         for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
           deltaLFn[cell_id][i_wave] = 0;
@@ -469,28 +458,12 @@ class EulerKineticSolver2D
                 return std::numeric_limits<size_t>::max();
               };
 
-              auto local_node_number_in_cell = [&](NodeId node_number) {
-                for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
-                  if (node_number == cell_to_node[i_node]) {
-                    return i_node;
-                  }
-                }
-                return std::numeric_limits<size_t>::max();
-              };
-
               const size_t i_face_node = local_face_number_in_node(face_id);
-              const size_t i_node_cell = local_node_number_in_cell(node_id);
 
-              TinyVector<Dimension> n_face = nlr(face_id, i_node_face);
-              TinyVector<Dimension> n_node = njr(cell_id, i_node_cell);
+              const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
 
-              if (dot(n_face, n_node) < 0) {
-                for (size_t d = 0; d < Dimension; ++d) {
-                  n_face[d] = -n_face[d];
-                }
-              }
+              const double li_Nlr = sign * dot(m_lambda_vector[i_wave], Nlr(face_id, i_node_face));
 
-              const double li_Nlr = dot(m_lambda_vector[i_wave], n_face);
               deltaLFn[cell_id][i_wave] += Fr[node_id][i_face_node][i_wave] * li_Nlr;
             }
           }