diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index a6e35649ef1401bd3eda2b885d6f199be473b7e4..b01d11a792bcd250857758c52d1c4885b78d896f 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -247,6 +247,7 @@ _computeEjkMean(const QuadratureFormula<2>& quadrature,
       }
 
       for (size_t i_y = 1; i_y <= degree; ++i_y) {
+#warning store y_row index and size
         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) {
           inv_Vj_wq_detJ_ek[k++] = y_yj * inv_Vj_wq_detJ_ek[begin_i_y_1 + l];
@@ -271,6 +272,29 @@ _computeEjkMean(const QuadratureFormula<3>& quadrature,
                 SmallArray<double>& inv_Vj_wq_detJ_ek,
                 SmallArray<double>& mean_of_ejk)
 {
+#warning Compute this one for all and rework design
+  SmallArray<size_t> m_yz_row_index((degree + 2) * (degree + 1) / 2 + 1);
+  SmallArray<size_t> m_z_triangle_index(degree + 1);
+
+  {
+    size_t i_z  = 0;
+    size_t i_yz = 0;
+
+    m_yz_row_index[i_yz++] = 0;
+    for (ssize_t n = degree + 1; n >= 1; --n) {
+      m_z_triangle_index[i_z++] = i_yz - 1;
+      for (ssize_t i = n; i >= 1; --i) {
+        m_yz_row_index[i_yz] = m_yz_row_index[i_yz - 1] + i;
+        ++i_yz;
+      }
+    }
+  }
+
+  SmallArray<size_t> m_yz_row_size{m_yz_row_index.size() - 1};
+  for (size_t i = 0; i < m_yz_row_size.size(); ++i) {
+    m_yz_row_size[i] = m_yz_row_index[i + 1] - m_yz_row_index[i];
+  }
+
   using Rd = TinyVector<3>;
   mean_of_ejk.fill(0);
 
@@ -291,68 +315,33 @@ _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 = 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";
+        const size_t begin_i_y_1 = m_yz_row_index[i_y - 1];
+        const size_t nb_monoms   = m_yz_row_size[i_y];
+        for (size_t l = 0; l < nb_monoms; ++l) {
           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 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_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) {
-            std::clog << "[" << k << "] <- z*[" << begin_i_yz_1 + l << "]\n";
+        const size_t nb_y      = m_yz_row_size[m_z_triangle_index[i_z]];
+        const size_t index_z   = m_z_triangle_index[i_z];
+        const size_t index_z_1 = m_z_triangle_index[i_z - 1];
+        for (size_t i_y = 0; i_y < nb_y; ++i_y) {
+          const size_t begin_i_yz_1 = m_yz_row_index[index_z_1 + i_y];
+          const size_t nb_monoms    = m_yz_row_size[index_z + i_y];
+          for (size_t l = 0; l < nb_monoms; ++l) {
             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) {
@@ -1316,29 +1305,6 @@ 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];
@@ -1346,21 +1312,11 @@ 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];
               }
             }
+
             for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
                  ++i_symmetry) {
               auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
@@ -1456,22 +1412,6 @@ PolynomialReconstruction::_build(
             }
           }
         } else {
-#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';
-            //     ;
-            //   }
-            // }
-          }
           Givens::solveCollectionInPlace(A, X, B);
         }
 
diff --git a/tests/test_PolynomialReconstruction_degree_3.cpp b/tests/test_PolynomialReconstruction_degree_3.cpp
index 4f50a7e46ba2729ff379fcd34879e7ad98fe7724..d8badda38a09584a34696b7311120f7b859e4b19 100644
--- a/tests/test_PolynomialReconstruction_degree_3.cpp
+++ b/tests/test_PolynomialReconstruction_degree_3.cpp
@@ -191,15 +191,9 @@ TEST_CASE("PolynomialReconstruction_degree_3", "[scheme]")
                 auto& mesh  = *p_mesh;
 
                 std::array<std::function<R3(const R1&)>, 3> vector_exact   //
-                  = {[&](const R1& x) -> R3 {
-                       return R3{p1(x), p2(x), p3(x)};
-                     },
-                     [&](const R1& x) -> R3 {
-                       return R3{p5(x), p7(x), p0(x)};
-                     },
-                     [&](const R1& x) -> R3 {
-                       return R3{p9(x), p8(x), p4(x)};
-                     }};
+                  = {[&](const R1& x) -> R3 { return R3{p1(x), p2(x), p3(x)}; },
+                     [&](const R1& x) -> R3 { return R3{p5(x), p7(x), p0(x)}; },
+                     [&](const R1& x) -> R3 { return R3{p9(x), p8(x), p4(x)}; }};
 
                 DiscreteFunctionP0Vector Vh = test_only::exact_projection(mesh, degree, vector_exact);
 
@@ -464,15 +458,14 @@ TEST_CASE("PolynomialReconstruction_degree_3", "[scheme]")
             const double y3  = y2 * y;
             const double z3  = z2 * z;
 
-            return yz2;
-
-            // return 2.3 + 1.7 * x - 1.3 * y + 2.1 * z                             //
-            //        + 1.7 * x2 + 1.4 * y2 + 1.7 * z2                              //
-            //        - 2.3 * xy + 1.6 * xz - 1.9 * yz                              //
-            //        + 1.2 * x3 - 2.1 * y3 - 1.1 * z3                              //
-            //        - 1.7 * x2y - 1.3 * x2z + 0.9 * xy2 - 0.7 * y2z + 1.5 * xz2   //- 0.7 * yz2
-            //        + 2.8 * xyz                                                   //
-            //   ;
+            return 2.3 + 1.7 * x - 1.3 * y + 2.1 * z     //
+                   + 1.7 * x2 + 1.4 * y2 + 1.7 * z2      //
+                   - 2.3 * xy + 1.6 * xz - 1.9 * yz      //
+                   + 1.2 * x3 - 2.1 * y3 - 1.1 * z3      //
+                   - 1.7 * x2y - 1.3 * x2z + 0.9 * xy2   //
+                   - 0.7 * y2z + 1.5 * xz2 - 0.7 * yz2   //
+                   + 2.8 * xyz                           //
+              ;
           };
 
           // auto p1 = [](const R3& x) {
@@ -542,48 +535,12 @@ TEST_CASE("PolynomialReconstruction_degree_3", "[scheme]")
                 auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
                 auto dpk_fh          = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
 
-                {
-                  // PolynomialReconstructionDescriptor descriptor2{IntegrationMethodType::element, 2};
-                  // DiscreteFunctionP0 f2h = test_only::exact_projection(mesh, 2, std::function(R_exact));
-                  // auto reconstructions2  = PolynomialReconstruction{descriptor2}.build(fh);
-                  // auto dpk_fh2           = reconstructions2[0]->get<DiscreteFunctionDPk<3, const double>>();
-                  // double max_error       = test_only::max_reconstruction_error(mesh, dpk_fh2,
-                  // std::function(R_exact)); REQUIRE(parallel::allReduceMax(max_error) ==
-                  // Catch::Approx(0).margin(1E-12));
-
-                  auto cell_type   = mesh.connectivity().cellType();
-                  auto cell_number = mesh.connectivity().cellNumber();
-
-                  for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
-                    std::clog << cell_id << "(" << cell_number[cell_id] << ")"
-                              << ":---- " << name(cell_type[cell_id]) << "\n2| ";
-                    // {
-                    //   auto coeffs = dpk_fh2.coefficients(cell_id);
-                    //   for (size_t i = 0; i < coeffs.size(); ++i) {
-                    //     if (std::abs(coeffs[i]) > 1e-12) {
-                    //       std::clog << ' ' << i << ':' << coeffs[i];
-                    //     }
-                    //   }
-                    // }
-                    std::clog << '\n';
-                    std::clog << "3| ";
-                    {
-                      auto coeffs = dpk_fh.coefficients(cell_id);
-                      for (size_t i = 0; i < coeffs.size(); ++i) {
-                        if (std::abs(coeffs[i]) > 1e-12) {
-                          std::clog << ' ' << i << ':' << coeffs[i];
-                        }
-                      }
-                    }
-                    std::clog << '\n';
-                    break;
-                  }
-                }
-
                 double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact));
                 REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
               }
             }
+
+#warning NOT FINISHED
           }
 
           // SECTION("R^3 data")