diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index b9b789cfa8e3d7bc5328062c6085217dbe142120..6b775344aefb0a10e1c45b7df73bcd64bd5ba010 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -65,6 +65,9 @@
#include <utils/checkpointing/WriteIQuadratureDescriptor.hpp>
#include <utils/checkpointing/WriteVariableBCDescriptor.hpp>
+#warning remove
+#include <scheme/test_developed_reconstruction.hpp>
+
#include <memory>
SchemeModule::SchemeModule()
@@ -756,6 +759,17 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("test_developed_reconstruction",
+ std::function(
+
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E) -> void {
+ test_developed_reconstruction(rho, u, E);
+ }
+
+ ));
+
MathFunctionRegisterForVh{*this};
}
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 8c174aab69d2c588b143d316f24fc10357c36598..5ff646da8e42d42fa3e289158fd0bc0ce76c2e63 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -16,6 +16,8 @@ add_library(
PolynomialBasisData.cpp
PolynomialBasisDataManager.cpp
PolynomialReconstruction.cpp
+
+ test_developed_reconstruction.cpp
)
target_link_libraries(
diff --git a/src/scheme/test_developed_reconstruction.cpp b/src/scheme/test_developed_reconstruction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6a4e4e6a030724dc671e90e534acb61c4d948b71
--- /dev/null
+++ b/src/scheme/test_developed_reconstruction.cpp
@@ -0,0 +1,967 @@
+#include <scheme/test_developed_reconstruction.hpp>
+
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/PolynomialBasisData.hpp>
+#include <scheme/PolynomialBasisDataManager.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <scheme/PolynomialReconstructionDescriptor.hpp>
+#include <utils/SmallArray.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
+
+#include <fstream>
+
+template <MeshConcept MeshType, typename DataType>
+void computeFGBar(const std::shared_ptr<const MeshType>& p_mesh,
+ const DiscreteFunctionDPk<MeshType::Dimension, const double>& f_bar,
+ const DiscreteFunctionDPk<MeshType::Dimension, const DataType>& g_bar,
+ DiscreteFunctionDPk<MeshType::Dimension, DataType>& fg_bar);
+
+template <typename DataType>
+void
+computeFGBar(const std::shared_ptr<const Mesh<1>>& p_mesh,
+ const DiscreteFunctionDPk<1, const double>& f_bar,
+ const DiscreteFunctionDPk<1, const DataType>& g_bar,
+ DiscreteFunctionDPk<1, DataType>& fg_bar)
+{
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto f_bar_coeffs = f_bar.coefficients(cell_id);
+ auto g_bar_coeffs = g_bar.coefficients(cell_id);
+ auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+
+ for (size_t i = 0; i < fg_bar_coeffs.size(); ++i) {
+ DataType sum;
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ sum = 0;
+ } else {
+ sum = zero;
+ }
+ for (size_t i2 = 0; i2 <= i; ++i2) {
+ sum += f_bar_coeffs[i - i2] * g_bar_coeffs[i2];
+ }
+ fg_bar_coeffs[i] = sum;
+ }
+ });
+};
+
+template <typename DataType>
+void
+computeFGBar(const std::shared_ptr<const Mesh<2>>& p_mesh,
+ const DiscreteFunctionDPk<2, const double>& f_bar,
+ const DiscreteFunctionDPk<2, const DataType>& g_bar,
+ DiscreteFunctionDPk<2, DataType>& fg_bar)
+{
+ Assert((f_bar.degree() == g_bar.degree()) and (f_bar.degree() == fg_bar.degree()));
+ const size_t degree = f_bar.degree();
+
+ SmallArray<size_t> y_row_index(degree + 1);
+ {
+ size_t i_y = 0;
+
+ y_row_index[i_y++] = 0;
+ for (ssize_t n = degree + 1; n > 1; --n) {
+ y_row_index[i_y] = y_row_index[i_y - 1] + n;
+ ++i_y;
+ }
+ }
+ auto index = [&](size_t i, size_t j) { return y_row_index[j] + i; };
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto f_bar_coeffs = f_bar.coefficients(cell_id);
+ auto g_bar_coeffs = g_bar.coefficients(cell_id);
+ auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+
+ for (size_t j = 0; j <= degree; ++j) {
+ for (size_t i = 0; i <= degree - j; ++i) {
+ DataType sum;
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ sum = 0;
+ } else {
+ sum = zero;
+ }
+
+ for (size_t j2 = 0; j2 <= j; ++j2) {
+ for (size_t i2 = 0; i2 <= i; ++i2) {
+ sum += f_bar_coeffs[index(i - i2, j - j2)] * g_bar_coeffs[index(i2, j2)];
+ }
+ }
+
+ fg_bar_coeffs[index(i, j)] = sum;
+ }
+ }
+ });
+}
+
+template <typename DataType>
+void
+computeFGBar(const std::shared_ptr<const Mesh<3>>& p_mesh,
+ const DiscreteFunctionDPk<3, const double>& f_bar,
+ const DiscreteFunctionDPk<3, const DataType>& g_bar,
+ DiscreteFunctionDPk<3, DataType>& fg_bar)
+{
+ Assert((f_bar.degree() == g_bar.degree()) and (f_bar.degree() == fg_bar.degree()));
+
+ const size_t degree = f_bar.degree();
+
+ SmallArray<size_t> yz_row_index((degree + 2) * (degree + 1) / 2 + 1);
+ SmallArray<size_t> z_triangle_index(degree + 1);
+
+ {
+ size_t i_z = 0;
+ size_t i_yz = 0;
+
+ yz_row_index[i_yz++] = 0;
+ for (ssize_t n = degree + 1; n >= 1; --n) {
+ z_triangle_index[i_z++] = i_yz - 1;
+ for (ssize_t i = n; i >= 1; --i) {
+ yz_row_index[i_yz] = yz_row_index[i_yz - 1] + i;
+ ++i_yz;
+ }
+ }
+ }
+
+ SmallArray<size_t> yz_row_size{yz_row_index.size() - 1};
+ for (size_t i = 0; i < yz_row_size.size(); ++i) {
+ yz_row_size[i] = yz_row_index[i + 1] - yz_row_index[i];
+ }
+
+ auto index = [&](size_t i, size_t j, size_t k) { return yz_row_index[z_triangle_index[k] + j] + i; };
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto f_bar_coeffs = f_bar.coefficients(cell_id);
+ auto g_bar_coeffs = g_bar.coefficients(cell_id);
+ auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+
+ for (size_t k = 0; k <= degree; ++k) {
+ for (size_t j = 0; j <= degree - k; ++j) {
+ for (size_t i = 0; i <= degree - j - k; ++i) {
+ DataType sum;
+ if constexpr (std::is_arithmetic_v<DataType>) {
+ sum = 0;
+ } else {
+ sum = zero;
+ }
+
+ for (size_t k2 = 0; k2 <= k; ++k2) {
+ for (size_t j2 = 0; j2 <= j; ++j2) {
+ for (size_t i2 = 0; i2 <= i; ++i2) {
+ sum += f_bar_coeffs[index(i - i2, j - j2, k - k2)] * g_bar_coeffs[index(i2, j2, k2)];
+ }
+ }
+ }
+
+ fg_bar_coeffs[index(i, j, k)] = sum;
+ }
+ }
+ }
+ });
+}
+
+template <MeshConcept MeshType, typename DataType>
+void computeKappaBar(const PolynomialReconstructionDescriptor descriptor,
+ const std::shared_ptr<const MeshType>& p_mesh,
+ const DiscreteFunctionDPk<MeshType::Dimension, const DataType>& u_bar,
+ const DiscreteFunctionP0<const double>& kappa_j,
+ DiscreteFunctionDPk<MeshType::Dimension, double>& kappa_bar);
+
+template <typename DataType>
+void
+computeKappaBar(const PolynomialReconstructionDescriptor descriptor,
+ const std::shared_ptr<const Mesh<1>>& p_mesh,
+ const DiscreteFunctionDPk<1, const DataType>& u_bar,
+ const DiscreteFunctionP0<const double>& kappa_j,
+ DiscreteFunctionDPk<1, double>& kappa_bar)
+{
+ auto polynomial_basis_means = PolynomialBasisDataManager::instance()
+ .getPolynomialBasisData(*p_mesh)
+ .getPolynomialBasisMeans(descriptor.integrationMethodType(), descriptor.degree());
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto u_bar_coeffs = u_bar.coefficients(cell_id);
+ auto kappa_bar_coeffs = kappa_bar.coefficients(cell_id);
+ auto basis_means = polynomial_basis_means[cell_id];
+
+ for (size_t i = 0; i < kappa_bar_coeffs.size(); ++i) {
+ double sum = 0;
+ for (size_t i2 = 0; i2 <= i; ++i2) {
+ sum += dot(u_bar_coeffs[i - i2], u_bar_coeffs[i2]);
+ }
+
+ kappa_bar_coeffs[i] = 0.5 * sum;
+ }
+
+ double kappa_bar_mean = 0;
+ for (size_t i = 0; i < basis_means.size(); ++i) {
+ kappa_bar_mean += basis_means[i] * kappa_bar_coeffs[i];
+ }
+
+ kappa_bar_coeffs[0] += kappa_j[cell_id] - kappa_bar_mean;
+ });
+}
+
+template <typename DataType>
+void
+computeKappaBar(const PolynomialReconstructionDescriptor descriptor,
+ const std::shared_ptr<const Mesh<2>>& p_mesh,
+ const DiscreteFunctionDPk<2, const DataType>& u_bar,
+ const DiscreteFunctionP0<const double>& kappa_j,
+ DiscreteFunctionDPk<2, double>& kappa_bar)
+{
+ Assert((u_bar.degree() == kappa_bar.degree()) and (u_bar.degree() == descriptor.degree()));
+ const size_t degree = descriptor.degree();
+
+ SmallArray<size_t> y_row_index(degree + 1);
+ {
+ size_t i_y = 0;
+
+ y_row_index[i_y++] = 0;
+ for (ssize_t n = degree + 1; n > 1; --n) {
+ y_row_index[i_y] = y_row_index[i_y - 1] + n;
+ ++i_y;
+ }
+ }
+ auto index = [&](size_t i, size_t j) { return y_row_index[j] + i; };
+
+ auto polynomial_basis_means = PolynomialBasisDataManager::instance()
+ .getPolynomialBasisData(*p_mesh)
+ .getPolynomialBasisMeans(descriptor.integrationMethodType(), descriptor.degree());
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto u_bar_coeffs = u_bar.coefficients(cell_id);
+ auto kappa_bar_coeffs = kappa_bar.coefficients(cell_id);
+ auto basis_means = polynomial_basis_means[cell_id];
+
+ for (size_t j = 0; j <= degree; ++j) {
+ for (size_t i = 0; i <= degree - j; ++i) {
+ double sum = 0;
+
+ for (size_t j2 = 0; j2 <= j; ++j2) {
+ for (size_t i2 = 0; i2 <= i; ++i2) {
+ sum += dot(u_bar_coeffs[index(i - i2, j - j2)], u_bar_coeffs[index(i2, j2)]);
+ }
+ }
+
+ kappa_bar_coeffs[index(i, j)] = 0.5 * sum;
+ }
+ }
+
+ double kappa_bar_mean = 0;
+ for (size_t i = 0; i < basis_means.size(); ++i) {
+ kappa_bar_mean += basis_means[i] * kappa_bar_coeffs[i];
+ }
+
+ kappa_bar_coeffs[0] += kappa_j[cell_id] - kappa_bar_mean;
+ });
+}
+
+template <typename DataType>
+void
+computeKappaBar(const PolynomialReconstructionDescriptor descriptor,
+ const std::shared_ptr<const Mesh<3>>& p_mesh,
+ const DiscreteFunctionDPk<3, const DataType>& u_bar,
+ const DiscreteFunctionP0<const double>& kappa_j,
+ DiscreteFunctionDPk<3, double>& kappa_bar)
+{
+ Assert((u_bar.degree() == kappa_bar.degree()) and (u_bar.degree() == descriptor.degree()));
+
+ const size_t degree = descriptor.degree();
+
+ SmallArray<size_t> yz_row_index((degree + 2) * (degree + 1) / 2 + 1);
+ SmallArray<size_t> z_triangle_index(degree + 1);
+
+ {
+ size_t i_z = 0;
+ size_t i_yz = 0;
+
+ yz_row_index[i_yz++] = 0;
+ for (ssize_t n = degree + 1; n >= 1; --n) {
+ z_triangle_index[i_z++] = i_yz - 1;
+ for (ssize_t i = n; i >= 1; --i) {
+ yz_row_index[i_yz] = yz_row_index[i_yz - 1] + i;
+ ++i_yz;
+ }
+ }
+ }
+
+ SmallArray<size_t> yz_row_size{yz_row_index.size() - 1};
+ for (size_t i = 0; i < yz_row_size.size(); ++i) {
+ yz_row_size[i] = yz_row_index[i + 1] - yz_row_index[i];
+ }
+
+ auto index = [&](size_t i, size_t j, size_t k) { return yz_row_index[z_triangle_index[k] + j] + i; };
+
+ auto polynomial_basis_means = PolynomialBasisDataManager::instance()
+ .getPolynomialBasisData(*p_mesh)
+ .getPolynomialBasisMeans(descriptor.integrationMethodType(), descriptor.degree());
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto u_bar_coeffs = u_bar.coefficients(cell_id);
+ auto kappa_bar_coeffs = kappa_bar.coefficients(cell_id);
+ auto basis_means = polynomial_basis_means[cell_id];
+
+ for (size_t k = 0; k <= degree; ++k) {
+ for (size_t j = 0; j <= degree - k; ++j) {
+ for (size_t i = 0; i <= degree - j - k; ++i) {
+ double sum = 0;
+
+ for (size_t k2 = 0; k2 <= j; ++k2) {
+ for (size_t j2 = 0; j2 <= j; ++j2) {
+ for (size_t i2 = 0; i2 <= i; ++i2) {
+ sum += dot(u_bar_coeffs[index(i - i2, j - j2, k - k2)], u_bar_coeffs[index(i2, j2, k2)]);
+ }
+ }
+ }
+
+ kappa_bar_coeffs[index(i, j, k)] = 0.5 * sum;
+ }
+ }
+ }
+
+ double kappa_bar_mean = 0;
+ for (size_t i = 0; i < basis_means.size(); ++i) {
+ kappa_bar_mean += basis_means[i] * kappa_bar_coeffs[i];
+ }
+
+ kappa_bar_coeffs[0] += kappa_j[cell_id] - kappa_bar_mean;
+ });
+}
+
+template <MeshConcept MeshType>
+void
+test_developed_reconstruction(const std::shared_ptr<const MeshType>& p_mesh,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_v)
+{
+ static constexpr size_t Dimension = MeshType::Dimension;
+ using Rd = TinyVector<Dimension>;
+ const size_t degree = 2;
+
+ PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::boundary, degree};
+
+ auto rho_h = rho_v->get<DiscreteFunctionP0<const double>>();
+ auto u_h = u_v->get<DiscreteFunctionP0<const Rd>>();
+ auto E_h = E_v->get<DiscreteFunctionP0<const double>>();
+
+ auto rho_E_h = rho_h * E_h;
+
+ DiscreteFunctionP0<double> kappa_h{p_mesh};
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { kappa_h[cell_id] = 0.5 * dot(u_h[cell_id], u_h[cell_id]); });
+
+ DiscreteFunctionP0<double> epsilon_h{p_mesh};
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { epsilon_h[cell_id] = E_h[cell_id] - kappa_h[cell_id]; });
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(rho_v, u_v, epsilon_h, rho_E_h);
+
+ auto rho_E_bar = reconstructions[3]->get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ auto polynomial_basis_means = PolynomialBasisDataManager::instance()
+ .getPolynomialBasisData(*p_mesh)
+ .getPolynomialBasisMeans(descriptor.integrationMethodType(), descriptor.degree());
+
+ auto remove_mean = [&](const auto& f_h, const auto& g_h, auto fg_bar) {
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+ auto basis_means = polynomial_basis_means[cell_id];
+ using DataType_fg = std::decay_t<decltype(fg_bar_coeffs[0])>;
+
+ DataType_fg fg_bar_mean;
+ if constexpr (std::is_arithmetic_v<DataType_fg>) {
+ fg_bar_mean = 0;
+ } else {
+ fg_bar_mean = zero;
+ }
+
+ for (size_t i = 0; i < basis_means.size(); ++i) {
+ fg_bar_mean += basis_means[i] * fg_bar_coeffs[i];
+ }
+
+ fg_bar_coeffs[0] += f_h[cell_id] * g_h[cell_id] - fg_bar_mean;
+ });
+ };
+
+ auto compute_HP = [&](const auto& f_bar, const auto& fg_bar, const auto& g_h, auto& HP_fg) {
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+ auto f_bar_coeffs = f_bar.coefficients(cell_id);
+
+ auto HP_fg_coeffs = HP_fg.coefficients(cell_id);
+
+ const auto& gj = g_h[cell_id];
+
+ for (size_t i = 0; i < HP_fg_coeffs.size(); ++i) {
+ HP_fg_coeffs[i] = fg_bar_coeffs[i] - f_bar_coeffs[i] * gj;
+ }
+ });
+ };
+
+ auto rho_bar = reconstructions[0]->get<DiscreteFunctionDPk<Dimension, const double>>();
+ auto u_bar = reconstructions[1]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+
+ DiscreteFunctionDPk<Dimension, Rd> rhou_bar{p_mesh, rho_bar.degree()};
+
+ computeFGBar(p_mesh, rho_bar, u_bar, rhou_bar);
+ // compute_fg_bar(rho_bar, u_bar, rhou_bar);
+ remove_mean(rho_h, u_h, rhou_bar);
+
+ DiscreteFunctionDPk<Dimension, Rd> HP_rhou{p_mesh, rho_bar.degree()};
+
+ compute_HP(rho_bar, rhou_bar, u_h, HP_rhou);
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ Rd HP_rhou_mean = zero;
+ auto HP_fg_coeffs = HP_rhou.coefficients(cell_id);
+ auto basis_means = polynomial_basis_means[cell_id];
+
+ for (size_t i = 0; i < basis_means.size(); ++i) {
+ HP_rhou_mean += basis_means[i] * HP_fg_coeffs[i];
+ }
+ });
+
+ auto u_bar_R = [&](const CellId cell_id, const Rd& x) {
+ return u_h[cell_id] + (1. / rho_bar[cell_id](x)) * HP_rhou[cell_id](x);
+ };
+
+ auto epsilon_bar = reconstructions[2]->get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ DiscreteFunctionDPk<Dimension, double> rhoepsilon_bar{p_mesh, rho_bar.degree()};
+
+ computeFGBar(p_mesh, rho_bar, epsilon_bar, rhoepsilon_bar);
+ remove_mean(rho_h, epsilon_h, rhoepsilon_bar);
+
+ DiscreteFunctionDPk<Dimension, double> HP_rhoepsilon{p_mesh, rho_bar.degree()};
+
+ compute_HP(rho_bar, rhoepsilon_bar, epsilon_h, HP_rhoepsilon);
+
+ auto epsilon_bar_R = [&](const CellId cell_id, const Rd& x) {
+ return epsilon_h[cell_id] + (1. / rho_bar[cell_id](x)) * HP_rhoepsilon[cell_id](x);
+ };
+
+ DiscreteFunctionDPk<Dimension, const double> kappa_bar = [&] {
+ DiscreteFunctionDPk<Dimension, double> tmp_kappa_bar{p_mesh, rho_bar.degree()};
+ computeKappaBar(descriptor, p_mesh, u_bar, kappa_h, tmp_kappa_bar);
+ return tmp_kappa_bar;
+ }();
+
+ DiscreteFunctionDPk<Dimension, double> rhokappa_bar{p_mesh, rho_bar.degree()};
+
+ computeFGBar(p_mesh, rho_bar, kappa_bar, rhokappa_bar);
+ remove_mean(rho_h, kappa_h, rhokappa_bar);
+
+ DiscreteFunctionDPk<Dimension, double> HP_rhokappa{p_mesh, rho_bar.degree()};
+
+ compute_HP(rho_bar, rhokappa_bar, kappa_h, HP_rhokappa);
+
+ auto kappa_bar_R = [&](const CellId cell_id, const Rd& x) {
+ return kappa_h[cell_id] + (1. / rho_bar[cell_id](x)) * HP_rhokappa[cell_id](x);
+ };
+
+ auto epsilon_smile_R = [&](const CellId cell_id, const Rd& x) {
+ auto u_bar_R_x = u_bar_R(cell_id, x);
+ return epsilon_bar_R(cell_id, x) + kappa_bar_R(cell_id, x) - 0.5 * dot(u_bar_R_x, u_bar_R_x);
+ };
+
+ if constexpr (Dimension == 1) {
+ std::cout << "creating output...";
+ {
+ std::ofstream fout("comp-u");
+ fout.precision(15);
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ auto x0 = xr[cell_to_node_matrix[cell_id][0]];
+ auto x1 = xr[cell_to_node_matrix[cell_id][1]];
+ for (size_t i = 0; i <= 100; ++i) {
+ auto x = (1 - i / 100.) * x0 + i / 100. * x1;
+ fout << x[0] << ' ' << u_bar[cell_id](x)[0] << ' ' << u_bar_R(cell_id, x)[0] << '\n';
+ }
+ fout << "\n\n";
+ }
+ }
+ {
+ std::ofstream fout("comp-epsilon");
+ fout.precision(15);
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ auto x0 = xr[cell_to_node_matrix[cell_id][0]];
+ auto x1 = xr[cell_to_node_matrix[cell_id][1]];
+ for (size_t i = 0; i <= 100; ++i) {
+ auto x = (1 - i / 100.) * x0 + i / 100. * x1;
+ fout << x[0] << ' ' << epsilon_bar[cell_id](x) << ' ' << epsilon_bar_R(cell_id, x) << ' '
+ << epsilon_smile_R(cell_id, x) << '\n';
+ }
+ fout << "\n\n";
+ }
+ }
+ {
+ std::ofstream fout("comp-epsilon");
+ fout.precision(15);
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ auto x0 = xr[cell_to_node_matrix[cell_id][0]];
+ auto x1 = xr[cell_to_node_matrix[cell_id][1]];
+ for (size_t i = 0; i <= 100; ++i) {
+ auto x = (1 - i / 100.) * x0 + i / 100. * x1;
+ fout << x[0] << ' ' << epsilon_bar[cell_id](x) << ' ' << epsilon_bar_R(cell_id, x) << ' '
+ << epsilon_smile_R(cell_id, x) << '\n';
+ }
+ fout << "\n\n";
+ }
+ }
+
+ {
+ std::ofstream fout("comp-rho_E");
+ fout.precision(15);
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ auto x0 = xr[cell_to_node_matrix[cell_id][0]];
+ auto x1 = xr[cell_to_node_matrix[cell_id][1]];
+ for (size_t i = 0; i <= 100; ++i) {
+ auto x = (1 - i / 100.) * x0 + i / 100. * x1;
+ fout << x[0] << ' ' << rho_E_bar[cell_id](x) << ' '
+ << rho_bar[cell_id](x) *
+ (epsilon_smile_R(cell_id, x) + 0.5 * dot(u_bar_R(cell_id, x), u_bar_R(cell_id, x)))
+ << '\n';
+ }
+ fout << "\n\n";
+ }
+ }
+
+ std::cout << "done\n";
+
+ if constexpr (Dimension == 1) {
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+ double sum = 0;
+ auto qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{10});
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ LineTransformation<Dimension> t{xr[cell_to_node_matrix[cell_id][0]], xr[cell_to_node_matrix[cell_id][1]]};
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ auto x = t(qf.point(i_q));
+ sum +=
+ std::abs(qf.weight(i_q) * (rho_E_bar[cell_id](x) -
+ rho_bar[cell_id](x) * (epsilon_smile_R(cell_id, x) +
+ 0.5 * dot(u_bar_R(cell_id, x), u_bar_R(cell_id, x)))));
+ }
+ }
+
+ std::cout << "L1 error = " << sum << '\n';
+ }
+ }
+}
+
+template <MeshConcept MeshType>
+void
+test_developed_reconstruction2(const std::shared_ptr<const MeshType>& p_mesh,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_v)
+{
+ static constexpr size_t Dimension = MeshType::Dimension;
+ using Rd = TinyVector<Dimension>;
+ const size_t degree = 3;
+
+ PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::boundary, degree};
+
+ auto rho_h = rho_v->get<DiscreteFunctionP0<const double>>();
+ auto u_h = u_v->get<DiscreteFunctionP0<const Rd>>();
+ auto E_h = E_v->get<DiscreteFunctionP0<const double>>();
+
+ auto rho_E_h = rho_h * E_h;
+
+ DiscreteFunctionP0<double> kappa_h{p_mesh};
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { kappa_h[cell_id] = 0.5 * dot(u_h[cell_id], u_h[cell_id]); });
+
+ DiscreteFunctionP0<double> epsilon_h{p_mesh};
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { epsilon_h[cell_id] = E_h[cell_id] - kappa_h[cell_id]; });
+
+ auto rho_epsilon_h = rho_h * epsilon_h;
+ auto rho_kappa_h = rho_h * kappa_h;
+ auto rho_u_h = rho_h * u_h;
+
+ auto reconstructions =
+ PolynomialReconstruction{descriptor}.build(rho_v, rho_u_h, rho_epsilon_h, rho_kappa_h, rho_E_h, epsilon_h, u_v);
+
+ auto rho_bar = reconstructions[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+ auto rho_u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ auto rho_epsilon_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+ auto rho_kappa_bar = reconstructions[3]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+ auto rho_E_bar = reconstructions[4]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+ auto epsilon_bar = reconstructions[5]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+ auto u_bar = reconstructions[6]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+
+ auto polynomial_basis_means = PolynomialBasisDataManager::instance()
+ .getPolynomialBasisData(*p_mesh)
+ .getPolynomialBasisMeans(descriptor.integrationMethodType(), descriptor.degree());
+
+ auto remove_mean = [&](const auto& f_h, const auto& g_h, auto fg_bar) {
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+ auto basis_means = polynomial_basis_means[cell_id];
+ using DataType_fg = std::decay_t<decltype(fg_bar_coeffs[0])>;
+
+ DataType_fg fg_bar_mean;
+ if constexpr (std::is_arithmetic_v<DataType_fg>) {
+ fg_bar_mean = 0;
+ } else {
+ fg_bar_mean = zero;
+ }
+
+ for (size_t i = 0; i < basis_means.size(); ++i) {
+ fg_bar_mean += basis_means[i] * fg_bar_coeffs[i];
+ }
+
+ fg_bar_coeffs[0] += f_h[cell_id] * g_h[cell_id] - fg_bar_mean;
+ });
+ };
+
+ auto compute_HP = [&](const auto& f_bar, const auto& fg_bar, const auto& g_h, auto& HP_fg) {
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+ auto f_bar_coeffs = f_bar.coefficients(cell_id);
+
+ auto HP_fg_coeffs = HP_fg.coefficients(cell_id);
+
+ const auto& gj = g_h[cell_id];
+
+ for (size_t i = 0; i < HP_fg_coeffs.size(); ++i) {
+ HP_fg_coeffs[i] = fg_bar_coeffs[i] - f_bar_coeffs[i] * gj;
+ }
+ });
+ };
+
+ // auto rho_bar = reconstructions[0]->get<DiscreteFunctionDPk<Dimension, const double>>();
+ // auto u_bar = reconstructions[1]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
+
+ // DiscreteFunctionDPk<Dimension, Rd> rhou_bar{p_mesh, rho_bar.degree()};
+
+ // computeFGBar(p_mesh, rho_bar, u_bar, rhou_bar);
+ // compute_fg_bar(rho_bar, u_bar, rhou_bar);
+ // remove_mean(rho_h, u_h, rhou_bar);
+
+ DiscreteFunctionDPk<Dimension, Rd> HP_rhou{p_mesh, rho_bar.degree()};
+
+ compute_HP(rho_bar, rho_u_bar, u_h, HP_rhou);
+
+ // parallel_for(
+ // p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ // Rd HP_rhou_mean = zero;
+ // auto HP_fg_coeffs = HP_rhou.coefficients(cell_id);
+ // auto basis_means = polynomial_basis_means[cell_id];
+
+ // for (size_t i = 0; i < basis_means.size(); ++i) {
+ // HP_rhou_mean += basis_means[i] * HP_fg_coeffs[i];
+ // }
+ // });
+
+ auto u_bar_R = [&](const CellId cell_id, const Rd& x) {
+ return u_h[cell_id] + (1. / rho_bar[cell_id](x)) * HP_rhou[cell_id](x);
+ };
+
+ // auto epsilon_bar = reconstructions[2]->get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ // DiscreteFunctionDPk<Dimension, double> rhoepsilon_bar{p_mesh, rho_bar.degree()};
+
+ // computeFGBar(p_mesh, rho_bar, epsilon_bar, rhoepsilon_bar);
+ // remove_mean(rho_h, epsilon_h, rhoepsilon_bar);
+
+ DiscreteFunctionDPk<Dimension, double> HP_rhoepsilon{p_mesh, rho_bar.degree()};
+
+ compute_HP(rho_bar, rho_epsilon_bar, epsilon_h, HP_rhoepsilon);
+
+ auto epsilon_bar_R = [&](const CellId cell_id, const Rd& x) {
+ return epsilon_h[cell_id] + (1. / rho_bar[cell_id](x)) * HP_rhoepsilon[cell_id](x);
+ };
+
+ // DiscreteFunctionDPk<Dimension, const double> kappa_bar = [&] {
+ // DiscreteFunctionDPk<Dimension, double> tmp_kappa_bar{p_mesh, rho_bar.degree()};
+ // computeKappaBar(descriptor, p_mesh, u_bar, kappa_h, tmp_kappa_bar);
+ // return tmp_kappa_bar;
+ // }();
+
+ // DiscreteFunctionDPk<Dimension, double> rhokappa_bar{p_mesh, rho_bar.degree()};
+
+ // computeFGBar(p_mesh, rho_bar, kappa_bar, rhokappa_bar);
+ // remove_mean(rho_h, kappa_h, rhokappa_bar);
+
+ DiscreteFunctionDPk<Dimension, double> HP_rhokappa{p_mesh, rho_bar.degree()};
+
+ compute_HP(rho_bar, rho_kappa_bar, kappa_h, HP_rhokappa);
+
+ auto kappa_bar_R = [&](const CellId cell_id, const Rd& x) {
+ return kappa_h[cell_id] + (1. / rho_bar[cell_id](x)) * HP_rhokappa[cell_id](x);
+ };
+
+ auto epsilon_smile_R = [&](const CellId cell_id, const Rd& x) {
+ auto u_bar_R_x = u_bar_R(cell_id, x);
+ return epsilon_bar_R(cell_id, x) + kappa_bar_R(cell_id, x) - 0.5 * dot(u_bar_R_x, u_bar_R_x);
+ };
+
+ if constexpr (Dimension == 1) {
+ std::cout << "creating output...";
+ {
+ std::ofstream fout("comp-u");
+ fout.precision(15);
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ auto x0 = xr[cell_to_node_matrix[cell_id][0]];
+ auto x1 = xr[cell_to_node_matrix[cell_id][1]];
+ for (size_t i = 0; i <= 100; ++i) {
+ auto x = (1 - i / 100.) * x0 + i / 100. * x1;
+ fout << x[0] << ' ' << u_bar[cell_id](x)[0] << ' ' << u_bar_R(cell_id, x)[0] << '\n';
+ }
+ fout << "\n\n";
+ }
+ }
+ {
+ std::ofstream fout("comp-epsilon");
+ fout.precision(15);
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ auto x0 = xr[cell_to_node_matrix[cell_id][0]];
+ auto x1 = xr[cell_to_node_matrix[cell_id][1]];
+ for (size_t i = 0; i <= 100; ++i) {
+ auto x = (1 - i / 100.) * x0 + i / 100. * x1;
+ fout << x[0] << ' ' << epsilon_bar[cell_id](x) << ' ' << epsilon_bar_R(cell_id, x) << ' '
+ << epsilon_smile_R(cell_id, x) << '\n';
+ }
+ fout << "\n\n";
+ }
+ }
+ {
+ std::ofstream fout("comp-epsilon");
+ fout.precision(15);
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ auto x0 = xr[cell_to_node_matrix[cell_id][0]];
+ auto x1 = xr[cell_to_node_matrix[cell_id][1]];
+ for (size_t i = 0; i <= 100; ++i) {
+ auto x = (1 - i / 100.) * x0 + i / 100. * x1;
+ fout << x[0] << ' ' << epsilon_bar[cell_id](x) << ' ' << epsilon_bar_R(cell_id, x) << ' '
+ << epsilon_smile_R(cell_id, x) << '\n';
+ }
+ fout << "\n\n";
+ }
+ }
+
+ {
+ std::ofstream fout("comp-rho_E");
+ fout.precision(15);
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ auto x0 = xr[cell_to_node_matrix[cell_id][0]];
+ auto x1 = xr[cell_to_node_matrix[cell_id][1]];
+ for (size_t i = 0; i <= 100; ++i) {
+ auto x = (1 - i / 100.) * x0 + i / 100. * x1;
+ fout << x[0] << ' ' << rho_E_bar[cell_id](x) << ' '
+ << rho_bar[cell_id](x) *
+ (epsilon_smile_R(cell_id, x) + 0.5 * dot(u_bar_R(cell_id, x), u_bar_R(cell_id, x)))
+ << '\n';
+ }
+ fout << "\n\n";
+ }
+ }
+
+ std::cout << "done\n";
+
+ if constexpr (Dimension == 1) {
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+ double sum = 0;
+ auto qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{10});
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ LineTransformation<Dimension> t{xr[cell_to_node_matrix[cell_id][0]], xr[cell_to_node_matrix[cell_id][1]]};
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ auto x = t(qf.point(i_q));
+ sum +=
+ std::abs(qf.weight(i_q) * (rho_E_bar[cell_id](x) -
+ rho_bar[cell_id](x) * (epsilon_smile_R(cell_id, x) +
+ 0.5 * dot(u_bar_R(cell_id, x), u_bar_R(cell_id, x)))));
+ }
+ }
+
+ std::cout << "L1 error = " << sum << '\n';
+ }
+ }
+}
+
+template <MeshConcept MeshType>
+void
+test_developed_reconstruction3(const std::shared_ptr<const MeshType>& p_mesh,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_v)
+{
+ static constexpr size_t Dimension = MeshType::Dimension;
+ using Rd = TinyVector<Dimension>;
+ const size_t degree = 3;
+
+ PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::boundary, degree};
+
+ auto rho_h = rho_v->get<DiscreteFunctionP0<const double>>();
+ auto u_h = u_v->get<DiscreteFunctionP0<const Rd>>();
+ auto E_h = E_v->get<DiscreteFunctionP0<const double>>();
+
+ auto rho_u_h = rho_h * u_h;
+ auto rho_E_h = rho_h * E_h;
+
+ auto reconstructions = PolynomialReconstruction{descriptor}.build(rho_v, u_v, rho_u_h);
+
+ auto rho_u_bar_cons = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+
+ auto polynomial_basis_means = PolynomialBasisDataManager::instance()
+ .getPolynomialBasisData(*p_mesh)
+ .getPolynomialBasisMeans(descriptor.integrationMethodType(), descriptor.degree());
+
+ auto remove_mean = [&](const auto& f_h, const auto& g_h, auto fg_bar) {
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+ auto basis_means = polynomial_basis_means[cell_id];
+ using DataType_fg = std::decay_t<decltype(fg_bar_coeffs[0])>;
+
+ DataType_fg fg_bar_mean;
+ if constexpr (std::is_arithmetic_v<DataType_fg>) {
+ fg_bar_mean = 0;
+ } else {
+ fg_bar_mean = zero;
+ }
+
+ for (size_t i = 0; i < basis_means.size(); ++i) {
+ fg_bar_mean += basis_means[i] * fg_bar_coeffs[i];
+ }
+
+ fg_bar_coeffs[0] += f_h[cell_id] * g_h[cell_id] - fg_bar_mean;
+ });
+ };
+
+ auto compute_HP = [&](const auto& f_bar, const auto& fg_bar, const auto& g_h, auto& HP_fg) {
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+ auto f_bar_coeffs = f_bar.coefficients(cell_id);
+
+ auto HP_fg_coeffs = HP_fg.coefficients(cell_id);
+
+ const auto& gj = g_h[cell_id];
+
+ for (size_t i = 0; i < HP_fg_coeffs.size(); ++i) {
+ HP_fg_coeffs[i] = fg_bar_coeffs[i] - f_bar_coeffs[i] * gj;
+ }
+ });
+ };
+
+ auto rho_bar_cons = reconstructions[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+ auto u_bar_cons = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+
+ DiscreteFunctionDPk<Dimension, const double> rho_bar = [&] {
+ auto v = copy(rho_bar_cons);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto v_coeffs = v.coefficients(cell_id);
+ v_coeffs[0] = rho_h[cell_id];
+ });
+ return v;
+ }();
+
+ DiscreteFunctionDPk<Dimension, const Rd> u_bar = [&] {
+ auto v = copy(u_bar_cons);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ auto u_bar_coeffs = v.coefficients(cell_id);
+ u_bar_coeffs[0] = u_h[cell_id];
+ });
+ return v;
+ }();
+
+ DiscreteFunctionDPk<Dimension, Rd> rhou_bar{p_mesh, rho_bar.degree()};
+
+ computeFGBar(p_mesh, rho_bar, u_bar, rhou_bar);
+ remove_mean(rho_h, u_h, rhou_bar);
+
+ DiscreteFunctionDPk<Dimension, Rd> HP_rhou{p_mesh, rho_bar.degree()};
+
+ compute_HP(rho_bar_cons, rhou_bar, u_h, HP_rhou);
+
+ auto u_bar_R = [&](const CellId cell_id, const Rd& x) {
+ return u_h[cell_id] + (1. / rho_bar_cons[cell_id](x)) * HP_rhou[cell_id](x);
+ };
+
+ if constexpr (Dimension == 1) {
+ std::cout << "done\n";
+
+ if constexpr (Dimension == 1) {
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+ {
+ double sum = 0;
+ auto qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{10});
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ LineTransformation<Dimension> t{xr[cell_to_node_matrix[cell_id][0]], xr[cell_to_node_matrix[cell_id][1]]};
+ for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ auto x = t(qf.point(i_q));
+ auto diff = (rho_u_bar_cons[cell_id](x) - rhou_bar[cell_id](x));
+ sum += std::abs(qf.weight(i_q) * std::abs(diff[0]));
+ }
+ }
+ std::cout << "rho u: L1 error = " << sum << '\n';
+ }
+
+ // {
+ // double sum = 0;
+ // auto qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{10});
+ // for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ // LineTransformation<Dimension> t{xr[cell_to_node_matrix[cell_id][0]], xr[cell_to_node_matrix[cell_id][1]]};
+ // for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+ // auto x = t(qf.point(i_q));
+ // auto diff = (rho_u_bar_cons[cell_id](x) - rho_bar_cons[cell_id](x) * (u_bar_R(cell_id, x)));
+ // sum += std::abs(qf.weight(i_q) * std::abs(diff[0]));
+ // }
+ // }
+ // std::cout << "rho u: L1 error = " << sum << '\n';
+ // }
+ }
+ }
+}
+
+void
+test_developed_reconstruction(const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_v)
+{
+ auto mesh_v = getCommonMesh({rho_v, u_v, E_v});
+
+ std::visit([&](auto&& p_mesh) { test_developed_reconstruction3(p_mesh, rho_v, u_v, E_v); }, mesh_v->variant());
+}
diff --git a/src/scheme/test_developed_reconstruction.hpp b/src/scheme/test_developed_reconstruction.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..625506fbe605702c4cb200650b341f3ec6d5f62f
--- /dev/null
+++ b/src/scheme/test_developed_reconstruction.hpp
@@ -0,0 +1,10 @@
+#ifndef TEST_DEVELOPED_RECONSTRUCTION_HPP
+#define TEST_DEVELOPED_RECONSTRUCTION_HPP
+
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+void test_developed_reconstruction(const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_v);
+
+#endif // TEST_DEVELOPED_RECONSTRUCTION_HPP