diff --git a/src/scheme/Order2AcousticSolver.cpp b/src/scheme/Order2AcousticSolver.cpp index 6bdf83320555525311950b519e4f9192deda6d04..216d40088e9b3ae0ee14b4accd20eec3eaf4d3a3 100644 --- a/src/scheme/Order2AcousticSolver.cpp +++ b/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));