diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 386d741d6f1b73028e50e2d71d87a22cb4183764..a6e35649ef1401bd3eda2b885d6f199be473b7e4 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -291,29 +291,68 @@ _computeEjkMean(const QuadratureFormula<3>& quadrature,
const double y_yj = X_Xj[1];
const double z_zj = X_Xj[2];
+ size_t n = degree;
+ std::clog << "DEGREE = " << degree << '\n';
+ Array<size_t> yz_begining((n + 2) * (n + 1) / 2 + 1);
+ Array<size_t> z_index(n + 1);
+
+ {
+ size_t i_z = 0;
+ size_t i_yz = 0;
+
+ yz_begining[i_yz++] = 0;
+ for (ssize_t m = n + 1; m >= 1; --m) {
+ z_index[i_z++] = i_yz - 1;
+ for (ssize_t i = m; i >= 1; --i) {
+ yz_begining[i_yz] = yz_begining[i_yz - 1] + i;
+ ++i_yz;
+ }
+ }
+ }
+
+ Array<size_t> yz_nb_monomes{yz_begining.size() - 1};
+ for (size_t i = 0; i < yz_nb_monomes.size(); ++i) {
+ yz_nb_monomes[i] = yz_begining[i + 1] - yz_begining[i];
+ }
+
{
size_t k = 0;
inv_Vj_wq_detJ_ek[k++] = wq * detT / Vi;
for (; k <= degree; ++k) {
inv_Vj_wq_detJ_ek[k] = x_xj * inv_Vj_wq_detJ_ek[k - 1];
+ std::clog << "[" << k << "] <- x*[" << k - 1 << "]\n";
}
for (size_t i_y = 1; i_y <= degree; ++i_y) {
- const size_t begin_i_y_1 = ((i_y - 1) * (2 * degree - i_y + 4)) / 2;
- for (size_t l = 0; l <= degree - i_y; ++l) {
+ const size_t begin_i_y_1 = yz_begining[i_y - 1];
+ std::clog << "COMP 1: " << degree - i_y << " vs " << yz_nb_monomes[i_y] - 1 << '\n';
+ for (size_t l = 0; l < yz_nb_monomes[i_y]; ++l) {
+ std::clog << "[" << k << "] <- y*[" << begin_i_y_1 + l << "]\n";
inv_Vj_wq_detJ_ek[k++] = y_yj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
}
}
for (size_t i_z = 1; i_z <= degree; ++i_z) {
- const size_t begin_i_z_1 = ((i_z - 1) * (3 * (degree + 2) * (degree + 3) + (i_z - (3 * degree + 8)) * i_z)) / 6;
+ const size_t index_z = z_index[i_z];
+ const size_t index_z_1 = z_index[i_z - 1];
+ std::clog << "index_z_1 = " << index_z_1 << '\n';
+ std::clog << "COMP 2: " << degree - i_z << " vs " << yz_nb_monomes[index_z] - 1 << '\n';
for (size_t i_y = 0; i_y <= degree - i_z; ++i_y) {
- const size_t begin_i_y_1 = (i_y * (2 * degree - i_y + 3)) / 2 + begin_i_z_1;
+ const size_t begin_i_yz_1 = yz_begining[index_z_1 + i_y];
+
+ std::clog << "COMP 3: " << degree - i_y - i_z << " vs " << yz_nb_monomes[index_z + i_y] - 1 << '\n';
for (size_t l = 0; l <= degree - i_y - i_z; ++l) {
- inv_Vj_wq_detJ_ek[k++] = z_zj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
+ std::clog << "[" << k << "] <- z*[" << begin_i_yz_1 + l << "]\n";
+ inv_Vj_wq_detJ_ek[k++] = z_zj * inv_Vj_wq_detJ_ek[begin_i_yz_1 + l];
}
}
}
+
+ std::clog << "z_begining = " << z_index << '\n';
+ std::clog << "yz_begining = " << yz_begining << '\n';
+ std::clog << "yz_nb_monomes = " << yz_nb_monomes << '\n';
+
+ std::exit(0);
}
for (size_t k = 1; k < dimension; ++k) {
@@ -1277,6 +1316,29 @@ PolynomialReconstruction::_build(
_computeEjkMean(Xj, xr, cell_to_node_matrix[cell_j_id], cell_type[cell_j_id], Vj[cell_j_id],
m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek, mean_j_of_ejk);
+ auto exact_xz2 = [](const Rd& X) {
+ const double x0 = X[0] - 0.5;
+ const double y0 = X[1] - 0.5;
+ const double z0 = X[2] - 0.5;
+ const double x1 = X[0] + 0.5;
+ const double y1 = X[1] + 0.5;
+ const double z1 = X[2] + 0.5;
+ return (0.5 * x1 * x1 - 0.5 * x0 * x0) * (1. / 3 * z1 * z1 * z1 - 1. / 3 * z0 * z0 * z0);
+ };
+ auto exact_yz2 = [](const Rd& X) {
+ const double x0 = X[0] - 0.5;
+ const double y0 = X[1] - 0.5;
+ const double z0 = X[2] - 0.5;
+ const double x1 = X[0] + 0.5;
+ const double y1 = X[1] + 0.5;
+ const double z1 = X[2] + 0.5;
+ return (0.5 * y1 * y1 - 0.5 * y0 * y0) * (1. / 3 * z1 * z1 * z1 - 1. / 3 * z0 * z0 * z0);
+ };
+ if (cell_j_id == 0) {
+ std::clog << "- mean_j_of_ejk[xz2] = " << mean_j_of_ejk[16] << " vs " << exact_xz2(Xj) << '\n';
+ std::clog << "- mean_j_of_ejk[yz2] = " << mean_j_of_ejk[17] << " vs " << exact_yz2(Xj) << '\n';
+ }
+
size_t index = 0;
for (size_t i = 0; i < stencil_cell_list.size(); ++i, ++index) {
const CellId cell_i_id = stencil_cell_list[i];
@@ -1284,6 +1346,17 @@ PolynomialReconstruction::_build(
_computeEjkMean(Xj, xr, cell_to_node_matrix[cell_i_id], cell_type[cell_i_id], Vj[cell_i_id],
m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek, mean_i_of_ejk);
+ const Rd& Xi = xj[cell_i_id];
+
+ if (cell_j_id == 0) {
+ std::clog << "- mean_i_of_ejk[" << rang::fgB::cyan << "xz2" << rang::fg::reset
+ << "] = " << mean_i_of_ejk[16] << " vs " << exact_xz2(Xi)
+ << "\tdelta=" << std::abs(mean_i_of_ejk[16] - exact_xz2(Xi)) << '\n';
+ std::clog << "- mean_i_of_ejk[" << rang::fgB::yellow << "yz2" << rang::fg::reset
+ << "] = " << mean_i_of_ejk[17] << " vs " << exact_yz2(Xi)
+ << "\tdelta=" << std::abs(mean_i_of_ejk[17] - exact_yz2(Xi)) << '\n';
+ }
+
for (size_t l = 0; l < basis_dimension - 1; ++l) {
A(index, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
}
@@ -1383,23 +1456,23 @@ PolynomialReconstruction::_build(
}
}
} else {
- Givens::solveCollectionInPlace(A, X, B);
#warning REMOVE
if (cell_j_id == 0) {
- for (size_t l = 0; l < B.numberOfRows(); ++l) {
- for (size_t i = 0; i < B.numberOfColumns(); ++i) {
- std::clog << "B(" << l << ", " << i << ") =" << B(l, i) << '\n';
- ;
- }
- }
-
- for (size_t l = 0; l < X.numberOfRows(); ++l) {
- for (size_t i = 0; i < X.numberOfColumns(); ++i) {
- std::clog << "X(" << l << ", " << i << ") =" << X(l, i) << '\n';
- ;
- }
- }
+ // for (size_t l = 0; l < B.numberOfRows(); ++l) {
+ // for (size_t i = 0; i < B.numberOfColumns(); ++i) {
+ // std::clog << "B(" << l << ", " << i << ") =" << B(l, i) << '\n';
+ // ;
+ // }
+ // }
+
+ // for (size_t l = 0; l < X.numberOfRows(); ++l) {
+ // for (size_t i = 0; i < X.numberOfColumns(); ++i) {
+ // std::clog << "X(" << l << ", " << i << ") =" << X(l, i) << '\n';
+ // ;
+ // }
+ // }
}
+ Givens::solveCollectionInPlace(A, X, B);
}
column_begin = 0;
diff --git a/tests/test_PolynomialReconstruction_degree_3.cpp b/tests/test_PolynomialReconstruction_degree_3.cpp
index d04ac14b7282e33e9b35ed2dc7fee18ed9aafb54..4f50a7e46ba2729ff379fcd34879e7ad98fe7724 100644
--- a/tests/test_PolynomialReconstruction_degree_3.cpp
+++ b/tests/test_PolynomialReconstruction_degree_3.cpp
@@ -464,7 +464,7 @@ TEST_CASE("PolynomialReconstruction_degree_3", "[scheme]")
const double y3 = y2 * y;
const double z3 = z2 * z;
- return (y - 0.5) * (z - 0.5) * (z - 0.5);
+ return yz2;
// return 2.3 + 1.7 * x - 1.3 * y + 2.1 * z //
// + 1.7 * x2 + 1.4 * y2 + 1.7 * z2 //
@@ -498,10 +498,10 @@ TEST_CASE("PolynomialReconstruction_degree_3", "[scheme]")
for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
SECTION(named_mesh.name())
{
- auto p_initial_mesh =
- CartesianMeshBuilder{TinyVector<3>{0, 0, 0}, TinyVector<3>{4, 4, 4}, TinyVector<3, size_t>{4, 4, 4}}
- .mesh()
- ->get<Mesh<3>>(); // named_mesh.mesh()->get<Mesh<3>>();
+ auto p_initial_mesh = CartesianMeshBuilder{TinyVector<3>{-0.5, -0.5, -0.5},
+ TinyVector<3>{3.5, 3.5, 3.5}, TinyVector<3, size_t>{4, 4, 4}}
+ .mesh()
+ ->get<Mesh<3>>(); // named_mesh.mesh()->get<Mesh<3>>();
auto& initial_mesh = *p_initial_mesh;
constexpr double theta = 0;