Add a space order 2 Fjr computation function

src/scheme/Order2AcousticSolver.cpp

@@ -322,22 +322,27 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
   _computeFjr(const MeshType& mesh,
               const NodeValuePerCell<const Rdxd>& Ajr,
               const NodeValue<const Rd>& ur,
-              const DiscreteVectorFunction& u,
-              const DiscreteScalarFunction& p) const
+              DiscreteFunctionDPk<Dimension,const Rd> DPk_uh,
+              DiscreteFunctionDPk<Dimension,const double> DPk_ph) const
   {
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
 
     const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+    const NodeValue<const Rd>& xr = mesh.xr();
 
-    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
 
     NodeValuePerCell<Rd> F{mesh.connectivity()};
     parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-        const auto& cell_nodes = cell_to_node_matrix[j];
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        const auto& node_to_cell                   = node_to_cell_matrix[r];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
 
-        for (size_t r = 0; r < cell_nodes.size(); ++r) {
-          F(j, r) = Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + p[j] * Cjr(j, r);
+        for(size_t j = 0; j < node_to_cell.size(); ++j){
+          const CellId J       = node_to_cell[j];
+          const unsigned int R = node_local_number_in_its_cells[j];
+          F(J,R) = Ajr(J,R) * (DPk_uh[J](xr[r]) - ur[r]) + DPk_ph[J](xr[r]) * Cjr(J,R);
         }
       });
 
@@ -400,7 +405,7 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     synchronize(br);
 
     NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
-    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
+    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, DPk_uh, DPk_ph);
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));