diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 5101d5da271894ed113fc0f57f5b19effe5f29ac..77fed256505a4862382e9e05691433a40ec3fa7e 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -977,9 +977,10 @@ PolynomialReconstruction::_build(
                       B(index, kB) = qi_qj[k];
                     }
                   } else if constexpr (is_tiny_matrix_v<DataType>) {
-                    for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
-                      for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
-                        B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+                    for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                      const size_t kB = column_begin + p * DataType::NumberOfColumns;
+                      for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                        B(index, kB + q) = qi_qj(p, q);
                       }
                     }
                   }
@@ -1018,9 +1019,9 @@ PolynomialReconstruction::_build(
 
                         const DataType& qi    = discrete_function[cell_i_id];
                         const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
-                        for (size_t k = 0; k < DataType::NumberOfRows; ++k) {
-                          for (size_t l = 0; l < DataType::NumberOfColumns; ++l) {
-                            B(index, column_begin + k * DataType::NumberOfColumns + l) = qi_qj(k, l);
+                        for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                          for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                            B(index, column_begin + p * DataType::NumberOfColumns + q) = qi_qj(p, q);
                           }
                         }
                       } else {
@@ -1114,12 +1115,55 @@ PolynomialReconstruction::_build(
                     }
                   }
                 } else if constexpr (is_tiny_matrix_v<DataType>) {
+                  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];
+                    for (size_t l = 0; l < qj_vector.size(); ++l) {
+                      const DataType& qj    = qj_vector[l];
+                      const DataType& qi    = discrete_function[cell_i_id][l];
+                      const DataType& qi_qj = qi - qj;
+
+                      for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                        const size_t kB = column_begin + l * DataType::Dimension + p * DataType::NumberOfColumns;
+                        for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                          B(index, kB + q) = qi_qj(p, q);
+                        }
+                      }
+                    }
+                  }
+
                   for (size_t i_symmetry = 0; i_symmetry < stencil_array.symmetryBoundaryStencilArrayList().size();
                        ++i_symmetry) {
-                    auto& ghost_stencil  = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
-                    auto ghost_cell_list = ghost_stencil[cell_j_id];
-                    if (ghost_cell_list.size() > 0) {
-                      throw NotImplementedError("TinyMatrix symmetries for reconstruction of DiscreteFunctionP0Vector");
+                    if constexpr ((DataType::NumberOfRows == MeshType::Dimension) and
+                                  (DataType::NumberOfColumns == MeshType::Dimension)) {
+                      auto& ghost_stencil = stencil_array.symmetryBoundaryStencilArrayList()[i_symmetry].stencilArray();
+                      auto ghost_cell_list = ghost_stencil[cell_j_id];
+
+                      const Rd& normal = symmetry_normal_list[i_symmetry];
+
+                      for (size_t i = 0; i < ghost_cell_list.size(); ++i, ++index) {
+                        const CellId cell_i_id = ghost_cell_list[i];
+
+                        for (size_t l = 0; l < qj_vector.size(); ++l) {
+                          const DataType& qj    = qj_vector[l];
+                          const DataType& qi    = discrete_function[cell_i_id][l];
+                          const DataType& qi_qj = symmetrize_matrix(normal, qi) - qj;
+
+                          for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                            const size_t kB = column_begin + l * DataType::Dimension + p * DataType::NumberOfColumns;
+                            for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                              B(index, kB + q) = qi_qj(p, q);
+                            }
+                          }
+                        }
+                      }
+                    } else {
+                      // LCOV_EXCL_START
+                      std::stringstream error_msg;
+                      error_msg << "cannot symmetrize vectors of dimension " << DataType::Dimension
+                                << " using a mesh of dimension " << MeshType::Dimension;
+                      throw UnexpectedError(error_msg.str());
+                      // LCOV_EXCL_STOP
                     }
                   }
                 }
@@ -1466,6 +1510,29 @@ PolynomialReconstruction::_build(
                           }
                         }
                         column_begin += DataType::Dimension;
+                      } else if constexpr (is_tiny_matrix_v<DataType>) {
+                        if (m_descriptor.degree() > 1) {
+                          auto& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+                          for (size_t i = 0; i < basis_dimension - 1; ++i) {
+                            auto& dpk_j_0 = dpk_j[component_begin];
+                            for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                              for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                                dpk_j_0(p, q) -=
+                                  X(i, column_begin + p * DataType::NumberOfColumns + q) * mean_j_of_ejk[i];
+                              }
+                            }
+                          }
+                        }
+
+                        for (size_t i = 0; i < basis_dimension - 1; ++i) {
+                          auto& dpk_j_ip1 = dpk_j[component_begin + i + 1];
+                          for (size_t p = 0; p < DataType::NumberOfRows; ++p) {
+                            for (size_t q = 0; q < DataType::NumberOfColumns; ++q) {
+                              dpk_j_ip1(p, q) = X(i, column_begin + p * DataType::NumberOfColumns + q);
+                            }
+                          }
+                        }
+                        column_begin += DataType::Dimension;
                       } else {
                         // LCOV_EXCL_START
                         throw UnexpectedError("unexpected data type");