Fix bugs in eucclhyd flux computation

src/scheme/EulerKineticSolver2D.cpp

@@ -384,26 +384,25 @@ class EulerKineticSolver2D
 
               const size_t i_node_cell = local_node_number_in_cell(node_id);
 
-              const TinyVector<Dimension> n_face = nlr(face_id, i_node_face);
-              const TinyVector<Dimension> n_node = njr(cell_id, i_node_cell);
+              TinyVector<Dimension> n_face = nlr(face_id, i_node_face);
+              TinyVector<Dimension> n_node = njr(cell_id, i_node_cell);
 
-              if (dot(n_face, n_node) > 0) {
-                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;
+              if (dot(n_face, n_node) < 0) {
+                for (size_t d = 0; d < Dimension; ++d) {
+                  n_face[d] = -n_face[d];
                 }
-              } else {
-                double li_nlr = -dot(m_lambda_vector[i_wave], n_face);
+              }
+              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;
               }
             }
-            }
+            if (sum != 0) {
               Fr[node_id][i_face][i_wave] /= sum;
             }
           }
+        }
       });
 
     return Fr;
@@ -441,12 +440,14 @@ 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();
 
     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;
@@ -455,8 +456,8 @@ class EulerKineticSolver2D
             const FaceId face_id     = cell_to_face[i_face_cell];
             const auto& face_to_node = face_to_node_matrix[face_id];
 
-            for (size_t i_node = 0; i_node < face_to_node.size(); ++i_node) {
-              const NodeId node_id     = face_to_node[i_node];
+            for (size_t i_node_face = 0; i_node_face < face_to_node.size(); ++i_node_face) {
+              const NodeId node_id     = face_to_node[i_node_face];
               const auto& node_to_face = node_to_face_matrix[node_id];
 
               auto local_face_number_in_node = [&](FaceId face_number) {
@@ -468,8 +469,28 @@ 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 double li_Nlr      = dot(m_lambda_vector[i_wave], Nlr(face_id, i_node));
+              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);
+
+              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 = dot(m_lambda_vector[i_wave], n_face);
               deltaLFn[cell_id][i_wave] += Fr[node_id][i_face_node][i_wave] * li_Nlr;
             }
           }