diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp index 0e10c1bf8c32f402e19eed70b48f434234d42632..788c6df649c37090a22efcfbd8a202be22d45cd9 100644 --- a/src/scheme/PolynomialReconstruction.cpp +++ b/src/scheme/PolynomialReconstruction.cpp @@ -278,6 +278,17 @@ PolynomialReconstruction::_build( } } + for (size_t i = 0; i < stencil_cell_list.size(); ++i) { + const CellId cell_i_id = stencil_cell_list[i]; + const double weight = 1. / l2Norm(xj[cell_i_id] - Xj); + for (size_t l = 0; l < basis_dimension - 1; ++l) { + A(i, l) *= weight; + } + for (size_t l = 0; l < number_of_columns; ++l) { + B(i, l) *= weight; + } + } + const SmallMatrix<double>& X = X_pool[t]; Givens::solveCollectionInPlace(A, X, B); diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp index 454e77c9d748399adcdb2fe2bbfef3e7174b4a3f..78f7bc5bad29c8da74249c21647d62763e9e6ee9 100644 --- a/src/scheme/PolynomialReconstruction.hpp +++ b/src/scheme/PolynomialReconstruction.hpp @@ -100,6 +100,15 @@ class PolynomialReconstruction } } + for (size_t i = 0; i < stencil_cell_list.size(); ++i) { + const CellId cell_i_id = stencil_cell_list[i]; + const double weight = 1. / l2Norm(xj[cell_i_id] - Xj); + for (size_t l = 0; l < MeshType::Dimension; ++l) { + A(i, l) *= weight; + } + b[i] *= weight; + } + SmallVector<double> x(MeshType::Dimension); Givens::solveInPlace(A, x, b); @@ -150,6 +159,17 @@ class PolynomialReconstruction } } + for (size_t i = 0; i < stencil_cell_list.size(); ++i) { + const CellId cell_i_id = stencil_cell_list[i]; + const double weight = 1. / l2Norm(xj[cell_i_id] - Xj); + for (size_t l = 0; l < MeshType::Dimension - 1; ++l) { + A(i, l) *= weight; + } + for (size_t l = 0; l < DataType::Dimension; ++l) { + B(i, l) *= weight; + } + } + TinyMatrix<MeshType::Dimension, DataType::Dimension> X; Givens::solveCollectionInPlace(A, X, B); @@ -205,6 +225,17 @@ class PolynomialReconstruction } } + for (size_t i = 0; i < stencil_cell_list.size(); ++i) { + const CellId cell_i_id = stencil_cell_list[i]; + const double weight = 1. / l2Norm(xj[cell_i_id] - Xj); + for (size_t l = 0; l < MeshType::Dimension - 1; ++l) { + A(i, l) *= weight; + } + for (size_t l = 0; l < DataType::Dimension; ++l) { + B(i, l) *= weight; + } + } + TinyMatrix<MeshType::Dimension, DataType::Dimension> X; Givens::solveCollectionInPlace(A, X, B); diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp index 7ce59a9109f1e13d040e46a69e3b2c2b76596790..e019cc0162e0ccf52a583020f8d8a8c9346809b9 100644 --- a/tests/test_PolynomialReconstruction.cpp +++ b/tests/test_PolynomialReconstruction.cpp @@ -506,7 +506,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1})); } - REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-13)); } } }